*______________________________________________________________________________
        PROGRAM NPSTAR_FIT
*______________________________________________________________________________
c
c R.T. AZUAH -- April 1999
c    fitting function to PIMC calculations of n*(p)
c
c______________________________________________________________________________
        implicit none
        external npstar
c Variables and common blocks
	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
c        integer*4 nptot,np,inp(50),nd,nrd,i,j,in_len,res_len,mn
	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
c        real*8 p(50),x(1000),y(1000),ysig(1000)
	real*8 p(50),pmin(50),pmax(50),x(1000),y(1000),ysig(1000)
c        real*8 rx(1000),ry(1000),er,qq,Sq,temp,mom,wbar,m(7),ix(1000),iy(1000)






c	logical okay
c
c        common /dats/ nd,x,y,ysig
c        common /pars/ np,inp,p
c        common /nams/ nptot,nam
c        common /labs/ xlab,ylab,xcap,ycap,title,subtitle


c
c
c
c	integer*4 text_len,lpt,strim,nrd,nid
c	logical okay,resfile,moms,comp,firstime
c	integer*4 nptot,nd,i,j,in_len,res_len,mn,fsf_len
c	real*8 p(50),pmin(50),pmax(50),x(1000),y(1000),ysig(1000)
c        real*8 rx(1000),ry(1000),er,qq,Sq,temp,mom,wbar,m(7),ix(1000),iy(1000)
c
c        common /dats/ y
c        common /rdats/ nrd,rx,ry
c        common /qval/ qq,Sq,temp
c        common /logf/ moms
c        common /logc/ firstime,comp,resfile
c        common /file/ filein,in_len,fileres,res_len



c
        data nam(1) /'Bkgd Constant'/,
     >       nam(2) /'Bkgd Slope'/,
     >       nam(3) /'A  '/,
     >       nam(4) /'a2 [(1/A)**2]'/,
     >       nam(5) /'a4 [(1/A)**4]'/,
     >       nam(6) /'a6 [(1/A)**6]'/

        nptot = 6
        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 = 6

 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
c
c Set up labels and captions
c
        title = ' Data file = '//filein(1:in_len)
        xcap = ' |p| (\A) '
        ycap = ' n(|p|) (\A\u3\d) '
c        write (subtitle(1)(13:17), '(i5)') 1
c	subtitle(1)(17:55) = 'Momentum Distr. fits '
c        subtitle(2)(13:20) = 'MARI'
c        subtitle(2)(37:56) = 'RTA/PES/HRG/WGS'
c        write (subtitle(3)(13:17), '(F5.2)') QQ
c        subtitle(3)(37:56) = 'MAY 1992'
c
c Welcome Message
        write (*,*) '#'
        write (*,*) '# ******** n(p) Fitting Program ****** '
        write (*,*) '#'
        write (*,*) '# Type HELP for a list of available commands'
        write (*,*) '#'
        write (*,*) '# Details of fits are stored in /tmp/np.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) = 1.0
        p(4) = 0.90
        p(5) = 0.46
        p(6) = 0.4
c pass control to FRILLS
        open (unit=8,file='/tmp/npstar_fit.lpt',status='unknown')
c
c        CALL FRILLS (npstar)

	firstime = .true.                ! for benefit of fofs when used
        CALL FRILLS (npstar,
	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        SUBROUTINE npstar(ycal,nv,inv,p)
        SUBROUTINE npstar(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)

c General Variables
        real*8 A1,a2,a4,a6
        real*8 jqy,x1(1000)
        integer*4 i,j
c  General common blocks
        common /parm/ a2,a4,a6,A1
        common /jqys/ jqy
c
c Assign fitting parameters from Frills
c
        bc = p(1)
        bs = p(2)
	IF (p(3).NE.0) THEN
	   A1 = 1/p(3)
	ELSE
	   A1 = 1/(p(3)+0.001)
	END IF
        a2 = p(4)
        a4 = p(5)
        a6 = p(6)
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
           CALL nofp(x(j))
           ycal(j) = jqy
        END DO
        npk = 2
        RETURN
        END
c
c
*_____________________________________________________________________________
        SUBROUTINE nofp(Y)
*_____________________________________________________________________________
c
c Program to calculate n*(p) at as a FT of n*(s).
c ____ ie Uncondensed part of n(s) in the superfluid
c
c
        implicit none
        real*8 Qq,Sq,Y,Yv,te
        real*8 A1,a2,a4,a6,pi
        real*8 bound,epa,epr,wi(2000),err,wbar
        real*8 jqy
        real*8 F4
        real D01AMF
        integer*4 lw,liw,iw(500),inf,ifail,i,j,mn,lpt,nord,last,neval
        external F4
c
        common /parm/ a2,a4,a6,A1
        common /jqys/ jqy
        common /yval/ Yv
c
        yv = Y
c
        pi = 4.0*ATAN(1.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


        CALL DQAGI(F4,bound,inf,epa,epr,jqy,err,neval,ifail,liw,lw,last,iw,wi)
        jqy = A1*jqy/(2*pi**2)
c        write (*,*) ' p = ',yv,' np = ',jqy
        RETURN
        END
c
c
*_____________________________________________________________________________
        REAL FUNCTION F4*8(s)           
*_____________________________________________________________________________
	implicit none
        real*8 a2,a4,a6,A1
        real*8 s,arg,Y
        common /parm/ a2,a4,a6,A1
	common /yval/ y
        arg = -a2*s**2/2+a4*s**4/24-a6*s**6/720
	IF (s.EQ.0) THEN
           F4 = 0
        ELSE IF (ABS(arg).LE.60.0) THEN
           F4 = EXP(arg)
c	   F4 = F4*Cos(Y*s)
	   F4 = s**2*F4*Sin(Y*s)/(Y*s)
        ELSE
           F4 = 0.0
        END IF
        RETURN
        END
