*______________________________________________________________________________
        PROGRAM GEN_JQY
*______________________________________________________________________________
c
c # Calculate J(Q,Y) given parameters.
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 -- Sep 98

c
c______________________________________________________________________________
c
        implicit none
	character*30 f1 /'jqy.dat'/,f2 /'jqy_n0.dat'/,f3 /'jqy_fk.dat'/
        character*30 f4 /'jqy_unc.dat'/,dec*1
	integer*4 strim, u1,u2,u3,u4
	logical okay,moms,comp,firstime,resfile
	integer*4 i,j,out_len
c	real*8 x(1000),y(1000),ysig(1000)
        real*8 qq,Sq,temp,y,miny,maxy,dy
        real*8 n0,kc,A,bc,bs,a2,a4,a6,b3,b4,b5,b6,jqy,jqy1,jqy2,jqy3
c
        common /qval/ qq,Sq,temp
        common /logf/ moms
        common /logc/ firstime,comp,resfile
        common /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6
        common /jqys/ jqy,jqy1,jqy2,jqy3
c
        comp = .false.
        resfile = .false.
        firstime = .true.
        u1=10; u2=20; u3=30; u4=40
        write (*,'(T5,A,$)') 'Enter min, max and interval of Y ::'
        read (*,*) miny,maxy,dy
        write (*,'(T5,A,$)') ' Enter Q value, S(Q) and Temp of data ::'
        read (*,*) qq,Sq,temp
c        write (*,'(T5,A)') ' Decide now whether to keep components of J(Q,Y) separate '
        if (temp < 2.17) then
           write (*,'(T5,A,$)') 'Output Components? Default [N] ::'
           read (*,'(A)') dec
           IF (dec.EQ.'y'.OR.dec.EQ.'Y')  comp = .true.
        end if
        open (unit=u1, file=f1, status='unknown')
        if (comp) then
           open (unit=u2, file=f2, status='unknown')
           open (unit=u3, file=f3, status='unknown')
           open (unit=u4, file=f4, status='unknown')
        end if
c
        n0 = 0.07
        kc = 0.5
        a2 = 0.897  
        b3 = 3.75
        a4 = 0.46 
        b4 = 0.0
        b5 = 2500.0
        a6 = 0.38
        b6 = 200.0
        write (*,*) ' '
        write (*,'(T5,A,G10.4,T30,A,G10.4)') '*** n0 = ',n0,'*** kc = ',kc
        write (*,'(T5,A,G10.4,T30,A,G10.4)') '*** a2 = ',a2,'*** b3 = ',b3
        write (*,'(T5,A,G10.4,T30,A,G10.4)') '*** a4 = ',a4,'*** b5 = ',b5
        write (*,'(T5,A,G10.4,T30,A,G10.4)') '*** a6 = ',a6,'*** b6 = ',b6

        write (*,'(T5,A,$)') 'Enter n0 (zero for normal phase): '
        read (*,*) n0
        write (*,'(T5,A,$)') 'Alter anything else? (def=N): '
        read (*,'(A)') dec
        if (dec.EQ.'Y' .OR. dec.EQ.'y') then
           write (*,'(T5,A,$)') 'Enter kc: '
           read (*,*) kc
           write (*,'(T5,A,$)') 'Enter a2: '
           read (*,*) a2
           write (*,'(T5,A,$)') 'Enter b3: '
           read (*,*) b3
           write (*,'(T5,A,$)') 'Enter a4: '
           read (*,*) a4
           write (*,'(T5,A,$)') 'Enter b5: '
           read (*,*) b5
           write (*,'(T5,A,$)') 'Enter a6: '
           read (*,*) a6
           write (*,'(T5,A,$)') 'Enter b6: '
           read (*,*) b6
        end if
c     convert parameters to appropriate units
        a2 = a2                         ! [(1/A)**2]
        b3 = b3/QQ                      ! [(1/A)**3]
        a4 = a4                         ! [(1/A)**4]
        b4 = 0.0                        ! [(1/A)**4]
        b5 = b5/QQ**3                   ! [(1/A)**5]
        a6 = a6                         ! [(1/A)**6]
        b6 = b6/QQ**2                   ! [(1/A)**6]


c
	firstime = .true.
c
c Calculate J(Q,Y) 
c
        y = miny
        maxy = maxy - 0.5*dy
        do while (y .LT. maxy)
           IF (comp) THEN               ! Use option to separate J(Q,Y) comps
              CALL JQY_COMP(y)
              write (u1,'(2G11.4)') y,jqy
              write (u2,'(2G11.4)') y,jqy1
              write (u3,'(2G11.4)') y,jqy2
              write (u4,'(2G11.4)') y,jqy3
           ELSE                         ! Use option for total J(Q,Y)
              CALL JQY_TOT(y)
              write (u1,'(2G11.4)') y,jqy
           END IF
           y = y + dy
        end do
        close (u1)
        u1 = strim(f1)
        write (*,*) '##################################################################'
        write (*,'(A)') ' #####   J(Q,Y) stored in ::'//f1(1:u1)//'                           #####' 
        if (comp) then
           close (u2)
           close (u3)
           close (u4)
        u2 = strim(f2)
        write (*,'(A)') ' #####   n0  component of  J(Q,Y) stored in ::'//f2(1:u2)//'      #####'
        u3 = strim(f3)
        write (*,'(A)') ' #####   f(k) component of J(Q,Y) stored in ::'//f3(1:u3)//'      #####' 
        u4 = strim(f4)
        write (*,'(A)') ' #####   Uncondensed comp  J(Q,Y) stored in ::'//f4(1:u4)//'     #####' 
        write (*,*) '##################################################################'
        end if
c
        STOP
        END
c
c _________________________________________________________________________________________________
c
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,b3,b4,b5,b6,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,moms,init_fofs,resfile
c
        common /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6
        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
c       write (*,*) ' jqy_comp : Fofs(0) = ',Fofs(0.0)
c       write (*,*) ' jqy_comp : Fofs(2) = ',Fofs(2.0)
c       write (*,*) ' jqy_comp : Fofs(5) = ',Fofs(3.0)
	nord = 1


        CALL DQAGI(F1,bound,inf,epa,epr,jqy1,err,neval,ifail,liw,lw,last,iw,wi)
c        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)
c        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)
c        CALL D01AMF(F3,bound,inf,epa,epr,jqy3,err,wi,lw,iw,liw,ifail)
        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,b3,b4,b5,b6,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,b3,b4,b5,b6
        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)
c        CALL D01AMF(F4,bound,inf,epa,epr,jqy,err,wi,lw,iw,liw,ifail)
c       write (*,*) ' JQY_TOT : jqy = ',jqy
        jqy = SQ/pi*jqy
        RETURN
        END
c
c
*_____________________________________________________________________________
        REAL*8 FUNCTION F1(s)           ! alpha(s)*beta(s)
                                        ! Condensate Component
*_____________________________________________________________________________
	implicit none
        real*8 s
        real*8 ALPHA,BETA
        external ALPHA,BETA
        F1 = ALPHA(s)*BETA(s)
        RETURN
        END
c
c
*_____________________________________________________________________________
        REAL*8 FUNCTION F2(s)           ! fofs(s)*alpha(s)*beta(s)
                                        ! condensate induced 'peaking' component
*_____________________________________________________________________________
	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
*_____________________________________________________________________________
        REAL*8 FUNCTION F3(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*8 FUNCTION F4(s)           ! alpha*[n0*beta + n0*beta*f + A*gama]
                                        ! Total J(Q,s)
*_____________________________________________________________________________
	implicit none
        real*8 s,n0,kc,A,a2,a4,a6,b3,b4,b5,b6
        real*8 ALPHA,BETA,GAMA,FOFS
        external ALPHA,BETA,GAMA,FOFS
        common /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6
        F4 = ALPHA(s)*( n0*BETA(s)*(1+ FOFS(s)) + A*GAMA(s) )
        RETURN
        END
c
c_____________________________________________________________________________________
c
*____________________________________________________________________________
        REAL*8 FUNCTION GAMA(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,b3,b4,b5,b6
        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

*____________________________________________________________________________
        REAL*8 FUNCTION BETA(s)
*____________________________________________________________________________
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
        real*8 s,k,arg
        common /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6
        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*8 FUNCTION ALPHA(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 /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6
        common /yval/ Y
        ALPHA = COS(Y*s+b3*s**3/6-b5*s**5/120)
        RETURN
        END
c
c_______________________________________________________________________________
c
*--------------------------------------------------------------------   
        REAL*8 FUNCTION FOFS(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
        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
	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*8 FUNCTION FS(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*8 FUNCTION COTH(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*8 FUNCTION SINC(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








































