*_____________________________________________________________________________
        PROGRAM NOFS
*_____________________________________________________________________________
c
c Calculate n(s) for superfluid He4
c
c R.T. Azuah  -- March 1998
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,val
        real*8 I0,const1,const2,const3,const4,const5
        real*8 F4,FOFS,mins,maxs,s,ds,arg
        integer*4 i,j,nord
        external F4,FOFS
        logical firstime,comp,resfile,moms,init_fofs
        character*60 ofile
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 /ordr/ nord


	common /ft_fs/ init_fofs
c
c Prompt and get necessary input
c
	write (*,'(T5,A,$)') 'Enter a2 (1/A^2), a4 (1/A^4) and b6 (1/A^6) :'
	read (*,*) a2,a4,a6
c        a2 = 0.897; a4 = 0.46;  a6 = 0.38
        write (*,'(T5,A,$)') 'Enter n0: '
        read (*,*) n0
        kc=0.5; temp=1.6
c	write (*,'(T5,A,$)') 'Enter min, max and interval for s: '
c	read (*,*) mins,maxs,ds
        mins = 0.0; maxs=10.0; ds=0.1
	write (*,'(T5,A,$)') 'Give output data file :'
	read (*,'(A)') ofile 
c        ofile = 'ns.dat'
	j = (maxs-mins)/ds+1
	open (unit=10,file=ofile,status='unknown')
	INQUIRE (unit=10,name=ofile)

c        write (*,*) 'a6,temp,ds = ',a6,temp,ds
c
c  Define some fundamental constants
c
        qq = 20.0
        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   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
c        I0 = FOFS(arg)
        I0 = 0.0
        A = 1 - n0*(1+I0)
        kcc = kc
        A = 1 - n0*(1+I0)
c        write (*,*) 'I0,A = ',I0,A

c
c Calculate n(s)
c
c	write (*,*) ' j = ',j
        DO i = 1,j
	   s = mins + ds*(i-1)
           nord = 1
c           write (*,*) ' '
	   write (10,'(3G12.4)') s,F4(s),0.001
        END DO
	call strtrim(ofile,ofile,j)
	write (*,*) ' '
	write (*,'(A)') ' Output stored in '//ofile(1:j)

        STOP
        END
c
*_____________________________________________________________________________
        REAL*8 FUNCTION F4(s)           ! alpha*[n0*beta + n0*beta*f + A*gama]
                                        ! Total J(Q,s)
                                        ! Here, alpha(s)=beta(s) = 1 
*_____________________________________________________________________________
	implicit none
        real*8 s,n0,kc,A,a2,a4,a6
        real*8 GAMA,FOFS,t1,t2
        external GAMA,FOFS
        common /parm/ n0,kc,A,a2,a4,a6
c        t1 = FOFS(s)
        t2 = GAMA(s)
        F4 = n0*(1.0 + t1) + A*t2
c        F4 = n0 + n0*FOFS(s) + A*GAMA(s)
c        write (*,*) 's,fofs,gama,f4 = ',s,t1,t2,f4
        RETURN
        END
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
        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
