*______________________________________________________________________________
        PROGRAM RQY_FIT
*______________________________________________________________________________
c
c # Fits a model R(Q,Y) to experimental data.
c # R(Q,Y) obtained as FT of R(Q,s), the intermediate Final State function
c # R(Q,s) = [ (n0+n0*f(s)+n'(s)) ]*R(Q,s) 
c 
c R.T. AZUAH -- Oct 1998
c
c
c______________________________________________________________________________
	implicit none
	external rqy_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
	logical okay,resfile,moms,comp,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 er,qq,Sq,temp,mom,wbar,m(7)
c
        common /dats/ y
        common /qval/ qq,Sq
        common /logf/ moms
        common /logc/ firstime
        common /file/ filein,in_len,fileres,res_len
c

        data nam(1) /'Bk Const'/,
     >       nam(2) /'Bk Slope'/,
     >       nam(3) /'b3   [\A\u-3\d]'/,
     >       nam(4) /'b4   [\A\u-4\d]'/,
     >       nam(5) /'b5   [\A\u-5\d]'/,
     >       nam(6) /'b6   [\A\u-6\d]'/,
     >       nam(7) /'b7   [\A\u-7\d]'/,
     >       nam(8) /'b8   [\A\u-8\d]'/,
     >       nam(9) /'b9   [\A\u-9\d]'/,
     >       nam(10) /'b10  [\A\u-10\d]'/

        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 = 10
 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 value and S(Q) of data ::'
        read (*,*) qq,Sq
c
c Set up labels and captions
c
        write (title,'(A,F4.1,A)') 'Q=',QQ,'\A\u-1\d'
        xcap = ' Y (\A\u-1\d) '
        ycap = ' R(Q,Y) (\A) '
c        write (subtitle(1)(13:17), '(i5)') 1
	IF (comp) THEN
	   subtitle(1)(17:55) = 'CAN PLOT J(Q,Y) COMPONENTS '
	ELSE
	   subtitle(1)(17:55) = 'CANNOT PLOT J(Q,Y) COMPONENTS'
	END IF
        subtitle(2)(13:20) = 'MARI'
        subtitle(2)(37:56) = 'RTA/HRG'
        write (subtitle(3)(13:17), '(F5.2)') QQ
        subtitle(3)(37:56) = 'APRIL 1998'
c Welcome Message
        write (*,*) '#'
        write (*,*) '# ******** RQY_FIT Fitting Program ****** '
        write (*,*) '#'
        write (*,*) '# Type HELP for a list of available commands'
        write (*,*) '#'
        write (*,*) '# Details of fits are stored in /tmp/rqy_fit.lpt '
        write (*,*) '#'
        write (*,*) '#'
        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) = 3.75/(1.04257)**3/QQ
        p(4) = 0.0
        p(5) = 2500.0/1.04257**5/QQ**3
        p(6) = 200.0/1.04257**6/QQ**2
        p(7) = 0.0;   p(8) = 0.0;   p(9) = 0.0;    p(10) = 0.0
c pass control to FRILLS
        open (unit=8,file='/tmp/rqy_fit.lpt',status='unknown')
	firstime = .true.
c
c        write (*,*) ' JQYS: firstime,comp,resfile = ',firstime,comp,resfile
        CALL FRILLS (rqy_cal,
	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_______________________________________________________________________________
c
*____________________________________________________________________________
        SUBROUTINE RQY_CAL(ycal,nv,inv,nptot,p,x,npk,ypk)
*_____________________________________________________________________________
c
c Determines parameters for R(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(*),ypk(1000,20),ycal(1000),p(*),bc,bs
	integer*4 inv(*),lpt,nd,npk,nv,nptot,nrd
        real*8 ydat(1000)
        logical moms
c Frills common blocks
c        common /rdats/ nrd,rx,ry
        common /dats/ ydat
        common /logf/ moms
c General Variables
        real*8 b3,b4,b5,b6,b7,b8,b9,b10
        real*8 jqy,qq,Sq,temp
        real*8 x1(1000),temp1(1000),temp2(1000)
        integer*4 i,j,nord
        logical firstime
c  General common blocks
        common /parm/ b3,b4,b5,b6,b7,b8,b9,b10
        common /qval/ qq,Sq
        common /ordr/ nord
        common /jqys/ jqy
        common /logc/ firstime

c        write (*,*) ' JQY_CAL: firstime,comp,resfile = ',firstime,comp,resfile

c
c Assign fitting parameters from Frills
c
        bc = p(1)
        bs = p(2)
        b3 = p(3)
        b4 = p(4)
        b5 = p(5)
        b6 = p(6)
        b7 = p(7)
        b8 = p(8)
        b9 = p(9)
        b10 = p(10)
c
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 RQY_TOT(x(j))
           ycal(j) = jqy
           ypk(j,2) = jqy
           npk = 2
           temp1(i) = ydat(j)
           temp2(i) = ycal(j)
        END DO
c
c If requested, calculate and display Experimental and calculated moments
c
        IF (moms) THEN                  
           CALL RQY_MOMS(nv,x1,temp1,temp2)                
        END IF

        RETURN
        END
c
c
*_____________________________________________________________________________
        SUBROUTINE RQY_TOT(Y)
*_____________________________________________________________________________
c
c Program to calculate R(Q,Y) at a given (Q,Y) as a FT of R(Q,s).
c
        implicit none
        real*8 Qq,Sq,Y,Yv
        real*8 b3,b4,b5,b6,b7,b8,b9,b10
        real*8 bound,epa,epr,wi(2000),err,wbar,pi
        real*8 jqy,arg
        real*8 F4,FOFS
        integer*4 lw,liw,iw(500),inf,ifail,i,j,mn,lpt,nord,last,neval
        external F5
        logical firstime,moms,init_fofs
c
        common /parm/ b3,b4,b5,b6,b7,b8,b9,b10
        common /FTan/ lw,liw,ifail,inf,bound,epa,epr
        common /jqys/ jqy
        common /ordr/ nord
        common /qval/ qq,Sq
        common /yval/ Yv
c
        common /logf/ moms
        common /logc/ firstime
	common /ft_fs/ init_fofs
c
        yv = Y
        pi = 4.0*ATAN(1.0)
        IF (firstime) THEN              !setup some initialize parameters
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
        END IF                                  ! end of initialisation
        CALL DQAGI(F5,bound,inf,epa,epr,jqy,err,neval,ifail,liw,lw,last,iw,wi)
        jqy = SQ/pi*jqy
        RETURN
        END
c
c
*_____________________________________________________________________________
        REAL FUNCTION F5*8(s)           ! R(Q,s)
                                        
*_____________________________________________________________________________
	implicit none
        real*8 s,b3,b4,b5,b6,b7,b8,b9,b10
        real*8 p1,arg,Y
        common /parm/ b3,b4,b5,b6,b7,b8,b9,b10
        common /yval/ Y

        p1 = COS(Y*s + b3*s**3/6 - b5*s**5/120 + b7*s**7/5040 - b9*s**9/362880.)
        arg = b4*s**4/24 - b6*s**6/720 + b8*s**8/40320 - b10*s**10/3628800.
        IF (ABS(arg).LE.60.0) THEN
           F5 = EXP(arg)*p1
        ELSE
           F5 = 0.0
        END IF

        RETURN
        END
c
*____________________________________________________________________________
        SUBROUTINE RQY_MOMS(npts,x1,y_expt,ycal_nores)
*_____________________________________________________________________________
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,nord,l1,l2
        real*8 b3,b4,b5,b6,b7,b8,b9,b10
        real*8 x1(1000),y_expt(1000),ycal_nores(1000)
        real*8 M_expt(10),M_ycal_nores(10),M_ycal_res(10),M_exact(10)
        real*8 I0,I2,I4,t1,t2,t3,wbar,mom,error
        real*8 qq,Sq,temp,arg
        character*80 line0,line1,line2,file1*60,file2*60
        logical firstime,moms
c
        common /parm/ b3,b4,b5,b6,b7,b8,b9,b10
        common /ordr/ nord
        common /file/ file1,l1,file2,l2
        common /qval/ qq,Sq
c
        common /logf/ moms
        common /logc/ firstime
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
c        IF (resfile) THEN
c           DO i=1,7
c              mn = i-1
c              wbar = 0.0
c              CALL MOMENTS(x1,ycal_res,npts,mn,mom,error,wbar)
c              M_ycal_res(i) = mom
c           END DO
c        END IF
c
c Evaluate the expected (Exact) moments from the input parameters used to
c calculate the fitting function
c
c        nord = 1
c        arg = 0.0
c        I0 = FOFS(arg)
c        nord = 3
c        I2 = FOFS(arg)
c        nord = 5
c        I4 = FOFS(arg)
c
        M_exact(1) = 1.0                                        ! Zeroth moment
        M_exact(2) = 0.0                                        ! First moment
        M_exact(3) = 0.0                                        ! Second moment
        M_exact(4) = b3                                         ! Third moment
        M_exact(5) = b4                                         ! Fourth moment
        M_exact(6) = b5 + 10.0*b3                               ! Fifth moment
        t1 = b6 + 10*b3**2                                      !** neglect term n0I6/7 **
c        t2 = 15.0*b4*(n0*I2/3 + A*a2)
c        t3 = 15.0*A*a2*(a2**2 + a4)
        M_exact(7) = t1 !+t2+t3                                   ! Sixth moment 
c
        write (*,*) '#_____________________________________________________________________________'
        write (*,'(A)') ' # ** Data file = '//file1(1:l1)
c        write (*,'(A)') ' # ** Resn file = '//file2(1:l2)
        write (*,*) '# '
        write (*,'(2(A,F6.3))') ' # ** Q = ',QQ,'  ** S(Q) = ',Sq  
c	write (*,'(3(A,F6.3))') ' # ** n0= ',n0,'  **  Io  = ',I0,'  **    A  = ',A 
        write (*,*) '# '
        line0 = '#****** Moments evaluated are as follows  : '
        write (*,*) '# '
        line1 = '                  DATA        CALCULATION          EXACT   '
        line2 = '                  ----        -----------          -----   '
        write (*,'(A)') line0
        write (*,'(A)') line1
        write (*,'(A)') line2
        write (*,20) 'M0 = ',M_expt(1),M_ycal_nores(1),M_exact(1)
        write (*,20) 'M1 = ',M_expt(2),M_ycal_nores(2),M_exact(2)
        write (*,20) 'M2 = ',M_expt(3),M_ycal_nores(3),M_exact(3)
        write (*,20) 'M3 = ',M_expt(4),M_ycal_nores(4),M_exact(4)
        write (*,20) 'M4 = ',M_expt(5),M_ycal_nores(5),M_exact(5)
        write (*,20) 'M5 = ',M_expt(6),M_ycal_nores(6),M_exact(6)
        write (*,20) 'M6 = ',M_expt(7),M_ycal_nores(7),M_exact(7)
 20     FORMAT (T5,A5,5X,G12.4,3X,G12.4,G12.4)
        write (*,*) '# '
        write (*,*) '#_____________________________________________________________________________'
        RETURN
        END

c
































