      subroutine JY_AA(ycal, nv, inv, p)
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
      real*8 x,y,ysig, ypk, ycal(1000), p(50), h, w, c, bc, bs
      real*8 x1(1000),y1(1000),e1(1000),mass
      real*8 rx(1000),ry(1000),rysig(1000),sigi,wfac
      real*8 lot,hit,low,hiw,sig2,c4,qq,sq,lam2,c6,eps
      real*8 xtable(1000),ytable(1000),ztable(1000)
      real*8 Sia(1000),S1(1000),S2(1000),S3(1000)
      real*8 u2,u3,u4,u5,u6
      integer*4 inv(1000),nmesh,cal_rqt
      logical resfile
      common /dats/ nd, x(1000), y(1000), ysig(1000)
      common /peak/npk, ypk(1000,10)
      common /rdats/ nrd,rx,ry,rysig
      common /new1/ nmesh,lot,hit,low,hiw,sig2,c4,qq,sq,lam2,c6,eps,cal_rqt,mass
      common /new2/ u2,u3,u4,u5,bc,bs
      common /new3/ xtable,Sia,S1,S2,S3
      common /new5/ sigi,wfac,resfile
! Definition of parameters
	bc = p(1)
	bs = p(2)
	u2 = p(3)
	u3 = p(4)
	u4 = p(5)
	u5 = p(6)
! Loop over data points - only those points included in INV to be fitted
	do i = 1,nv
          j = inv(i)
	  x1(j)=x(j)
          ypk(j,1) = bc + bs * x1(j)	!Sloping background
	END DO
c
c Calculate J(Q,Y) using routine AAF1
c  returns frequency, real and imaginary part of J(Q,Y) in common block /new3/
c  need common blocks /new1/, /new2/ and /new3/
c
	CALL AAF1
	DO i=1,nmesh
	   ztable(i)=SQRT(ABS(Sia(i)))
	END DO
c
c rebin J(Q,Y) to obtain values at the points where exper't data are found
c
	CALL rebin_data(nmesh,xtable,Sia,ztable,nv,x1,y1,e1)
	DO i=1,nv
	   j=inv(i)
	   ypk(j,2)=y1(j)
	END DO
	CALL rebin_data(nmesh,xtable,S1,ztable,nv,x1,y1,e1)
	DO i=1,nv
	   j=inv(i)
	   ypk(j,3)=y1(j)
	END DO
	CALL rebin_data(nmesh,xtable,S2,ztable,nv,x1,y1,e1)
	DO i=1,nv
	   j=inv(i)
	   ypk(j,4)=y1(j)
	END DO
	CALL rebin_data(nmesh,xtable,S3,ztable,nv,x1,y1,e1)
	DO i=1,nv
	  j=inv(i)
	   ypk(j,5)=y1(j)
	END DO
	DO i=1,nv
	  j=inv(i)
          ycal(j) = ypk(j,1)+ypk(j,2)+ypk(j,3)+ypk(j,4)+ypk(j,5)
	end do
	npk = 5
	return
	end 
      SUBROUTINE AAF1
**************************************************************************      
* 	R.T. AZUAH OCTOBER 1992
**************************************************************************
      IMPLICIT NONE                                                     
      INTEGER NMESH,I,TYPE,J,mn
      INTEGER MAXDIM,IERR,SKIP,DATUNT,DATMAX,cal_rqt
      PARAMETER(MAXDIM=1000)         
      real*8 A,EPS,LOW,HIW,ERR,XERR,Y,SIGV,PI,LOT,HIT,mass
      real*8 XTABLE(1000),sqe(1000),step,ytable(1000),ztable(1000)
      real*8 WORK(2*MAXDIM+2),YMOM(0:MAXDIM),MOMENT(0:10)
      real*8 EXACT(0:10),EXERR(0:10),TT(0:6),LAM2,OTT(0:10)
      real*8 SQ,WR,wbar,U2,U3,U4,U5,U6,C4,C6,SIG2,WFAC,QQ,bc,bs
      real*8 rtab1(1000),rtab2(1000),rtab3(1000),rtab4(1000)
      real*8 Ri(1000),dtab(1000),wp,Y_Yp,wd,wpr
      real*8 Sia(maxdim),S1(maxdim),S2(maxdim),S3(maxdim)
      real*8 obsx(1000),obsy1(1000),obsy2(1000),obsy3(1000),obsy4(1000)
      real*8 u22,u33,u44,u55,u66,w,mom,error
      real*8 rx(1000),ry(1000),rysig(1000),sigi
      character dec,sqw_out*50,sia_out*50,s1_out*50,s2_out*50,s3_out*50,mom_out*50
      integer dec_sqw,dec_sia,dec_s1,dec_s2,dec_s3,dec_mom,nrd
      REAL*8 DPCHQA,const3,const4,EXTERNAL DPCHQA
      LOGICAL SPLINE,DATFLG,CONFLG,EXPFLG,resfile
      PARAMETER (SPLINE=.FALSE.)                     
      COMMON /PARM/ Y,TT,WR
      common /rdats/ nrd,rx,ry,rysig
      common /new1/ nmesh,lot,hit,low,hiw,sig2,c4,qq,sq,lam2,c6,eps,cal_rqt,mass
      common /new2/ u2,u3,u4,u5,bc,bs
      common /new3/ xtable,Sia,S1,S2,S3
      common /new5/ sigi,wfac,resfile
C
	const3=2.072194*1.008665*2.0
	const4=mass/(QQ*const3)			! Y=const4*(w-wr) ; J=S/const4
	PI=4*ATAN(1.0)
	lam2=sigi
C
	IF (cal_rqt.EQ.1) THEN
	  write (*,'(T5,A,$)') '######### Output J(Q,Y)? def[n] :: '
	  read (*,'(a)')  dec
	  IF (dec.EQ.'y') THEN
	    dec_sqw=1
	    OPEN(unit=20,file='jqy.dat',status='new',carriagecontrol='list')
	    inquire (unit=20,name=sqw_out)
	    write (*,'(T5,A)') 'J(Q,Y) output stored in '//sqw_out
	  END IF
	  write (*,'(T5,A,$)') '######### Output JIA(Q,Y)? def[n] :: '
	  read (*,'(a)')  dec
	  IF (dec.EQ.'y') THEN
	    dec_sia=1
	    OPEN(unit=40,file='JIA.dat',status='new',carriagecontrol='list')
	    inquire (unit=40,name=sia_out)
	    write (*,'(T5,A)') 'Jia(Q,Y) output stored in '//sia_out
	  END IF
	  write (*,'(T5,A,$)') '######### Output J1(Q,Y)? def[n] :: '
	  read (*,'(a)')  dec
	  IF (dec.EQ.'y') THEN
	    dec_s1=1
	    OPEN(unit=60,file='J1.dat',status='new',carriagecontrol='list')
	    inquire (unit=60,name=s1_out)
	    write (*,'(T5,A)') 'J1(Q,Y) output stored in '//s1_out
	  END IF
	  write (*,'(T5,A,$)') '######### Output J2(Q,Y)? def[n] :: '
	  read (*,'(a)')  dec
	  IF (dec.EQ.'y') THEN
	    dec_s2=1
	    OPEN(unit=90,file='J2.dat',status='new',carriagecontrol='list')
	    inquire (unit=90,name=s2_out)
	    write (*,'(T5,A)') 'J2(Q,Y) output stored in '//s2_out
	  END IF
	  write (*,'(T5,A,$)') '######### Output J3(Q,Y)? def[n] :: '
	  read (*,'(a)')  dec
	  IF (dec.EQ.'y') THEN
	    dec_s3=1
	    OPEN(unit=95,file='J3.dat',status='new',carriagecontrol='list')
	    inquire (unit=95,name=s3_out)
	    write (*,'(T5,A)') 'J3(Q,Y) output stored in '//s3_out
	  END IF
	  write (*,'(T5,A,$)') '######### Output MOMENTS? def[n] :: '
	  read (*,'(a)')  dec
	  IF (dec.EQ.'y') THEN
	    dec_mom=1
	    OPEN(unit=70,file='mom.dat',status='new',carriagecontrol='list')
	    inquire (unit=70,name=mom_out)
	    write (*,'(T5,A)') 'Moments  stored in '//mom_out
	  END IF
	END IF
c	OPEN(unit=10,file='sqt.dat',status='new',carriagecontrol='list')
	WR=QQ**2*WFAC         
	wpr=wr/sq
c
c Calculate Sia, S1 and S2 arrays
c
c
	IF (resfile) THEN
	   u22= U2*QQ**2
	ELSE
	   u22= u2*QQ**2+lam2**2
	END IF
	u33= U3*QQ**2
	u44= U4*QQ**4
	u55= u5*QQ**2
	step=(hiw-low)/(nmesh-1)
	DO i=1,nmesh
	  Y = low+(i-1)*step
	  xtable(i)=Y
	  wd = (Y/const4)/SQRT(u22)
	  Sia(i)= Sq/SQRT(2*pi*u22)*EXP(-(Y/const4)**2/(2*u22))/const4
	  S1(i) = -u33*wd/(2*U22**1.5)*(1 - wd**2/3)*Sia(i)
	  S2(i) = u44/(8*u22**2)*(1 - 2*wd**2 + wd**4/3)*Sia(i)
	  S3(i) = u55*wd/(8*u22**2.5)*(1 - 2/3*wd**2 + wd**4/15)*Sia(i)
	  sqe(i) = sia(i)+s1(i)+s2(i)  !+s3(i)
	END DO
c
c Convolute with instrument resolution function if required
c
	IF (resfile) THEN
 	   DO i=1,nmesh
	      Y = low+(i-1)*step
	      xtable(i)=Y
		DO j=1,nrd
		   Y_Yp = xtable(i)-rx(j)
		   wd = (Y_Yp/const4)/SQRT(u22)
	  	   Sia(j)= Sq/SQRT(2*pi*u22)*EXP(-(Y_Yp/const4)**2/(2*u22))/const4
	  	   S1(j) = -u33*wd/(2*U22**1.5)*(1 - wd**2/3)*Sia(j)
	  	   S2(j) = u44/(8*u22**2)*(1 - 2*wd**2 + wd**4/3)*Sia(j)
	  	   S3(j) = u55*wd/(8*u22**2.5)*(1 - 2/3*wd**2 + wd**4/15)*Sia(j)
		   obsy1(j)=Sia(j)*ry(j)
		   obsy2(j)=S1(j)*ry(j)
		   obsy3(j)=S2(j)*ry(j)
		   obsy4(j)=S3(j)*ry(j)
		END DO
		wbar=0
		mn=0
		CALL MOMENTS(rx,obsy1,nrd,mn,mom,error,wbar)
		rtab1(i)=mom
		CALL MOMENTS(rx,obsy2,nrd,mn,mom,error,wbar)
		rtab2(i)=mom
		CALL MOMENTS(rx,obsy3,nrd,mn,mom,error,wbar)
		rtab3(i)=mom
		CALL MOMENTS(rx,obsy4,nrd,mn,mom,error,wbar)
		rtab4(i)=mom

	   END DO

	   DO i=1,nmesh
	      Sia(i) = rtab1(i)
	      S1(i) = rtab2(i)
	      S2(i) = rtab3(i)
	      S3(i) = rtab4(i)
	      sqe(i) = sia(i)+s1(i)+s2(i)  !+s3(i)
	   END DO
	END IF
c
c Output exact moments to output file and compare with calculated 
c moments from FT data of J(Q,Y)
c
       IF (dec_mom.EQ.1) THEN
	write (70,'(T5,A,G10.4)') 'AA - Moments check for Q = ',QQ
	write (70,'(T5,A,G10.4)') 'Sq      = ',Sq
	write (70,'(T5,A,G10.4)') 'Sigmai  = ',lam2
	write (70,'(T5,A,G10.4)') 'Bconst  = ',bc
	write (70,'(T5,A,G10.4)') 'Bslope  = ',bs
	write (70,'(T5,A,G10.4)') 'U2      = ',u2
	write (70,'(T5,A,G10.4)') 'U3      = ',u3
	write (70,'(T5,A,G10.4)') 'U4      = ',u4
	write (70,'(T5,A,G10.4)') 'U5      = ',u5
	write (70,'(T5,A,G10.4)') 'Min Y   = ',low
	write (70,'(T5,A,G10.4)')  'Max Y   = ',hiw
	write (70,'(T5,A,I6)') 'No of points = ',nmesh
	write (70,*) 
	write (70,'(T5,A)') '######## Exact moments of J(Q,Y) ######'
	write (70,*) 
	IF (resfile) THEN
	   DO i=0,4    !5
	     mn=i
	     wbar=0.0
	     call moments(rx,ry,nrd,mn,mom,error,wbar)
	     exact(i)=mom
	     exerr(i)=abs(error)
	   END DO
	   exact(0)=sq*exact(0)
	   exact(1)=0.0+exact(1)
	   exact(2)=sq*u22*const4**2+exact(2)
	   exact(3)=sq*u33*const4**3+exact(3)
	   exact(4)=sq*(u44+3*u22**2)*const4**4+exact(4)
	   exact(5)=sq*(u55+10*((u22+lam2**2)*u33))*const4**4+exact(5)
	ELSE
	   exact(0)=sq	
	   exact(1)=0.0
	   exact(2)=sq*u22*const4**2
	   exact(3)=sq*u33*const4**3
	   exact(4)=sq*(u44+3*u22**2)*const4**4
	   exact(5)=sq*(u55+10*((u22+lam2**2)*u33))
	END IF
	write (70,'(T5,A,G10.4)') 'ZEROTH CENTRAL MOMENT  :: ',exact(0)
	write (70,'(T5,A,G10.4)') 'FIRST  CENTRAL MOMENT  :: ',exact(1)
	write (70,'(T5,A,G10.4)') 'SECOND CENTRAL MOMENT  :: ',exact(2)
	write (70,'(T5,A,G10.4)') 'THIRD  CENTRAL MOMENT  :: ',exact(3)
	write (70,'(T5,A,G10.4)') 'FOURTH CENTRAL MOMENT  :: ',exact(4)
	write (70,'(T5,A,G10.4)') 'FIFTH  CENTRAL MOMENT  :: ',exact(5)
	write (70,*) 
	DO i=0,4    !5
	  mn=i
	  wbar=0.0
	  call moments(xtable,sqe,nmesh,mn,mom,error,wbar)
	  exact(i)=mom
	  exerr(i)=abs(error)
	END DO
	write (70,'(T5,A)') '######## Calculated moments of J(Q,Y) ######'
	write (70,*) 
	write (70,'(T5,2(A,G10.4))') 'ZEROTH CENTRAL MOMENT  :: ',exact(0),' +/- ',exerr(0)
	write (70,'(T5,2(A,G10.4))') 'FIRST  CENTRAL MOMENT  :: ',exact(1),' +/- ',exerr(1)
	write (70,'(T5,2(A,G10.4))') 'SECOND CENTRAL MOMENT  :: ',exact(2),' +/- ',exerr(2)
	write (70,'(T5,2(A,G10.4))') 'THIRD  CENTRAL MOMENT  :: ',exact(3),' +/- ',exerr(3)
	write (70,'(T5,2(A,G10.4))') 'FOURTH CENTRAL MOMENT  :: ',exact(4),' +/- ',exerr(4)
	write (70,'(T5,2(A,G10.4))') 'FIFTH  CENTRAL MOMENT  :: ',exact(5),' +/- ',exerr(5)
	write (70,*) 
       END IF
c
c Do calculations one last time when exiting to write to output files
c
      IF (cal_rqt.EQ.1) THEN
	u22= u2*QQ**2
	u33= U3*QQ**2
	u44= U4*QQ**4
	u55= u5*QQ**2
	step=(hiw-low)/(nmesh-1)
	DO i=1,nmesh
	  Y = low+(i-1)*step
	  xtable(i)=Y
	  wd = (Y/const4)/SQRT(u22)
	  Sia(i)= Sq/SQRT(2*pi*u22)*EXP(-(Y/const4)**2/(2*u22))/const4
	  S1(i) = -u33*wd/(2*U22**1.5)*(1 - wd**2/3)*Sia(i)
	  S2(i) = u44/(8*u22**2)*(1 - 2*wd**2 + wd**4/3)*Sia(i)
	  S3(i) = u55*wd/(8*u22**2.5)*(1 - 2/3*wd**2 + wd**4/15)*Sia(i)
	  sqe(i) = sia(i)+s1(i)+s2(i)  !+s3(i)
	END DO
	IF (dec_sqw.EQ.1) write (20,'(3G15.4)') ((xtable(i),sqe(i),0.05*sqe(i)),i=1,nmesh)
	IF (dec_sia.EQ.1) write (40,'(3G15.4)') ((xtable(i),sia(i),0.05*sia(i)),i=1,nmesh)
	IF (dec_s1.EQ.1) write (60,'(3G15.4)') ((xtable(i),s1(i),0.05*s1(i)),i=1,nmesh)
	IF (dec_s2.EQ.1) write (90,'(3G15.4)') ((xtable(i),s2(i),0.05*s2(i)),i=1,nmesh)
	IF (dec_s3.EQ.1) write (95,'(3G15.4)') ((xtable(i),s3(i),0.05*s3(i)),i=1,nmesh)
      END IF
      RETURN                                                              
      END                                                               
