*____________________________________________________________________________
        SUBROUTINE JQY_CAL(ycal,nv,inv,nptot,p,x,npk,ypk)
*_____________________________________________________________________________
c
c Determines parameters for J(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(1000),ypk(1000,20),ycal(1000),p(50),bc,bs
	integer*4 inv(1000),lpt,nd,npk,nv,nptot,nrd
        real*8 rx(1000),ry(1000),ydat(1000)
        logical moms
c Frills common blocks
        common /rdats/ nrd,rx,ry
        common /dats/ ydat
        common /logf/ moms
c General Variables
        real*8 n0,kc,A,a2,a4,a6,b3,b4,b5,b6
        real*8 jqy,jqy1,jqy2,jqy3,qq,Sq,temp
        real*8 x1(1000),temp1(1000),temp2(1000),temp3(1000),temp4(1000),temp5(1000)
        integer*4 i,j,nord
        logical firstime,comp,resfile
c  General common blocks
        common /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6
        common /qval/ qq,Sq,temp
        common /ordr/ nord
        common /jqys/ jqy,jqy1,jqy2,jqy3
        common /logc/ firstime,comp,resfile
c
c Assign fitting parameters from Frills
c
        bc = p(1)
        bs = p(2)
        n0 = p(3)*0.01
        kc = p(4)
        a2 = p(5)
        b3 = p(6)
        a4 = p(7)
        b4 = p(8)
        b5 = p(9)
        a6 = p(10)
        b6 = p(11)
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
           temp5(i) = ydat(j)
           IF (comp) THEN               ! Use option to separate J(Q,Y) comps
              CALL JQY_COMP(x(j))
              temp1(i) = jqy1
              temp2(i) = jqy2
              temp3(i) = jqy3
              ycal(j)  = jqy            
              npk = 4
           ELSE                         ! Use option for total J(Q,Y)
              CALL JQY_TOT(x(j))
              ycal(j) = jqy
              ypk(j,2) = jqy
              temp1(i) = jqy
              npk = 2
           END IF
           temp4(i) = ycal(j)
        END DO
c
c Convolute Calculated function with Instrumental Resolution if present 
c
        IF (resfile) THEN
           IF (comp) THEN
              CALL CONVOLVE(nv,x1,temp1,nrd,rx,ry)
              CALL CONVOLVE(nv,x1,temp2,nrd,rx,ry)
              CALL CONVOLVE(nv,x1,temp3,nrd,rx,ry)
              DO i = 1,nv
		 j = inv(i)
                 ypk(j,2) = temp1(i)
                 ypk(j,3) = temp2(i)
                 ypk(j,4) = temp3(i)
                 ycal(j) = temp1(i)+temp2(i)+temp3(i)+ypk(j,1)
		 temp1(i) = temp1(i)+temp2(i)+temp3(i) 
              END DO
           ELSE
              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
        END IF
c
c If requested, calculate and display Experimental and calculated moments
c
        IF (moms) THEN                  
           CALL JQY_MOMS(nv,x1,temp5,temp4,temp1)                
        END IF

        RETURN
        END
c
c
*_____________________________________________________________________________
        SUBROUTINE JQY_COMP(Y)
*_____________________________________________________________________________
c
c Program to calculate J(Q,Y) at a given Q as a FT of S(Q,s).
c
c ** This version separates the 3 components of J(Q,Y) which are
c       1) Condensate
c       2) Condensate Induced Singularity
c       3) Uncondensed components
c
c Therefore these components can be plotted separately
c
        implicit none
        real*8 qq,Sq,Y,Yv
        real*8 n0,kc,A,a2,a4,a6,b3,b4,b5,b6,n00,kcc
        real*8 Temp,mass,Vs,rho,kB,pi,hbar,ec,amu,te
        real*8 bound,epa,epr,wi(2000),err,wbar
        real*8 I0,const1,const2,const3,const4,const5
        real*8 jqy,jqy1,jqy2,jqy3
        real*8 F1,F2,F3,FOFS,alpha,gama
        integer*4 lw,liw,iw(500),inf,ifail,i,j,mn,lpt,nord,last,neval
        external F1,F2,F3,FOFS,alpha,gama
        logical firstime,comp,resfile,moms,init_fofs
c
        common /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6
        common /fcon/ Temp,mass,Vs,rho,kB,pi,hbar,ec,amu
        common /ycons/ const1,const2,const3,const4,const5
        common /FTan/ lw,liw,ifail,inf,bound,epa,epr
        common /jqys/ jqy,jqy1,jqy2,jqy3
        common /ordr/ nord
        common /qval/ qq,Sq,te
        common /yval/ Yv
c
        common /logf/ moms
        common /logc/ firstime,comp,resfile
	common /ft_fs/ init_fofs
c
        Yv = Y
        temp = te
        IF (firstime) THEN              !setup some initial patameters
c
c  Define some fundamental constants
c
           mass = 4.0026
           const1 = 2.072914                       ! Ei = const1*ki**2
           const2 = const1*1.008665/mass           ! E  = const2*QQ**2
           const3 = const1*1.008665*2
           const4 = mass/(QQ*const3)               ! Y  = const4*(w-wr)
           pi = 4.0*ATAN(1.0)
           hbar = 1.0545919                        ! *10**(-34) Js
           ec   = 1.6021917                        ! *10**(-19) C
           amu  = 1.660531                         ! *10**(-27) Kg
           kB   = 1.380622                         ! *10**(-23) J/K
           rho  = 0.6022169/27.0                   ! 1/A**3         He4 at SVP
           Vs   = 240                              ! m/s            He4 at SVP
           const5 = (ec*mass*amu/hbar**2)**2/100.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   n(s) = n0[1+f(s)] + A*n*(s)
c      where A = 1 - n0[1+I0]   and   I0=f(s=0)
c
           nord = 1
	   init_fofs = .true.
           I0 = fofs(0.0)
           A = 1 - n0*(1+I0)
           kcc = kc
           firstime = .false.      
        END IF                                  ! end of initialisation
        A = 1 - n0*(1+I0)
c
c  Calculate new value for A whenever kc (I0s) changes
c
        IF (kcc.NE.kc) THEN
           nord = 1
	   init_fofs = .true.
           I0 = fofs(0.0)
           A = 1 - n0*(1+I0)
           kcc = kc
        END IF
c       write (*,*) ' jqy_comp : Fofs(0) = ',Fofs(0.0)
c       write (*,*) ' jqy_comp : Fofs(2) = ',Fofs(2.0)
c       write (*,*) ' jqy_comp : Fofs(5) = ',Fofs(3.0)
	nord = 1


        CALL DQAGI(F1,bound,inf,epa,epr,jqy1,err,neval,ifail,liw,lw,last,iw,iw)
c        CALL D01AMF(F1,bound,inf,epa,epr,jqy1,err,wi,lw,iw,liw,ifail)
        jqy1 = n0*SQ/pi*jqy1
        CALL DQAGI(F2,bound,inf,epa,epr,jqy2,err,neval,ifail,liw,lw,last,iw,iw)
c        CALL D01AMF(F2,bound,inf,epa,epr,jqy2,err,wi,lw,iw,liw,ifail)
        jqy2 = n0*SQ/pi*jqy2
        CALL DQAGI(F3,bound,inf,epa,epr,jqy3,err,neval,ifail,liw,lw,last,iw,iw)
c        CALL D01AMF(F3,bound,inf,epa,epr,jqy3,err,wi,lw,iw,liw,ifail)
        jqy3 = A*SQ/pi*jqy3
        jqy = jqy1+jqy2+jqy3
        RETURN
        END
c
c
*_____________________________________________________________________________
        SUBROUTINE JQY_TOT(Y)
*_____________________________________________________________________________
c
c Program to calculate J(Q,Y) at a given Q as a FT of S(Q,s).
c
c ** This version calculates the total J(Q,Y) ie does NOT separately calculates
c    the Condensate, Condensate Induced Singularity and Uncondensed components
c
c Therefore these components CANNOT be plotted separately
c
        implicit none
        real*8 Qq,Sq,Y,Yv,te
        real*8 n0,kc,A,a2,a4,a6,b3,b4,b5,b6,kcc
        real*8 Temp,mass,Vs,rho,kB,pi,hbar,ec,amu
        real*8 bound,epa,epr,wi(2000),err,wbar
        real*8 I0,const1,const2,const3,const4,const5
        real*8 jqy,jqy1,jqy2,jqy3
        real*8 F4,FOFS
c        real D01AMF
        integer*4 lw,liw,iw(500),inf,ifail,i,j,mn,lpt,nord,last,neval
        external F4,FOFS
        logical firstime,comp,resfile,moms,init_fofs
c
        common /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6
        common /fcon/ Temp,mass,Vs,rho,kB,pi,hbar,ec,amu
        common /ycons/ const1,const2,const3,const4,const5
        common /FTan/ lw,liw,ifail,inf,bound,epa,epr
        common /jqys/ jqy,jqy1,jqy2,jqy3
        common /ordr/ nord
        common /qval/ qq,Sq,te
        common /yval/ Yv
c
        common /logf/ moms
        common /logc/ firstime,comp,resfile
	common /ft_fs/ init_fofs
c
        yv = Y
        temp = te
        IF (firstime) THEN              !setup some initialize parameters
c
c  Define some fundamental constants
c
           mass = 4.0026
           const1 = 2.072914                       ! Ei = const1*ki**2
           const2 = const1*1.008665/mass           ! E  = const2*QQ**2
           const3 = const1*1.008665*2
           const4 = mass/(QQ*const3)               ! Y  = const4*(w-wr)
           pi = 4.0*ATAN(1.0)
           hbar = 1.0545919                        ! *10**(-34) Js
           ec   = 1.6021917                        ! *10**(-19) C
           amu  = 1.660531                         ! *10**(-27) Kg
           kB   = 1.380622                         ! *10**(-23) J/K
           rho  = 0.6022169/27.0                   ! 1/A**3         He4 at SVP
           Vs   = 240                              ! m/s            He4 at SVP
           const5 = (ec*mass*amu/hbar**2)**2/100.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   n(s) = n0[1+f(s)] + A*n*(s)
c      where A = 1 - n0[1+I0]   and   I0=f(s=0)
c
           nord = 1
	   init_fofs = .true.
           I0 = FOFS(0.0)
           A = 1 - n0*(1+I0)
           kcc = kc
           firstime = .false.      
        END IF                                  ! end of initialisation
        A = 1 - n0*(1+I0)
c
c  Calculate new value for A whenever kc (I0) changes
c
        IF (kcc.NE.kc) THEN
           nord = 1
	   init_fofs = .true.
           I0 = FOFS(0.0)
           A = 1 - n0*(1+I0)
           kcc = kc
        END IF
	nord = 1
        CALL DQAGI(F4,bound,inf,epa,epr,jqy,err,neval,ifail,liw,lw,last,iw,iw)
c        CALL D01AMF(F4,bound,inf,epa,epr,jqy,err,wi,lw,iw,liw,ifail)
c       write (*,*) ' JQY_TOT : jqy = ',jqy
        jqy = SQ/pi*jqy
        RETURN
        END
c
c
*_____________________________________________________________________________
        REAL FUNCTION F1*8(s)           ! alpha(s)*beta(s)
                                        ! Condensate Component
*_____________________________________________________________________________
	implicit none
        real*8 s
        real*8 ALPHA,BETA
        external ALPHA,BETA
        F1 = ALPHA(s)*BETA(s)
        RETURN
        END
c
c
*_____________________________________________________________________________
        REAL FUNCTION F2*8(s)           ! fofs(s)*alpha(s)*beta(s)
                                        ! condensate induced 'peaking' component
*_____________________________________________________________________________
	implicit none
        real*8 s
        real*8 ALPHA,BETA,FOFS
        external ALPHA,BETA,FOFS
        F2 = FOFS(s)*ALPHA(s)*BETA(s)
        RETURN
        END
c
c
*_____________________________________________________________________________
        REAL FUNCTION F3*8(s)           ! alpha(s)*gama(s)
                                        ! Uncondensed component
*_____________________________________________________________________________
	implicit none
        real*8 s
        real*8 ALPHA,GAMA
        external ALPHA,GAMA
        F3 = ALPHA(s)*GAMA(s)
        RETURN
        END
*_____________________________________________________________________________
        REAL FUNCTION F4*8(s)           ! alpha*[n0*beta + n0*beta*f + A*gama]
                                        ! Total J(Q,s)
*_____________________________________________________________________________
	implicit none
        real*8 s,n0,kc,A,a2,a4,a6,b3,b4,b5,b6
        real*8 ALPHA,BETA,GAMA,FOFS
        external ALPHA,BETA,GAMA,FOFS
        common /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6
        F4 = ALPHA(s)*( n0*BETA(s)*(1+ FOFS(s)) + A*GAMA(s) )
        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,nord,l1,l2
        real*8 n0,kc,A,a2,a4,a6,b3,b4,b5,b6
        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 I0,I2,I4,t1,t2,t3,wbar,mom,error
        real*8 qq,Sq,temp
        real*8 FOFS 
        external FOFS
        character*80 line0,line1,line2,file1*60,file2*60
        logical firstime,comp,resfile,moms
c
        common /parm/ n0,kc,A,a2,a4,a6,b3,b4,b5,b6
        common /ordr/ nord
        common /file/ file1,l1,file2,l2
        common /qval/ qq,Sq,temp
c
        common /logf/ moms
        common /logc/ firstime,comp,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
c       IF (resfile) THEN
        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       END IF
c
c Evaluate the expected (Exact) moments from the input parameters used to
c calculate the fitting function
c
        nord = 1
        I0 = FOFS(0.0)
        nord = 3
        I2 = FOFS(0.0)
        nord = 5
        I4 = FOFS(0.0)
c
        M_exact(1) = 1.0                                        ! Zeroth moment
        M_exact(2) = 0.0                                        ! First moment
        M_exact(3) = n0*I2/3.0 + A*a2                           ! Second moment
        M_exact(4) = b3                                         ! Third moment
        M_exact(5) = b4 + n0*I4/5 + A*(3*a2**2+a4)              ! Fourth moment
        M_exact(6) = b5 + 10.0*(n0*I2/3*b3 + A*a2*b3)           ! Fifth moment
        t1 = b6 + A*a6 + 10*b3**2       !** neglect term n0I6/7 **
        t2 = 15.0*b4*(n0*I2/3 + A*a2)
        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)
        write (*,'(A)') ' # ** Resn file = '//file2(1:l2)
        write (*,*) '# '
        write (*,'(3(A,F6.3))') ' # ** Q = ',QQ,'  ** S(Q) = ',Sq,'  ** Temp. = ',temp    
	write (*,'(3(A,F6.3))') ' # ** n0= ',n0,'  **  Io  = ',I0,'  **    A  = ',A 
        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 (*,'(A)') ' #** NOTE: A term (n0*I6/7) was neglected from the exact M6 '
        write (*,'(A)') ' #   _____ I6 is the sixth moment of f(s) '
        write (*,*) '#_____________________________________________________________________________'
        RETURN
        END














