*______________________________________________________________________________
        PROGRAM JQYJ_FIT
*______________________________________________________________________________
c
c # Fits a model J(Q,Y) to experimental data.
c # J(Q,Y) = Jia(Q,Y)oR(Q,Y)
c # Jia(Q,Y) obtained as FT of Jia(Q,s), the intermediate scattering function
c # Jia(Q,s) = n0 + n0*f(s) + n'(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,Y) = Final State function
c 
c         In this version, n'(s) and R(Q,Y) are given.
c         Only n0 and f(s) are to be determined.
c
c R.T. AZUAH -- April 1999
c______________________________________________________________________________
	implicit none
	external jqyj_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,filefsf
	integer*4 text_len,lpt,strim,nrd,nid
	logical okay,resfile,moms,comp,firstime
	integer*4 nptot,nd,i,j,in_len,res_len,mn,fsf_len
	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),ix(1000),iy(1000)
c
        common /dats/ y
        common /rdats/ nrd,rx,ry
        common /nofk/  nnp,npx,npy
        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) /'a4   [\A\u-4\d]'/,
     >       nam(7) /'a6   [\A\u-6\d]'/
        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 = 7
 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 value, S(Q) and Temp of data ::'
        read (*,*) qq,Sq,temp
c       Read in Instrumental Resolution function data
 20     write (*,'(T5,A,$)') ' Enter Instrument 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) ix(i),iy(i),er
           END DO
 25        nid = i - 1  
           close (18)
        ELSE
           write (*,*) '    ** ERROR - Instr Resolution file does not exist '
           GOTO 20
        END IF
c       Read in final state Function data
 26     write (*,'(T5,A,$)') ' Enter Final State Function file ::'
        read (*,'(A)') filefsf
        INQUIRE (file=filefsf,exist=okay)
        IF (okay) THEN
           open (unit=18,file=filefsf,status='old')
           fsf_len = strim(filefsf)
           DO i = 1,1000
              read (18,*,end=27) rx(i),ry(i)
           END DO
 27        nrd = i - 1  
           close (18)
        ELSE
           write (*,*) '    ** ERROR - Final State file does not exist '
           GOTO 26
        END IF
c       Convolute Final State function with Instrumental Resolution 
        call convolve(nrd,rx,ry,nid,ix,iy)
c        open (unit=10,file='tmp.dat',status='unknown')
c        do i = 1,nrd
c           write (10,'(3G12.4)') rx(i),ry(i),0.05*ry(i)
c        end do
c        close (10)
c        write (*,*) 'Convoluted R(Q,Y) stored in :: tmp.dat'
c
        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) '
c        write (subtitle(1)(13:17), '(i5)') 1
	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 (*,*) '# ******** JQYJ_FIT Fitting Program ****** '
        write (*,*) '#'
        write (*,*) '# Type HELP for a list of available commands'
        write (*,*) '#'
        write (*,*) '# Details of fits are stored in /tmp/jqyj_fit.lpt '
        write (*,*) '#'
        write (*,*) '#'
c        call norm (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) = 0.46
        p(7) = 0.38
c pass control to FRILLS
        open (unit=8,file='/tmp/jqyj_fit.lpt',status='unknown')
	firstime = .true.
c
c        write (*,*) ' JQYS: firstime,comp,resfile = ',firstime,comp,resfile
        CALL FRILLS (jqyj_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
*____________________________________________________________________________
        SUBROUTINE JQYJ_CAL(ycal,nv,inv,nptot,p,x,npk,ypk)
*_____________________________________________________________________________
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(1000),ypk(1000,20),ycal(1000),p(50),bc,bs
	integer*4 inv(1000),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,F1
        real*8 jqy,jqy1,jqy2,jqy3,qq,Sq,temp,y1(1000)
        real*8 x1(1000),temp1(1000),temp2(1000),temp3(1000),temp4(1000)
        integer*4 i,j,nord
        logical firstime,comp,resfile
        external F1
c  General common blocks
        common /parm/ n0,kc,A,a2,a4,a6
        common /qval/ qq,Sq,temp
        common /ordr/ nord
        common /jqys/ jqy,jqy1,jqy2,jqy3
        common /logc/ firstime,comp,resfile

c        write (*,*) ' JQY_CAL: firstime,comp,resfile = ',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)
        a4 = p(6)
        a6 = p(7)
c
c
c Calculate Function plus sloping background
c __ Calculation is made only at points where data is present
c
c        DO i = 1,nv
c           j = inv(i)
c           x1(i) = x(j)
c           y1(i) = ydat(j)
c        END DO
c        CALL REB_DAT1(nrd,rx,ry,temp4,nv,x1,temp1,temp4)   ! rebin R(Q,Y) at Y values of data
                                                           !  this is used for n0*R(Q,Y) comp of J(Q,Y)
c
        DO i = 1,nv
           j = inv(i)
           x1(i) = x(j)
           y1(i) = ydat(j)
           ypk(j,1) = bc + bs*x(j)      ! Sloping background
           temp1(i) = n0*F1(x1(i))      ! condensate comp {n0*R(Q,Y)}
           IF (comp) THEN               ! Use option to separate J(Q,Y) comps
              CALL JQY_COMP(x(j))
              temp2(i) = jqy2
              temp3(i) = jqy3
              ycal(j)  = jqy            
              npk = 4
           ELSE                         ! Use option for total J(Q,Y)
              CALL JQY_TOT(x(j))
              ycal(j) = jqy + temp1(i)
              ypk(j,2) = jqy + temp1(i)
              temp2(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
c              CALL CONVOLVE(nv,x1,temp1,nrd,rx,ry)      ! {conv not required for n0 comp.}
              CALL CONVOLVE(nv,x1,temp2,nrd,rx,ry)
              CALL CONVOLVE(nv,x1,temp3,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)
                 ycal(j) = temp1(i)+temp2(i)+temp3(i)+ypk(j,1)
		 temp1(i) = temp1(i)+temp2(i)+temp3(i) 
              END DO
           ELSE
              CALL CONVOLVE(nv,x1,temp2,nrd,rx,ry)
              DO i = 1,nv
		 j = inv(i)
                 ypk(j,2) = temp1(i) + temp2(i)
                 ycal(j) = ypk(j,2) + ypk(j,1)
                 temp1(i) = ycal(j)
              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,y1,temp4,temp1)                
        END IF

        RETURN
        END
c
c
*_____________________________________________________________________________
        SUBROUTINE JQY_COMP(Y)
*_____________________________________________________________________________
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,n00,kcc
        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
        real*8 F1,F2,F3,FOFS,alpha,gama,arg
        integer*4 lw,liw,iw(500),inf,ifail,i,j,mn,lpt,nord,last,neval
        external F1,F2,F3,FOFS,alpha,gama
        logical firstime,comp,resfile,moms,init_fofs
c
        common /parm/ n0,kc,A,a2,a4,a6
        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
        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)
           kcc = kc
           firstime = .false.      
        END IF                                  ! end of initialisation
        A = 1 - n0*(1+I0)
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)
           kcc = kc
        END IF
	nord = 1


c        CALL DQAGI(F1,bound,inf,epa,epr,jqy1,err,neval,ifail,liw,lw,last,iw,wi)
c        jqy1 = n0*F1(Y)
        CALL DQAGI(F2,bound,inf,epa,epr,jqy2,err,neval,ifail,liw,lw,last,iw,wi)
        jqy2 = n0*SQ/pi*jqy2
        CALL DQAGI(F3,bound,inf,epa,epr,jqy3,err,neval,ifail,liw,lw,last,iw,wi)
        jqy3 = A*SQ/pi*jqy3
        jqy = jqy1+jqy2+jqy3
        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,kcc
        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
        real*8 F4,FOFS
c        real D01AMF
        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
        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
        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
c           write (*,*) ' JQY_TOT : Entering first time initialisation '
           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)
           kcc = kc
           firstime = .false.      
        END IF                                  ! end of initialisation
        A = 1 - n0*(1+I0)
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)
           kcc = kc
        END IF
	nord = 1
        CALL DQAGI(F4,bound,inf,epa,epr,jqy,err,neval,ifail,liw,lw,last,iw,wi)
        jqy = SQ/pi*jqy
        RETURN
        END
c
c
*_____________________________________________________________________________
        REAL FUNCTION F1*8(Y)           ! find value of R(Q,Y) at Y by interpolation
*_____________________________________________________________________________
        implicit none
        real*8 Y
c Interpolation routine variables and common blocks
        logical firstime /.true./,spline /.false./
        integer MSP,lwrk,ierr,nrd,nout
        real*8 rx(1000),ry(1000),re(1000),xout(1),yout(1),SM,wrk(2000),deri(1000),rqy
        character*1 start
        common /rdats/ nrd,rx,ry

c        write (*,*) 'F1: Calling F1'
        if (firstime) then
c           write (*,*) 'F1: nrd = ',nrd
           lwrk = 2*nrd
           CALL DPCHEZ(nrd,rx,ry,deri,spline,wrk,lwrk,ierr)
           firstime = .false.
        end if
        CALL DPCHEV (nrd,rx,ry,deri,1,Y,rqy,wrk,ierr)
        F1 = rqy

        return
        end
c
*_____________________________________________________________________________
        REAL FUNCTION F2*8(s)           ! fofs(s)*alpha(s)
                                        ! condensate induced 'peaking' component
*_____________________________________________________________________________
	implicit none
        real*8 s
        real*8 ALPHA,BETA,FOFS
        external ALPHA,BETA,FOFS
        F2 = FOFS(s)*ALPHA(s)
        RETURN
        END
c
c
*_____________________________________________________________________________
        REAL FUNCTION F3*8(s)           ! alpha(s)*gama(s)
                                        ! Uncondensed component
*_____________________________________________________________________________
	implicit none
        real*8 s
        real*8 ALPHA,GAMA
        external ALPHA,GAMA
        F3 = ALPHA(s)*GAMA(s)
        RETURN
        END
*_____________________________________________________________________________
        REAL FUNCTION F4*8(s)           ! alpha*[n0*f + A*gama]            {n0 comp included later}
                                        ! Total J(Q,s)
*_____________________________________________________________________________
	implicit none
        real*8 s,n0,kc,A,a2,a4,a6
        real*8 ALPHA,BETA,GAMA,FOFS
        external ALPHA,BETA,GAMA,FOFS
        logical first /.true./
        common /parm/ n0,kc,A,a2,a4,a6
c        if (first) then
c           first = .false.
c           write (*,*) 'F4: n0,A = ',n0,A
c        end if
        F4 = ALPHA(s)*(n0* FOFS(s) + A*GAMA(s) )
        RETURN
        END
c
c
*_____________________________________________________________________________
        SUBROUTINE JQY_MOMS(npts,x1,y_expt,ycal_nores,ycal_res)
*_____________________________________________________________________________
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
        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
        real*8 FOFS 
        external FOFS
        character*80 line0,line1,line2,file1*60,file2*60
        logical firstime,comp,resfile,moms
c
        common /parm/ n0,kc,A,a2,a4,a6
        common /ordr/ nord
        common /file/ file1,l1,file2,l2
        common /qval/ qq,Sq,temp
c
        common /logf/ moms
        common /logc/ firstime,comp,resfile

        b3=0.0; b4=0.0; b5=0.0; b6=0.0
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
c Evaluate the expected (Exact) moments from the input parameters used to
c calculate the fitting function
c
        nord = 1
        arg = 0.0
        I0 = FOFS(arg)
        nord = 3
        I2 = FOFS(arg)
        nord = 5
        I4 = 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
        write (*,*) '#_____________________________________________________________________________'
        write (*,'(A)') ' # ** Data file = '//file1(1:l1)
        write (*,'(A)') ' # ** Resn file = '//file2(1:l2)
        write (*,*) '# '
        write (*,'(3(A,F6.3))') ' # ** Q = ',QQ,'  ** S(Q) = ',Sq,'  ** Temp. = ',temp    
	write (*,'(3(A,F6.3))') ' # ** n0= ',n0,'  **  Io  = ',I0,'  **    A  = ',A 
        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 (*,*) '#_____________________________________________________________________________'
        RETURN
        END

c
c_____________________________________________________________________________________
c
*____________________________________________________________________________
        REAL FUNCTION GAMA*8(s)
*____________________________________________________________________________
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
        real*8 s,k,arg
        common /parm/ n0,kc,A,a2,a4,a6
        arg = -a2*s**2/2 + a4*s**4/24 - a6*s**6/720
        IF (ABS(arg).LE.60.0) THEN
           GAMA = EXP(arg)
        ELSE
           GAMA = 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
        real*8 s,k,arg,Y
        common /yval/ Y
        ALPHA = COS(Y*s)
        RETURN
        END
c
c_______________________________________________________________________________
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
        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
	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
c Set up variables for Infinite integration
c
c        liw = 100
c        lw = liw*4
c        ifail = 0
c        inf = 1
c        bound = 0.0
c        epa = 0.0
c        epr = 0.0001
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
c           write (*,*) ' FOFs : mass,amu,Vs,hbar= ',mass,amu,Vs,hbar
           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
c
c           write (*,*) ' FOFs : B,D,F pii= ',B,D,F,pii
c
c Make a sample of FOFS for s = 0 to 100 in steps of 0.12
c Set up a cubic spline for this from which FOFS can be deterimined for any s
c
	   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
c              write (*,'(3G12.4)') x(i),y(i),0.001
           END DO
c Produce a smoothing cubic spline function
           lwrk = 2*MSP
           CALL DPCHEZ(MSP,X,Y,DERI,SPLINE,WRK,LWRK,IERR)           
c           CALL E02BEF(start,MSP,X,Y,W,SM,NEST,NSP,KSP,CSP,FP,WRK,LWRK,IWRK,IFA)
        end if
c        write (*,*) 'MSP,Y(1),Y(10) = ',MSP,y(1),y(10)
        CALL DPCHEV (MSP,X,Y,DERI,1,s,res,wrk,IERR)
c        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
*--------------------------------------------------------------------   
        REAL FUNCTION COTH*8(X)
*--------------------------------------------------------------------   
        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
*--------------------------------------------------------------------   
        REAL FUNCTION SINC*8(X)
*--------------------------------------------------------------------   
        implicit none
        real*8 x
c
        IF (x.EQ.0.0) THEN
           SINC = 1.0
           RETURN
        END IF
        SINC = SIN(x)/x
        RETURN
        END








































