*_____________________________________________________________________________
       PROGRAM GEN_NP
*_____________________________________________________________________________
c
c Program to calculate n(p) .
c
        implicit none
        real*8 Qq,Sq,k,te,B,D,F,sv,ppi,jy
        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(4000),err,wbar
        real*8 I0,const1,const2,const3,const4,const5
        real*8 nk,mink,maxk,dk,fkk,s,ke,I2
        real*8 F4,FOFS,FK
        integer*4 lw,liw,iw(1000),inf,ifail,i,j,mn,lpt,nord,last,neval,strim
        character*60 ofile,ofile1,dec*1
        external F4,FOFS,FK
        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 /ordr/ nord
        common /qval/ qq,Sq,te
        common /yval/ k
c
        common /fsc1/ ppi,B,D,F,sv
	common /ft_fs/ init_fofs
c
        temp = 1.6
        kc = 0.5
        n0 = 0.06
        a2 = 0.897
        a4 = 0.46
        a6 = 0.38
        qq = 29.0
        sq = 1.0
        
        write (*,*) '********** Default n(p) Parameters ****************'
        write (*,'(T5,A,G10.4,T30,A,G10.4)') '*** n0 = ',n0,'* temp = ',temp
        write (*,'(T5,A,G10.4,T30,A,G10.4)') '*** kc = ',kc,'*** a2 = ',a2
        write (*,'(T5,A,G10.4,T30,A,G10.4)') '*** a4 = ',a4,'*** a6 = ',a6
        write (*,*) ' '
        write (*,'(T5,A,$)') 'Enter n0 and temp: '
        read (*,*) n0,temp
        write (*,'(T5,A,$)') 'Alter anything else? (def=N): '
        read (*,'(A)') dec
        if (dec.EQ.'Y' .OR. dec.EQ.'y') then
           write (*,'(T5,A,$)') 'Enter Sq: '
           read (*,*) Sq
           write (*,'(T5,A,$)') 'Enter kc: '
           read (*,*) kc
           write (*,'(T5,A,$)') 'Enter a2: '
           read (*,*) a2
           write (*,'(T5,A,$)') 'Enter a4: '
           read (*,*) a4
           write (*,'(T5,A,$)') 'Enter a6: '
           read (*,*) a6
        end if

	write (*,'(T5,A,$)') 'Enter min, max and interval for p:'
	read (*,*) mink,maxk,dk
	write (*,'(T5,A,$)') 'Give output data file :'
	read (*,'(A)') ofile 
	open (unit=10,file=ofile,status='unknown')
	INQUIRE (unit=10,name=ofile)

c        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
        ppi = pi
c     
c     Set up variables for Infinite integration
c     
        liw = 100
        lw = liw*4
        ifail = 0
        inf = 1
        bound = 0.0
        epa = 0.00001
        epr = 0.001
c
c   n(s) = n0[1+f(s)] + A*n*(s)
c      where A = 1 - n0[1+I0]   and   I0=f(s=0)
c
        s = 0.0
        nord = 3
        init_fofs = .true.
        I2 = FOFS(s)
        nord = 1
        init_fofs = .true.
        s = 0.0
        kcc = kc
        I0 = FOFS(s)
        A = 1 - n0*(1+I0)
        write (*,*) ' ############# A = ',A
        ke = A*3/2*1.043*11.6045*a2 + n0*I2
        write (*,*) ' ############# <KE> = ',ke
        maxk = maxk - 0.5*dk
        k = mink
        DO WHILE (k .LE. maxk)
           CALL DQAGI(F4,bound,inf,epa,epr,nk,err,neval,ifail,liw,lw,last,iw,wi)
           nk = A/(2*pi*pi) * nk 
           write (10,'(3G12.4)') k,nk,0.0001
c           if (k.LE.mink) write (*,*) ' k = ',k,'     n(k) = ',nk
           k = k + dk
        END DO
        close (10)
	j = strim(ofile)
	write (*,*) ' '
	write (*,'(A)') ' *** n*(k) stored in '//ofile(1:j)

        write (*,'(T5,A,$)') 'Output f(k)? [N]:'
        read (*,'(A)') dec
        if (dec.EQ.'Y' .OR. dec.EQ.'y') then
           ofile1 = 'fk.dat'
           open (unit=20,file=ofile1,status='unknown')
           dk = 0.01
           k = 0.0
           do i=1,150
              k = k + dk
              write (20,'(3G12.4)') k,n0*FK(k),0.0001
           end do
           k=2.0; write (20,'(3G12.4)') k,n0*FK(k),0.0001
           k=3.0; write (20,'(3G12.4)') k,n0*FK(k),0.0001
           k=4.0; write (20,'(3G12.4)') k,n0*FK(k),0.0001
           close(20)
           j = strim(ofile1)
           write (*,'(A)') ' *** f(k) stored in '//ofile1(1:j)
        end if

        RETURN
        END
c
c
*_____________________________________________________________________________
        REAL FUNCTION F4*8(s)           ! alpha*[n0*beta + n0*beta*f + A*gama]
                                        ! Total J(Q,s)
*_____________________________________________________________________________
	implicit none
        real*8 s,n0,kc,A,a2,a4,a6,ks
        real*8 FOFS,gam,k,arg,SINC
        external SINC
        common /yval/ k
        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
           F4 = EXP(arg)
	   F4 = s**2*F4*SINC(k*s)
        ELSE
           F4 = 0.0
        END IF

        RETURN
        END
c
*--------------------------------------------------------------------   
	REAL FUNCTION FK*8(k) 
*--------------------------------------------------------------------   
	implicit none
	REAL*8 k,pi,B,D,F,sv,COTH
        external COTH
        common /fsc1/ pi,B,D,F,sv
        logical first /.true./
	IF (ABS(k).LE.1.0e-10) THEN
	   Fk=0.0
	ELSE 
	   Fk=B/k*COTH(D*k)*EXP(-k**2/F)
	END IF
	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
        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
        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
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



















































