*______________________________________________________________________________
        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
	character*60 filein,fileres,text*80,dec*1,bkgd_name*80
	integer*4 text_len,lpt,strim,nrd
	logical okay,resfile,moms,firstime,yScaledToMass4
	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 xbkgd(1000),ybkgd(1000),dummy,mass,m4,m3,cen_2
        integer*4 lbkgd,bkgd_len,answer
c
c        common /dats/ y
        common /rdats/ rx,ry,nrd
        common /qval/ qq,mass,cen_2
c        common /logf/ moms
        common /logc/ resfile,yScaledToMass4
c        common /file/ filein,in_len,fileres,res_len
c
	common /bkgd/ xbkgd,ybkgd,lbkgd ! bkgd data

        data nam(1) /'Bkgd Constant'/,
     >       nam(2) /'Bkgd Slope'/,
     >       nam(3) /'Exp Bkgd Scale'/,
     >       nam(4) /'u2 [\A\u-2\d]'/,
     >       nam(5) /'u3 [\A\u-3\d]'/,
     >       nam(6) /'u4 [\A\u-4\d]'/,
     >       nam(7) /'Weight (central peak)'/,
     >       nam(8) /'u2 (2nd peak)'/,
     >       nam(9) /'u3 (2nd peak)'/,
     >       nam(10) /'u4 (2nd peak)'/,
     >       nam(11) /'Weight (2nd peak)'/,
     >       nam(12) /'Shift (2nd peak)'/,
     >       nam(13) /'Mass Ratio (2nd peak)'/
        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 = 13
 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 of data ::'
        read (*,*) qq
        write (*,'(T5,A,$)') ' Indicate y-scale mass used for data 1-He4, 2-He3::'
        read (*,*) answer
        yScaledToMass4 = .false.
        if (answer .eq. 1) yScaledToMass4 = .true.
        
 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
	okay = .false.
	do while (.not.okay)
	   write (*,'(T5,A,$)') 'Enter background data filename::'
	   read (*,'(A)') bkgd_name
	   bkgd_len = strim(bkgd_name)
	   INQUIRE(file=bkgd_name,exist=okay)
	end do
	open (unit=19,file=bkgd_name,status='old')
	do i=1,5000
	   read(19,*,end=29) xbkgd(i),ybkgd(i),dummy
	end do
 29	lbkgd = i-1
	close (19)
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 (*,*) '# ******** He3/4 AA Fitting Program ****** '
        write (*,*) '#'
        write (*,*) '# Type HELP for a list of available commands'
        write (*,*) '#'
        write (*,*) '# Details of fits are stored in /tmp/he34.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
c        p(4) = 1.0*QQ**2        !1.0
c        p(5) = 3.0*QQ**2        !3.75/(1.04257)**3/QQ
c        p(6) = 0.5*QQ**4        !0.5
c        p(7) = 2000*QQ**2       !2500.0/1.04257**5/QQ**3
        p(4) = 1.0              !1.0
        p(5) = 3.0/(1.04257)**3/QQ
        p(6) = 0.5              !0.5
        !p(7) = 2000.0/1.04257**5/QQ**3
        p(7) = Sq*0.5
        p(8) = p(4)
        p(9) = p(5)
        p(10) = p(6)
        mass = 4.0026
        m4 = 4.0026
        m3 = 3.016
        if (yScaledToMass4) then
           cen_2 = 0.5*QQ*(m4 - m3)/m3
           mass = m4
           p(7) = 0.5*Sq
           p(11) = 0.15*Sq
           p(13) = m4/m3
        else
           cen_2 = 0.5*QQ*(m3 - m4)/m4
           mass = m3
           p(11) = 0.5*Sq
           p(7) = 0.15*Sq
           p(13) = m3/m4
        endif
        p(12) = cen_2

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)
        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,m4,m3,cen_2
      real*8 rx(1000),ry(1000)
      real*8 qq,sq,sq_2
      real*8 jia(1000),j1(1000),j2(1000),j3(1000),buf(1000)
      real*8 jia_2(1000),j1_2(1000),j2_2(1000),j3_2(1000),buf_2(1000)
      real*8 y,wd,bc,bkgd_scale,bs,wd_2,u2,u3,u4,u5,u2_2,u3_2,u4_2,u5_2,mass_rat
      real*8 const3,const4,pi
      integer*4 inv(*),nmesh,nrd,i,j,nv,nptot,npk
      real*8 xbkgd(1000),ybkgd(1000)
      integer*4 lbkgd
      logical resfile,yScaledToMass4
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,mass,cen_2
c     common /jqys/ jqy
        common /logc/ resfile,yScaledToMass4
c
	common /bkgd/ xbkgd,ybkgd,lbkgd ! bkgd data
 
	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)
        bkgd_scale = p(3)
	u2 = p(4)
	u3 = p(5)
	u4 = p(6)
        SQ = p(7)
	u2_2 = p(8)
	u3_2 = p(9)
	u4_2 = p(10)
        cen_2 = p(12)
        mass_rat = p(13)
	SQ_2 = p(11) * mass_rat
! 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 + bkgd_scale*ybkgd(j) !Sloping+ (real) exp 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
           wd = Y/SQRT(u2)
           Jia(i)= 1/SQRT(2*pi*u2)*DEXP(-Y**2/(2*u2))
           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) = Jia(i)+J1(i)+J2(i)!+J3(i)
           Y = (Y - cen_2)/mass_rat        ! correction for using wrong mass in y-scaling
           wd_2 = Y/SQRT(u2_2)
           Jia_2(i)= 1/SQRT(2*pi*u2_2)*DEXP(-Y**2/(2*u2_2))
           J1_2(i) = -u3_2*wd_2/(2*U2_2**1.5)*(1 - wd_2**2/3)*Jia_2(i)
           J2_2(i) = u4_2/(8*u2_2**2)*(1 - 2*wd_2**2 + wd_2**4/3)*Jia_2(i)
           !J3_2(i) = u5_2*wd_2/(8*u2_2**2.5)*(1 - 2/3*wd_2**2 + wd_2**4/15)*Jia_2(i)
           buf_2(i) = Jia_2(i)+J1_2(i)+J2_2(i)!+J3_2(i)
           ycal(j) = ypk(j,1)+Sq*buf(i)+Sq_2*buf_2(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,5) = SQ_2*buf_2(i)
           ypk(j,6) = buf(i)    ! save unbroadened calculation (main peak) in comp 6
           ypk(j,7) = buf_2(i)    ! save unbroadened calculation (2nd peak) in comp 7
	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)
           CALL CONVOLVE(nv,x1,buf_2,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,5) = SQ_2*buf_2(i)
              ycal(j) = ypk(j,1)+ypk(j,2)+ypk(j,3)+ypk(j,4)+ypk(j,5)
           END DO
        END IF

        return
        end
