*______________________________________________________________________________
        PROGRAM JQYS_AER_FIT
*______________________________________________________________________________
c
c # Fits a model J(Q,Y) to experimental data.
c # J(Q,Y) obtained as FT of J(Q,s), the intermediate scattering function
c # J(Q,s) = [ (n0+n0*f(s)+n'(s)) ]*R(Q,s) 
c # where  * n0 = Condensate fraction
c          * f(s) = condensate induced peaking function (for low momentum)
c          * n'(s) = broad (uncondensed) component of momentum distribution
c          * R(Q,s) = Final State function
c 
c R.T. AZUAH -- Feb 1996
c          * Adapted for confined geometries by including an extra broad term scattering
c          * J(Q,s) = J(Q,s)_old + Gaussian
c R.T. AZUAH -- May 2001

c
c______________________________________________________________________________
c        implicit none
c        external jqy_cal
c Variables and common blocks
c        character*80 xlab*5,ylab*5,text,xcap*40,ycap*40,title,nam(50)*20
c        character*60 filein,fileres,fun*20,dec*1,subtitle(3)
c        integer*4 nptot,np,inp(50),nd,nrd,i,j,in_len,res_len,mn
c        real*8 p(50),x(1000),y(1000),ysig(1000),rx(1000),ry(1000),er 
c        real*8 mom,wbar,m(7),
c        real*8 qq,Sq,temp
c        logical resfile,moms,comp,firstime,okay

        implicit none
        external jqy_cal 
c
        character xlab*5,ylab*5,xcap*60,ycap*60,subtitle(3)*60,title*80,nam(50)*20
        character*60 filein,fileres,text*80,dec*1
        integer*4 text_len,lpt,strim,nrd
        logical okay,resfile,moms,comp,firstime
        integer*4 nptot,nd,i,j,in_len,res_len,mn
        real*8 p(50),pmin(50),pmax(50),x(1000),y(1000),ysig(1000)
        real*8 rx(1000),ry(1000),er,qq,Sq,temp,mom,wbar,m(7)
c
        common /dats/ y
        common /rdats/ nrd,rx,ry
        common /qval/ qq,Sq,temp
        common /logf/ moms
        common /logc/ firstime,comp,resfile
        common /file/ filein,in_len,fileres,res_len
c

        data nam(1) /'Bk Const'/,
     >       nam(2) /'Bk Slope'/,
     >       nam(3) /'n0   [%]'/,
     >       nam(4) /'kc   [\A\u-1\d]'/,
     >       nam(5) /'a2   [\A\u-2\d]'/,
     >       nam(6) /'b3   [\A\u-3\d]'/,
     >       nam(7) /'a4   [\A\u-4\d]'/,
     >       nam(8) /'b4   [\A\u-4\d]'/,
     >       nam(9) /'b5   [\A\u-5\d]'/,
     >       nam(10) /'a6   [\A\u-6\d]'/,
     >       nam(11) /'b6   [\A\u-6\d]'/,
     >       nam(12) /'Amp [%]'/,
     >       nam(13) /'a2b [\A\u-2\d]'/,
     >       nam(14) /'Weight'/
       data xlab, ylab /'X', 'Y'/
        data subtitle(1) /'FUNCTION TYPE : xxxxx'/
        data subtitle(2) /'Instrument: xxxxxxxx          User: xxxxxxxxxxxxxxxxxxxx'/
        data subtitle(3) /'Q value is: xxxxx     DATE OF EXPT: xxxxxxxxxxxxxxxxxxxx'/
        nptot = 14
 10     write (*,'(T5,A,$)') ' Enter data filename ::'
        read (*,'(A)') filein
        INQUIRE (file=filein,exist=okay)
        IF (okay) THEN
           open (unit=18,file=filein,status='old')
           in_len = strim(filein)
           DO i = 1,1000
              read (18,*,end=15) x(i),y(i),ysig(i)
           END DO
 15        nd = i - 1   
           close (18)
        ELSE
           write (*,*) '    ** ERROR - Data file does not exist '
           GOTO 10
        END IF
        write (*,'(T5,A,$)') ' Enter Q and Temp of data ::'
        read (*,*) qq,temp
 20     write (*,'(T5,A,$)') ' Enter corresponding Instr. Resolution file ::'
        read (*,'(A)') fileres
        INQUIRE (file=fileres,exist=okay)
        IF (okay) THEN
           open (unit=18,file=fileres,status='old')
           res_len = strim(fileres)
           resfile = .true.
           DO i = 1,1000
              read (18,*,end=25) rx(i),ry(i),er
           END DO
 25        nrd = i - 1  
           close (18)
        ELSE
           write (*,*) '    ** ERROR - Instr Resolution file does not exist '
           GOTO 20
        END IF
        write (*,'(T5,A)') ' Decide now whether to keep components of J(Q,Y) separate '
        write (*,'(T5,A,$)') ' Components (slower)? Default [N] ::'
        read (*,'(A)') dec
        IF (dec.EQ.'y'.OR.dec.EQ.'Y')  comp = .true.
c
c Set up labels and captions
c
        write (title,'(A,F4.1,A,F3.1,A)') 'Q=',QQ,'\A\u-1\d  T=',temp,'K'
        xcap = ' Y (\A\u-1\d) '
        ycap = ' J(Q,Y) (\A) '
	IF (comp) THEN
	   subtitle(1)(17:55) = 'CAN PLOT J(Q,Y) COMPONENTS '
	ELSE
	   subtitle(1)(17:55) = 'CANNOT PLOT J(Q,Y) COMPONENTS'
	END IF
        subtitle(2)(13:20) = 'MARI'
        subtitle(2)(37:56) = 'RTA/HRG'
        write (subtitle(3)(13:17), '(F5.2)') QQ
        subtitle(3)(37:56) = 'APRIL 1998'
c Welcome Message
        write (*,*) '#'
        write (*,*) '# ******** JQYS Fitting Program ****** '
        write (*,*) '#'
        write (*,*) '# Type HELP for a list of available commands'
        write (*,*) '#'
        write (*,*) '# Details of fits are stored in /tmp/jqys.lpt '
        write (*,*) '#'
        write (*,*) '#'
        Sq = 1
        call area (x,y,ysig,nd,Sq) 	! normalise exp't data to S(Q)
        p(1) = 0.0			! estimate initial fitting parameters
        p(2) = 0.0
        p(3) = 7.0
        p(4) = 0.5
        p(5) = 0.897
        p(6) = 2.75/QQ
        p(7) = 0.46
        p(8) = 0.0
        p(9) = 2200.0/QQ**3
        p(10) = 0.38
        p(11) = 175.0/QQ**2
        p(12) = 10.0
        p(13) = 2*p(5)
        p(14) = Sq
c pass control to FRILLS
        open (unit=8,file='/tmp/jqys.lpt',status='unknown')
	firstime = .true.
c
        CALL FRILLS (jqy_cal,
	1    nd,x,y,ysig,
	2    nptot,p,pmin,pmax,nam,
	3    xlab,ylab,xcap,ycap,subtitle,title)

        close (8)			! close fittong log file
        STOP
        END
c
c____________________________________________________________________________
c
c____________________________________________________________________________
        SUBROUTINE JQY_CAL(ycal,nv,inv,nptot,p,x,npk,ypk)
c____________________________________________________________________________
c
c Determines parameters for J(Q,Y) calculation and also checks the results for
c consistentcy between input parameters and calculated output.
c
        implicit none
c Frills Variables 
        real*8 x(*),ypk(5000,20),ycal(*),p(*),bc,bs
        integer*4 inv(*),lpt,nd,npk,nv,nptot,nrd
        real*8 rx(1000),ry(1000),ydat(1000)
        logical moms
c Frills common blocks
        common /rdats/ nrd,rx,ry
        common /dats/ ydat
        common /logf/ moms
c General Variables
        real*8 n0,kc,A,a2,a4,a6,b3,b4,b5,b6,Ab,A2b
        real*8 jqy,jqy1,jqy2,jqy3,jqy4,qq,Sq,temp
        real*8 x1(1000),temp1(1000),temp2(1000),temp3(1000),temp4(1000)
        real*8 temp5(1000),temp6(1000)
        integer*4 i,j,nord
        logical firstime,comp,resfile
c  General common blocks
        common /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6,Ab,A2b
        common /qval/ qq,Sq,temp
        common /ordr/ nord
        common /jqys/ jqy,jqy1,jqy2,jqy3,jqy4
        common /logc/ firstime,comp,resfile
c
c Assign fitting parameters from Frills
c
        bc = p(1)
        bs = p(2)
        n0 = p(3)*0.01
        kc = p(4)
        a2 = p(5)
        b3 = p(6)
        a4 = p(7)
        b4 = p(8)
        b5 = p(9)
        a6 = p(10)
        b6 = p(11)
        Ab = p(12)*0.01
        A2b = p(13)
        Sq = p(14)
c
c
c Calculate Function plus sloping background
c __ Calculation is made only at points where data is present
c
        DO i = 1,nv
           j = inv(i)
           x1(i) = x(j)
           ypk(j,1) = bc + bs*x(j)      ! Sloping background
           temp5(i) = ydat(j)
           IF (comp) THEN               ! Use option to separate J(Q,Y) comps
              CALL JQY_COMP(x(j))
              temp1(i) = jqy1
              temp2(i) = jqy2
              temp3(i) = jqy3
              temp6(i) = jqy4
              ycal(j)  = jqy            
              npk = 5
           ELSE                         ! Use option for total J(Q,Y)
              CALL JQY_TOT(x(j))
              ycal(j) = jqy
              ypk(j,2) = jqy
              temp1(i) = jqy
              npk = 2
           END IF
           temp4(i) = ycal(j)
        END DO
c
c Convolute Calculated function with Instrumental Resolution if present 
c
        IF (resfile) THEN
           IF (comp) THEN
              CALL CONVOLVE(nv,x1,temp1,nrd,rx,ry)
              CALL CONVOLVE(nv,x1,temp2,nrd,rx,ry)
              CALL CONVOLVE(nv,x1,temp3,nrd,rx,ry)
              CALL CONVOLVE(nv,x1,temp6,nrd,rx,ry)
              DO i = 1,nv
                 j = inv(i)
                 ypk(j,2) = temp1(i)
                 ypk(j,3) = temp2(i)
                 ypk(j,4) = temp3(i)
                 ypk(j,5) = temp6(i)
                 ycal(j) = temp1(i)+temp2(i)+temp3(i)+temp6(i)+ypk(j,1)
                 temp1(i) = temp1(i)+temp2(i)+temp3(i)+temp6(i) 
              END DO
           ELSE
              CALL CONVOLVE(nv,x1,temp1,nrd,rx,ry)
              DO i = 1,nv
                 j = inv(i)
                 ypk(j,2) = temp1(i)
                 ycal(j) = temp1(i) + ypk(j,1)
              END DO
           END IF
        END IF
c
c If requested, calculate and display Experimental and calculated moments
c
        IF (moms) THEN                  
           CALL JQY_MOMS(nv,x1,temp5,temp4,temp1)                
        END IF

        RETURN
        END
c
c
c_____________________________________________________________________________
        SUBROUTINE JQY_COMP(Y)
c_____________________________________________________________________________
c
c Program to calculate J(Q,Y) at a given Q as a FT of S(Q,s).
c
c ** This version separates the 3 components of J(Q,Y) which are
c       1) Condensate
c       2) Condensate Induced Singularity
c       3) Uncondensed components
c
c Therefore these components can be plotted separately
c
        implicit none
        real*8 qq,Sq,Y,Yv
        real*8 n0,kc,A,a2,a4,a6,b3,b4,b5,b6,n00,kcc,Ab,A2b
        real*8 Temp,mass,Vs,rho,kB,pi,hbar,ec,amu,te
        real*8 bound,epa,epr,wi(2000),err,wbar
        real*8 I0,const1,const2,const3,const4,const5
        real*8 jqy,jqy1,jqy2,jqy3,jqy4
        real*8 F1,F2,F3,FOFS,alpha,gama,arg,GAUSS
        integer*4 lw,liw,iw(500),inf,ifail,i,j,mn,lpt,nord,last,neval
        external F1,F2,F3,FOFS,alpha,gama,GAUSS
        logical firstime,comp,resfile,moms,init_fofs
c
        common /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6,Ab,A2b
        common /fcon/ Temp,mass,Vs,rho,kB,pi,hbar,ec,amu
        common /ycons/ const1,const2,const3,const4,const5
        common /FTan/ lw,liw,ifail,inf,bound,epa,epr
        common /jqys/ jqy,jqy1,jqy2,jqy3,jqy4
        common /ordr/ nord
        common /qval/ qq,Sq,te
        common /yval/ Yv
c
        common /logf/ moms
        common /logc/ firstime,comp,resfile
	common /ft_fs/ init_fofs
c
        Yv = Y
        temp = te
        IF (firstime) THEN              !setup some initial patameters
c
c  Define some fundamental constants
c
           mass = 4.0026
           const1 = 2.072914                       ! Ei = const1*ki**2
           const2 = const1*1.008665/mass           ! E  = const2*QQ**2
           const3 = const1*1.008665*2
           const4 = mass/(QQ*const3)               ! Y  = const4*(w-wr)
           pi = 4.0*ATAN(1.0)
           hbar = 1.0545919                        ! *10**(-34) Js
           ec   = 1.6021917                        ! *10**(-19) C
           amu  = 1.660531                         ! *10**(-27) Kg
           kB   = 1.380622                         ! *10**(-23) J/K
           rho  = 0.6022169/27.0                   ! 1/A**3         He4 at SVP
           Vs   = 240                              ! m/s            He4 at SVP
           const5 = (ec*mass*amu/hbar**2)**2/100.0
c
c Set up variables for Infinite integration
c
           liw = 500
           lw = liw*4
           ifail = 0
           inf = 1
           bound = 0.0
           epa = 0.0000001
           epr = 0.0001
c
c   n(s) = n0[1+f(s)] + A*n*(s)
c      where A = 1 - n0[1+I0]   and   I0=f(s=0)
c
           nord = 1
	   init_fofs = .true.
           arg = 0.0
           I0 = fofs(arg)
           A = 1 - n0*(1+I0) - Ab
           kcc = kc
           firstime = .false.      
        END IF                                  ! end of initialisation
        A = 1 - n0*(1+I0) - Ab
c
c  Calculate new value for A whenever kc (I0s) changes
c
        IF (kcc.NE.kc) THEN
           nord = 1
	   init_fofs = .true.
           arg = 0.0
           I0 = fofs(arg)
           A = 1 - n0*(1+I0) - Ab
           kcc = kc
        END IF
c
        CALL DQAGI(F1,bound,inf,epa,epr,jqy1,err,neval,ifail,liw,lw,last,iw,wi)
        !CALL D01AMF(F1,bound,inf,epa,epr,jqy1,err,wi,lw,iw,liw,ifail)
        jqy1 = n0*SQ/pi*jqy1
        CALL DQAGI(F2,bound,inf,epa,epr,jqy2,err,neval,ifail,liw,lw,last,iw,wi)
        !CALL D01AMF(F2,bound,inf,epa,epr,jqy2,err,wi,lw,iw,liw,ifail)
        jqy2 = n0*SQ/pi*jqy2
        CALL DQAGI(F3,bound,inf,epa,epr,jqy3,err,neval,ifail,liw,lw,last,iw,wi)
        !CALL D01AMF(F3,bound,inf,epa,epr,jqy3,err,wi,lw,iw,liw,ifail)
        jqy3 = A*SQ/pi*jqy3
        CALL DQAGI(GAUSS,bound,inf,epa,epr,jqy4,err,neval,ifail,liw,lw,last,iw,wi)
        !CALL D01AMF(F3,bound,inf,epa,epr,jqy3,err,wi,lw,iw,liw,ifail)
        jqy4 = Ab*SQ/pi*jqy4
        jqy = jqy1+jqy2+jqy3+jqy4
        RETURN
        END
c
c
*_____________________________________________________________________________
        SUBROUTINE JQY_TOT(Y)
*_____________________________________________________________________________
c
c Program to calculate J(Q,Y) at a given Q as a FT of S(Q,s).
c
c ** This version calculates the total J(Q,Y) ie does NOT separately calculates
c    the Condensate, Condensate Induced Singularity and Uncondensed components
c
c Therefore these components CANNOT be plotted separately
c
        implicit none
        real*8 Qq,Sq,Y,Yv,te
        real*8 n0,kc,A,a2,a4,a6,b3,b4,b5,b6,kcc,Ab,A2b
        real*8 Temp,mass,Vs,rho,kB,pi,hbar,ec,amu
        real*8 bound,epa,epr,wi(2000),err,wbar
        real*8 I0,const1,const2,const3,const4,const5
        real*8 jqy,jqy1,jqy2,jqy3,arg,jqy4
        real*8 F4,FOFS
        integer*4 lw,liw,iw(500),inf,ifail,i,j,mn,lpt,nord,last,neval
        external F4,FOFS
        logical firstime,comp,resfile,moms,init_fofs
c
        common /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6,Ab,A2b
        common /fcon/ Temp,mass,Vs,rho,kB,pi,hbar,ec,amu
        common /ycons/ const1,const2,const3,const4,const5
        common /FTan/ lw,liw,ifail,inf,bound,epa,epr
        common /jqys/ jqy,jqy1,jqy2,jqy3,jqy4
        common /ordr/ nord
        common /qval/ qq,Sq,te
        common /yval/ Yv
c
        common /logf/ moms
        common /logc/ firstime,comp,resfile
	common /ft_fs/ init_fofs
c
        yv = Y
        temp = te
        IF (firstime) THEN              !setup some initialize parameters
c
c  Define some fundamental constants
c
           mass = 4.0026
           const1 = 2.072914                       ! Ei = const1*ki**2
           const2 = const1*1.008665/mass           ! E  = const2*QQ**2
           const3 = const1*1.008665*2
           const4 = mass/(QQ*const3)               ! Y  = const4*(w-wr)
           pi = 4.0*ATAN(1.0)
           hbar = 1.0545919                        ! *10**(-34) Js
           ec   = 1.6021917                        ! *10**(-19) C
           amu  = 1.660531                         ! *10**(-27) Kg
           kB   = 1.380622                         ! *10**(-23) J/K
           rho  = 0.6022169/27.0                   ! 1/A**3         He4 at SVP
           Vs   = 240                              ! m/s            He4 at SVP
           const5 = (ec*mass*amu/hbar**2)**2/100.0
c
c Set up variables for Infinite integration
c
           liw = 500
           lw = liw*4
           ifail = 0
           inf = 1
           bound = 0.0
           epa = 0.0000001
           epr = 0.0001
c
c   n(s) = n0[1+f(s)] + A*n*(s)
c      where A = 1 - n0[1+I0]   and   I0=f(s=0)
c
           nord = 1
	   init_fofs = .true.
           arg = 0.0
           I0 = FOFS(arg)
           A = 1 - n0*(1+I0) - Ab
           kcc = kc
           firstime = .false.      
        END IF                                  ! end of initialisation
        A = 1 - n0*(1+I0) - Ab
c
c  Calculate new value for A whenever kc (I0) changes
c
        IF (kcc.NE.kc) THEN
           nord = 1
	   init_fofs = .true.
           arg = 0.0
           I0 = FOFS(arg)
           A = 1 - n0*(1+I0) - Ab
           kcc = kc
        END IF
c	nord = 1

        CALL DQAGI(F4,bound,inf,epa,epr,jqy,err,neval,ifail,liw,lw,last,iw,wi)
        !CALL D01AMF(F4,bound,inf,epa,epr,jqy,err,wi,lw,iw,liw,ifail)
        jqy = SQ/pi*jqy
        RETURN
        END
c
c
c_____________________________________________________________________________
        REAL FUNCTION F1*8(s)           ! alpha(s)*beta(s)
                                        ! Condensate Component
c_____________________________________________________________________________
	implicit none
        real*8 s
        real*8 ALPHA,BETA
        external ALPHA,BETA
        F1 = ALPHA(s)*BETA(s)
        RETURN
        END
c
c
c_____________________________________________________________________________
        REAL FUNCTION F2*8(s)           ! fofs(s)*alpha(s)*beta(s)
                                        ! condensate induced 'peaking' component
c_____________________________________________________________________________
	implicit none
        real*8 s
        real*8 ALPHA,BETA,FOFS
        external ALPHA,BETA,FOFS
        F2 = FOFS(s)*ALPHA(s)*BETA(s)
        RETURN
        END
c
c
c_____________________________________________________________________________
        REAL FUNCTION F3*8(s)           ! alpha(s)*gama(s)
                                        ! Uncondensed component
c_____________________________________________________________________________
	implicit none
        real*8 s
        real*8 ALPHA,GAMA
        external ALPHA,GAMA
        F3 = ALPHA(s)*GAMA(s)
        RETURN
        END
c_____________________________________________________________________________
        REAL FUNCTION F4*8(s)           ! alpha*[n0*beta + n0*beta*f + A*gama]
                                        ! Total J(Q,s)
c_____________________________________________________________________________
	implicit none
        real*8 s,n0,kc,A,a2,a4,a6,b3,b4,b5,b6,Ab,A2b
        real*8 ALPHA,BETA,GAMA,FOFS,GAUSS
        external ALPHA,BETA,GAMA,FOFS,GAUSS
        common /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6,Ab,A2b
        F4 = ALPHA(s)*( n0*BETA(s)*(1+ FOFS(s)) + A*GAMA(s) )+Ab*GAUSS(s)
        RETURN
        END
c
c
c_____________________________________________________________________________
        SUBROUTINE JQY_MOMS(npts,x1,y_expt,ycal_nores,ycal_res)
c_____________________________________________________________________________
c
c Program to calculate moments of 
c       1) Experimental data           (numerical)
c       2) Calculated J(Q,Y) lineshape (numerical)
c       3) and Expected moments from input parameters
c
c ** This serves as a consistency check between the INPUTS and OUTPUT to the
c    program
c
        implicit none
        integer*4 i,j,mn,npts,nord,l1,l2
        real*8 n0,kc,A,a2,a4,a6,b3,b4,b5,b6,Ab,A2b
        real*8 x1(1000),y_expt(1000),ycal_nores(1000),ycal_res(1000)
        real*8 M_expt(10),M_ycal_nores(10),M_ycal_res(10),M_exact(10)
        real*8 I0,I2,I4,t1,t2,t3,wbar,mom,error
        real*8 qq,Sq,temp,arg,ke,a2_ave,a2_eff,kefac
        real*8 FOFS 
        external FOFS
        character*80 line0,line1,line2,file1*60,file2*60
        logical firstime,comp,resfile,moms,ft_fofs
c
        common /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6,Ab,A2b
        common /ordr/ nord
        common /file/ file1,l1,file2,l2
        common /qval/ qq,Sq,temp
c
        common /logf/ moms
        common /logc/ firstime,comp,resfile
        common /ft_fs/ ft_fofs
c
        kefac = 1.5*1.043*11.6045
c
c Calculate Moments of experimental data
c
        DO i=1,7
           mn = i-1
           wbar = 0.0
           CALL MOMENTS(x1,y_expt,npts,mn,mom,error,wbar)
           M_expt(i) = mom
        END DO
c
c Calculate Moments of Evaluated function without resolution
c
        DO i=1,7
           mn = i-1
           wbar = 0.0
           CALL MOMENTS(x1,ycal_nores,npts,mn,mom,error,wbar)
           M_ycal_nores(i) = mom
        END DO
c
c Calculate Moments of Evaluated function after convolution with instr. Reso.
c
c       IF (resfile) THEN
        DO i=1,7
           mn = i-1
           wbar = 0.0
           CALL MOMENTS(x1,ycal_res,npts,mn,mom,error,wbar)
           M_ycal_res(i) = mom
        END DO
c       END IF
c
c Evaluate the expected (Exact) moments from the input parameters used to
c calculate the fitting function
c
        arg = 0.0
        nord = 3
        ft_fofs = .true.
        I2 = FOFS(arg)
        nord = 5
        ft_fofs = .true.
        I4 = FOFS(arg)
        nord = 1
        ft_fofs = .true.
        I0 = FOFS(arg)
c
        M_exact(1) = 1.0                                        ! Zeroth moment
        M_exact(2) = 0.0                                        ! First moment
        M_exact(3) = n0*I2/3.0 + A*a2                           ! Second moment
        M_exact(4) = b3                                         ! Third moment
        M_exact(5) = b4 + n0*I4/5 + A*(3*a2**2+a4)              ! Fourth moment
        M_exact(6) = b5 + 10.0*(n0*I2/3*b3 + A*a2*b3)           ! Fifth moment
        t1 = b6 + A*a6 + 10*b3**2       !** neglect term n0I6/7 **
        t2 = 15.0*b4*(n0*I2/3 + A*a2)
        t3 = 15.0*A*a2*(a2**2 + a4)
        M_exact(7) = t1+t2+t3                           ! Sixth moment 
c        ke = A*3/2*1.043*11.6045*a2 +              n0*I2
        a2_eff = A*a2+Ab*a2b
        a2_ave = a2_eff/(A+Ab)
        ke = kefac*( (A+Ab)*a2_ave + n0*I2 ) 
c
        write (*,*) '#_________________________________________________'
        write (*,*) '#_________________________________________________'
        write (*,'(A)') ' # ** Data file = '//file1(1:l1)
        write (*,'(A)') ' # ** Resn file = '//file2(1:l2)
        write (*,*) '# '
        write (*,'(3(A,F7.3))') ' # ** Q  = ',QQ,'  ** S(Q) = ',Sq,'  ** Temp. = ',temp    
	write (*,'(4(A,F7.3))') ' # ** n0 = ',n0,'  **  Io  = ',I0,'  **    A  = ',A,'  ** Ab = ',Ab 
        write (*,'(2(A,F6.3))') ' # ** I2 = ',I2,'  **  KE  = ',ke
        write (*,'(3(A,F6.3))') ' # ** a2 = ',a2,'  ** a2b  = ',a2b,'  ** combined = ',a2_ave
        write (*,*) '# '
        line0 = ' #****** Moments evaluated are as follows  : '
        write (*,*) '# '
        line1 = '                  DATA        CALCULATION    CALN+INST RES        EXACT   '
        line2 = '                  ----        -----------    -------------        -----   '
        write (*,'(A)') line0
        write (*,'(A)') line1
        write (*,'(A)') line2
        write (*,20) 'M0 = ',M_expt(1),M_ycal_nores(1),M_ycal_res(1),M_exact(1)
        write (*,20) 'M1 = ',M_expt(2),M_ycal_nores(2),M_ycal_res(2),M_exact(2)
        write (*,20) 'M2 = ',M_expt(3),M_ycal_nores(3),M_ycal_res(3),M_exact(3)
        write (*,20) 'M3 = ',M_expt(4),M_ycal_nores(4),M_ycal_res(4),M_exact(4)
        write (*,20) 'M4 = ',M_expt(5),M_ycal_nores(5),M_ycal_res(5),M_exact(5)
        write (*,20) 'M5 = ',M_expt(6),M_ycal_nores(6),M_ycal_res(6),M_exact(6)
        write (*,20) 'M6 = ',M_expt(7),M_ycal_nores(7),M_ycal_res(7),M_exact(7)
 20     FORMAT (T5,A5,5X,G12.4,3X,G12.4,4X,G12.4,7X,G12.4)
        write (*,*) '# '
        write (*,'(A)') ' #** NOTE: A term (n0*I6/7) was neglected from the exact M6 '
        write (*,'(A)') ' #   _____ I6 is the sixth moment of f(s) '
        write (*,*) '#__________________________________________________'
        write (*,*) '#__________________________________________________'
        write (*,*) '#'


        RETURN
        END


c
c____________________________________________________________________________
        REAL FUNCTION GAUSS*8(s)
c____________________________________________________________________________
c A Gaussian function to represent scattering from tightly bound helium layer atoms
c in confined geometries
c
        implicit none
        real*8 n0,kc,A,a2,a4,a6,b3,b4,b5,b6,Ab,A2b
        real*8 s,k,arg,Y
        common /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6,Ab,A2b
        common /yval/ Y
        arg = -A2b*s**2/2
        IF (ABS(arg).LE.60.0) THEN
           GAUSS = COS(Y*s)*EXP(arg)
        ELSE
           GAUSS = 0.0
        END IF
        RETURN
        END



c
c____________________________________________________________________________
        REAL FUNCTION GAMA*8(s)
c____________________________________________________________________________
c The even function component of the
c non-complex Exponential part of S(Q,s)
c
        implicit none
        real*8 n0,kc,A,a2,a4,a6,b3,b4,b5,b6,Ab,A2b
        real*8 s,k,arg
        common /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6,Ab,A2b
        arg = -a2*s**2/2+(a4+b4)*s**4/24-(a6+b6)*s**6/720
        IF (ABS(arg).LE.60.0) THEN
           GAMA = EXP(arg)
        ELSE
           GAMA = 0.0
        END IF
        RETURN
        END

c____________________________________________________________________________
        REAL FUNCTION BETA*8(s)
c____________________________________________________________________________
c The even function component of the
c non-complex Exponential part of R(Q,s)
c
        implicit none
        real*8 n0,kc,A,a2,a4,a6,b3,b4,b5,b6,Ab,A2b
        real*8 s,k,arg
        common /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6,Ab,A2b
        arg = b4*s**4/24-b6*s**6/720
        IF (ABS(arg).LE.60.0) THEN
           BETA = EXP(arg)
        ELSE
           BETA = 0.0
        END IF
        RETURN
        END

*____________________________________________________________________________
        REAL FUNCTION ALPHA*8(s)
*____________________________________________________________________________
c The even function component of the
c complex Exponential part of S(Q,s)
c
        implicit none
        real*8 n0,kc,A,a2,a4,a6,b3,b4,b5,b6,Ab,A2b
        real*8 s,k,arg,Y
        common /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6,Ab,A2b
        common /yval/ Y
        ALPHA = COS(Y*s+b3*s**3/6-b5*s**5/120)
        RETURN
        END

c
*--------------------------------------------------------------------   
        REAL FUNCTION FOFS*8(s)
*--------------------------------------------------------------------   
c
c To Calculate f(s).
c f(s) represents the condensate-induced singularity in the momentum dist.
c
c     n(s) = n0[1+f(s)] + A*n*(s)
c
c  variables
c   s = distance travelled by the struck particle.
c   k = 1-D atomic momentum
c 
        implicit none
        real*8 n0,kc,A,a2,a4,a6,b3,b4,b5,b6,Ab,A2b
        real*8 s,k
        real*8 Temp,mass,Vs,rho,kB,pi,hbar,ec,amu
        real*8 bound,epa,epr,wi(2000),err,wbar
        real*8 I0,FS
        integer*4 lw,liw,iw(500),inf,ifail,i,j,mn,lpt,nord,last,neval
        external FS
c  General common blocks
        common /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6,Ab,A2b
        common /fcon/ Temp,mass,Vs,rho,kB,pi,hbar,ec,amu
        common /FTan/ lw,liw,ifail,inf,bound,epa,epr 
        common /ordr/ nord
c  FS variables and common blocks
        real*8 sv,B,D,F,pii
        common /fsc1/ pii,B,D,F,sv
c Interpolation routine variables and common blocks
        logical ft_fofs,spline
        integer MSP,lwrk,ierr
        real*8 x(1000),y(1000),SM,wrk(2000),deri(1000)
        real*8 val,res,s_max,s_min
        character*1 start
        common /ft_fs/ ft_fofs
c
        s_min = 0.0
        s_max = 100.0
        spline = .false.
        ierr = 0
c
        IF (s .gt. s_max.OR.s.LT.s_min) THEN
           FOFS = 0.0    
           RETURN
        END IF        
        IF (ft_fofs) THEN
           ft_fofs = .false.
           pii = pi
           B = mass*amu*Vs/(2.0*hbar*(2.0*pi)**3*rho)/1000.0       ! A**2
           D = hbar*Vs/(2.0*kb*temp)/10.0                          ! A
           F = 2.0*kc**2                                           ! (1/A)**2
           !Make a sample of FOFS for s = 0 to 100 in steps of 0.12
           !Set up a cubic spline for this from which FOFS can be deterimined for any s
	   MSP = 260
           DO i = 1,MSP
              x(i) = (i-1)*0.40		! max s = 110.0
              sv = x(i)
              CALL DQAGI(FS,bound,inf,epa,epr,I0,err,neval,ifail,liw,lw,last,iw,wi)
              y(i) = I0
           END DO
c Produce a smoothing cubic spline function
           lwrk = 2*MSP
           CALL DPCHEZ(MSP,X,Y,DERI,SPLINE,WRK,LWRK,IERR)           
           !CALL E02BEF(start,MSP,X,Y,W,SM,NEST,NSP,KSP,CSP,FP,WRK,LWRK,IWRK,IFA)
        end if
        CALL DPCHEV (MSP,X,Y,DERI,1,s,res,wrk,IERR)
        !CALL E02BBF(NSP,KSP,CSP,s,res,ifa)
        FOFS = res
        RETURN
        END                                    
c
*--------------------------------------------------------------------   
        REAL FUNCTION FS*8(k) 
*--------------------------------------------------------------------   
        implicit none
        real*8 s,k,pi,B,D,F,COTH,SINC
        real*8 arg1,arg2,arg3
        integer*4 nord
        external COTH,SINC
        common /ordr/ nord
        COMMON /fsc1/ pi,B,D,F,s
c
        arg1 = D*k
	IF (F.EQ.0.0) THEN
	   FS = 0.0
	   RETURN
	ELSE 
           arg2 = -k**2/F
	END IF
        arg3 = k*s
        IF (k.EQ.0.0) THEN
           FS = 0
        ELSE IF (s.EQ.0.0) THEN
           FS = 4*pi*B*(k**nord)*COTH(arg1)*EXP(arg2)
        ELSE
           FS = 4*pi*B*(k**nord)*COTH(arg1)*EXP(arg2)*SINC(arg3)
        END IF
        RETURN
        END
c
c Define a HYPERBOLIC COT
c
c____________________________________________________________________
        REAL FUNCTION COTH*8(X)
c____________________________________________________________________
        implicit none
        real*8 x
c
        IF (x.GT.10.0) THEN
           COTH = 1.0
           RETURN
        END IF
        IF (x.EQ.0) THEN
           COTH = 0.0
        ELSE 
           COTH = (EXP(x)+EXP(-1*x))/(EXP(x)-EXP(-1*x))
        END IF
        RETURN
        END
c
c Define Sin(x)/x
c____________________________________________________________________
c____________________________________________________________________
        REAL FUNCTION SINC*8(X)
c____________________________________________________________________
        implicit none
        real*8 x
c
        IF (x.EQ.0.0) THEN
           SINC = 1.0
           RETURN
        END IF
        SINC = SIN(x)/x
        RETURN
        END








































