*______________________________________________________________________________
        PROGRAM JQY_FIT
*______________________________________________________________________________

c # Fits a model J(Q,Y) to experimental data.
c # J(Q,Y) obtained as FT of J(Q,s), the intermediate scattering function
c # J(Q,s) = J(s)*R(Q,s)
c # where  * J(s) is the OBDM in s-space
c          * R(Q,s) = Final State function
c  
c R.T. AZUAH -- Aug 1996
c
c______________________________________________________________________________
        implicit none
        external jqy_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
	integer*4 text_len,lpt,strim,nrd
	logical okay,resfile,moms,firstime
	integer*4 nptot,nd,i,j,in_len,res_len,mn
	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)
c
        common /dats/ y
        common /rdats/ nrd,rx,ry
        common /qval/ qq,Sq,temp
        common /logf/ moms
        common /logc/ resfile
        common /file/ filein,in_len,fileres,res_len
c
        data nam(1) /'Bkgd Constant'/,
     >       nam(2) /'Bkgd Slope'/,
     >       nam(3) /'u2 [\A\u-2\d]'/,
     >       nam(4) /'u3 [\A\u-3\d]'/,
     >       nam(5) /'u4 [\A\u-4\d]'/,
     >       nam(6) /'u5 [\A\u-5\d]'/,
     >       nam(7) /'u6 [\A\u-6\d]'/,
     >       nam(8) /'Weight'/
        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 = 8
 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 and Temp of data ::'
        read (*,*) qq,temp
 20     write (*,'(T5,A,$)') ' Enter corresponding Instr. 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) rx(i),ry(i),er
           END DO
 25        nrd = i - 1  
           close (18)
        ELSE
           write (*,*) '    ** ERROR - Instr Resolution file does not exist '
           GOTO 20
        END IF

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) '
	subtitle(1)(17:55) = 'Convolution Model Fit'
        subtitle(2)(13:20) = 'MARI'
        subtitle(2)(37:56) = 'RTA/HRG'
        write (subtitle(3)(13:17), '(F5.2)') QQ
        subtitle(3)(37:56) = 'Feb 2000'
c Welcome Message
        write (*,*) '#'
        write (*,*) '# ******** JQYS Fitting Program ****** '
        write (*,*) '#'
        write (*,*) '# Type HELP for a list of available commands'
        write (*,*) '#'
        write (*,*) '# Details of fits are stored in SYS$SCRATCH:JQYSN.LPT'
        write (*,*) '#'
        write (*,*) '#'
        call area (x,y,ysig,nd,Sq) 	! determine data area -- S(Q)
        p(1) = 0.0			! estimate initial fitting parameters
        p(2) = 0.0
        p(3) = 1.0
        p(4) = 3.75/(1.04257)**3/QQ
        p(5) = 0.5
        p(6) = 2500.0/1.04257**5/QQ**3
        p(7) = 0.4
        p(8) = Sq
c pass control to FRILLS
        open (unit=8,file='/tmp/jqy.lpt',status='unknown')
        firstime = .true.
        CALL FRILLS (jqy_cal,
     $       nd,x,y,ysig,
     $       nptot,p,pmin,pmax,nam,
     $       xlab,ylab,xcap,ycap,subtitle,title)
c
c Save fitted parameters to a file
c
c	filein = 'parameters.fit'
c	INQUIRE(file=filein,exist=okay)
c	IF (okay) THEN
c	   open (unit=27,file=filein,status='old',access='append')
c	ELSE
c	   open (unit=27,file=filein,status='new',carriagecontrol='list')
c	   text = '  Q     U2     dU2      U3     dU3     U4     dU4    U5     dU5     U6    dU6 '
c	   write(27,'(A)') text
c	END IF
c	write(27,'(1X,F4.1,5(2X,F6.3,1X,F6.4))') qq,(p(i),pfsig(i),i=3,7)
c	close (27)
        close (8)			! close fitting log file
        STOP
        END

*____________________________________________________________________________
        SUBROUTINE JQY_CAL(ycal,nv,inv,nptot,p,x,npk,ypk)
*_____________________________________________________________________________
c
c Calculates J(Q,Y) 
c
        implicit none
c Frills Variables 
        real*8 x(1000),ypk(5000,20),ycal(1000),p(50)
        integer*4 nv,inv(1000),nd,npk,nrd,nptot
	logical moms
c Frills common blocks
	common /logf/ moms
c General Variables
        real*8 rx(1000),ry(1000),ydat(1000),bc,bs
        real*8 u2,u3,u4,u5,u6
        real*8 jqy,qq,Sq,temp
        real*8 x1(1000),temp1(1000),temp2(1000),temp3(1000)
        integer*4 i,j
        logical resfile
c  General common blocks
        common /rdats/ nrd,rx,ry
        common /dats/ ydat
        common /parm/ u2,u3,u4,u5,u6
        common /qval/ qq,Sq,temp
        common /jqys/ jqy
        common /logc/ resfile
c
c Assign fitting parameters from Frills
c
        bc = p(1)
        bs = p(2)
        u2 = p(3)
        u3 = p(4)
        u4 = p(5)
        u5 = p(6)
        u6 = p(7)
        SQ = p(8)
c
c Calculate Function plus sloping background
c __ Calculation is made only at data points
c
        DO i = 1,nv
           j = inv(i)
           x1(i) = x(j)
           ypk(j,1) = bc + bs*x(j)      ! Sloping background
           temp3(i) = ydat(j)
           CALL JQY_TOT(x(j))
           ycal(j) = jqy
           ypk(j,2) = jqy
           temp1(i) = jqy
           npk = 2
           temp2(i) = ycal(j)
        END DO
c
c Convolute Calculated function with Instrumental Resolution if present 
c
        IF (resfile) THEN
           CALL CONVOLVE(nv,x1,temp1,nrd,rx,ry)
           DO i = 1,nv
	      j = inv(i)
              ypk(j,2) = temp1(i)
              ycal(j) = temp1(i) + ypk(j,1)
           END DO
        END IF
c
c If requested, calculate and display Experimental and calculated moments
c
        IF (moms) THEN                  
           CALL JQY_MOMS(nv,x1,temp3,temp2,temp1)                
        END IF
        RETURN
        END
c
c
*_____________________________________________________________________________
        SUBROUTINE JQY_TOT(Y)
*_____________________________________________________________________________
c
c Program to calculate J(Q,Y) at a given Q and Y as a FT of J(Q,s).
c
c
        implicit none
        real*8 Qq,Sq,Y,Yv
        real*8 u2,u4,u3,u5,u6
        real*8 bound,epa,epr,wi(2000),err,wbar
        real*8 jqy,pi,temp
        real*8 JOFS
        integer*4 lw,liw,iw(500),inf,ifail,i,j,mn,lpt,nord,last,neval
        external JOFS
        logical resfile
c
        common /parm/ u2,u3,u4,u5,u6
        common /jqys/ jqy
        common /qval/ qq,Sq,temp
        common /yval/ Yv
        common /logc/ resfile
c
        yv = Y
	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
c
c calculate J(Q,Y)
c
        !CALL D01AMF(JOFS,bound,inf,epa,epr,jqy,err,wi,lw,iw,liw,ifail)
        CALL DQAGI(JOFS,bound,inf,epa,epr,jqy,err,neval,ifail,liw,lw,last,iw,wi)
        jqy = SQ/pi*jqy
   
        RETURN
        END
c
*_____________________________________________________________________________
        REAL*8 FUNCTION JOFS(s)           
*_____________________________________________________________________________
	implicit none
        real*8 s,u2,u3,u4,u5,u6,y,arg
        common /parm/ u2,u3,u4,u5,u6
        common /yval/ Y
        arg = -u2*s**2/2+u4*s**4/24-u6*s**6/720
        IF (ABS(arg).LE.60.0) THEN
           JOFS = EXP(arg)*COS(Y*s+u3*s**3/6-u5*s**5/120)
        ELSE
           JOFS = 0.0
        END IF
	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,l1,l2
        real*8 u2,u3,u4,u5,u6
        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 wbar,mom,error
        real*8 qq,Sq,temp
        character*80 line0,line1,line2,file1*60,file2*60
        logical resfile,moms
c
        common /parm/ u2,u3,u4,u5,u6
        common /file/ file1,l1,file2,l2
        common /qval/ qq,Sq,temp
c
        common /logf/ moms
        common /logc/ resfile
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
        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
        M_exact(1) = 1.0                                        ! Zeroth moment
        M_exact(2) = 0.0                                        ! First moment
        M_exact(3) = u2			                        ! Second moment
        M_exact(4) = u3                                         ! Third moment
        M_exact(5) = 3*u2**2+u4			                ! Fourth moment
        M_exact(6) = u5 + 10.0*u2*u3				! Fifth moment
        M_exact(7) = 15*(u2**3+u2*u4)+10*u3**2+u6               ! Sixth moment
c
        write (*,*) '#_____________________________________________________________________________'
        write (*,'(A)') ' # ** Data file = '//file1(1:l1)
        write (*,'(A)') ' # ** Resn file = '//file2(1:l2)
        write (*,*) '# '
        write (*,'(2(A,F6.3))') ' # ** Q = ',QQ,'  ** S(Q) = ',Sq    
        write (*,'(A,F6.3)') '  **  KE  = ',3/2.0*1.043*11.6045*u2
        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 (*,*) '#_____________________________________________________________________________'
        RETURN
        END
