*______________________________________________________________________________
        PROGRAM JAA_FIT
*______________________________________________________________________________

c # Fits Additive Approach model J(Q,Y) to experimental data.
c # J(Q,y) = JIA + J1 + J2 + ...
c # where  * JIA is the (gaussian) IA component of model
c          * J1,J2, etc are additional corrections
c  
c R.T. AZUAH -- 2002
c
c______________________________________________________________________________
        implicit none
        external jaa_cal
c
	character xlab*5,ylab*5,xcap*60,ycap*60,subtitle(3)*60,title*80,nam(50)*20,ans*1
	character*60 filein,fileres,text*80,dec*1,f1,f2,f3
	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)
        real*8 ycal(1000),ypk(5000,20)
        integer*4 nv,npk,inv(1000)
c
c        common /dats/ y
        common /rdats/ rx,ry,nrd
        common /qval/ qq
c        common /logf/ moms
        common /logc/ resfile
c        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) /'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 = 7
        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)
              inv(i) = i
           END DO
 15        nd = i - 1   
           nv = nd
           close (18)
        ELSE
           write (*,*) '    ** ERROR - Data file does not exist '
           GOTO 10
        END IF
        write (*,'(T5,A,$)') ' Enter Q value of data ::'
        read (*,*) qq
 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)') 'Q=',QQ,'\A\u-1\d'
        xcap = ' Y (\A\u-1\d) '
        ycap = ' J(Q,Y) (\A) '
	subtitle(1)(17:55) = 'AA 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 (*,*) '# ******** JAA Fitting Program ****** '
        write (*,*) '#'
        write (*,*) '# Type HELP for a list of available commands'
        write (*,*) '#'
        write (*,*) '# Details of fits are stored in /tmp/jaa.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*QQ**2        !1.0
        p(4) = 3.0*QQ**2        !3.75/(1.04257)**3/QQ
        p(5) = 0.5*QQ**4        !0.5
        p(6) = 2000*QQ**2       !2500.0/1.04257**5/QQ**3
        p(7) = Sq
c pass control to FRILLS
        open (unit=8,file='/tmp/jaa.lpt',status='unknown')
        firstime = .true.
        CALL FRILLS (jaa_cal,
     $       nd,x,y,ysig,
     $       nptot,p,pmin,pmax,nam,
     $       xlab,ylab,xcap,ycap,subtitle,title)
c
c Save function components to files
c
 10     write (*,'(T5,A,$)') 'Save fitted AA function components (y or n)?::'
        read (*,'(A)') ans
        if ((ans .eq. 'y') .or. (ans .eq. 'Y')) then 
           write (*,*) p(3),p(4),p(5)
           call jaa_cal(ycal,nv,inv,nptot,p,x,npk,ypk)
           f1 = 'jg.dat'
           f2 = 'j1.dat'
           f3 = 'j2.dat'
           open (unit=31,file=f1,status='unknown')
           open (unit=32,file=f2,status='unknown')
           open (unit=33,file=f3,status='unknown')
           do i=1,nd
              write (31,'(3G11.4,X)') x(i),ypk(i,2),0.05*ypk(i,2)
              write (32,'(3G11.4,X)') x(i),ypk(i,3),0.05*ypk(i,3)
              write (33,'(3G11.4,X)') x(i),ypk(i,4),0.05*ypk(i,4)
           end do
           
           close(31)
           close(32)
           close(33)

           write (*,*) 'Gaussian (u2) component saved in "jg.dat"'
           write (*,*) 'u3 component saved in "j1.dat"'
           write (*,*) 'u4 component saved in "j2.dat"'
        end if
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
c===============================================================================
c===============================================================================
      SUBROUTINE jaa_cal(ycal,nv,inv,nptot,p,x,npk,ypk)


c
c Additive Approach Model function. Here J(Q,Y) is represented by additive
c terms, each representing a particular aspect of the scattering function.
c
      implicit none
      real*8 x(*),ycal(*),ypk(5000,20),p(*)
      real*8 x1(1000),mass
      real*8 rx(1000),ry(1000)
      real*8 qq,sq
      real*8 jia(1000),j1(1000),j2(1000),j3(1000),buf(1000)
      real*8 y,wd,bc,bs,u2,u3,u4,u5,u6
      real*8 const3,const4,pi
      integer*4 inv(*),nmesh,nrd,i,j,nv,nptot,npk
      logical resfile
c Frills common blocks
c	common /logf/ moms
c  General common blocks
        common /rdats/ rx,ry,nrd
c        common /dats/ ydat
c        common /parm/ u2,u3,u4,u5,u6
        common /qval/ qq
c        common /jqys/ jqy
        common /logc/ resfile


        mass = 4.0026
	const3=2.072194*1.008665*2.0
	const4=mass/(QQ*const3) ! Y=const4*(w-wr) ; Jqy=Sqw/const4
	PI=4*ATAN(1.0)

! Definition of parameters
	bc = p(1)
	bs = p(2)
	u2 = p(3)
	u3 = p(4)
	u4 = p(5)
	u5 = p(6)
        SQ = p(7)
! Loop over data points - only those points included in INV to be fitted
	do i = 1,nv
           j = inv(i)
           x1(j)=x(j)
           y = x(j)
           ypk(j,1) = bc + bs * y !Sloping background
           !u2= U2*QQ**2
           !u3= U3*QQ**2
           !u4= U4*QQ**4
           !u5= u5*QQ**2
           wd = (Y/const4)/SQRT(u2)
           Jia(i)= 1/SQRT(2*pi*u2)*EXP(-(Y/const4)**2/(2*u2))/const4
           J1(i) = -u3*wd/(2*U2**1.5)*(1 - wd**2/3)*Jia(i)
           J2(i) = u4/(8*u2**2)*(1 - 2*wd**2 + wd**4/3)*Jia(i)
           J3(i) = u5*wd/(8*u2**2.5)*(1 - 2/3*wd**2 + wd**4/15)*Jia(i)
           buf(i) = SQ*(Jia(i)+J1(i)+J2(i)+J3(i))
           ycal(j) = ypk(j,1)+buf(i)
	END DO
        
        
        IF (resfile) THEN
           CALL CONVOLVE(nv,x1,jia,nrd,rx,ry)
           CALL CONVOLVE(nv,x1,j1,nrd,rx,ry)
           CALL CONVOLVE(nv,x1,j2,nrd,rx,ry)
           CALL CONVOLVE(nv,x1,j3,nrd,rx,ry)
           DO i = 1,nv
              j = inv(i)
              ypk(j,2) = SQ*jia(i)
              ypk(j,3) = SQ*j1(i)
              ypk(j,4) = SQ*j2(i)
              ypk(j,5) = SQ*j3(i)
              ypk(j,6) = buf(i) ! save unbroadened calculation in comp 6 
              ycal(j) = ypk(j,1)+ypk(j,2)+ypk(j,3)+ypk(j,4)+ypk(j,5)
           END DO
        END IF

        return
        end
