C+
C NAME:
C	fit_fft2
C PURPOSE:
C	To make a Levenberg-Marquardt fit of the fft IPS power spectrum.
C
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	run fit_fft2
C INPUTS:
C	Main program
C OUTPUTS:
C	Reads in an IPS power spectrum and fits coefficient parameters and determines errors in the power spectrum fit.
C FUNCTIONS/SUBROUTINES:
C	mrqmin.f
C	mrqcof.f
C	gaussj.f
C	covsrt.f
C	funcs.f
C	make_model.f
C	make_derivatives.f	
C	read_fft.f
C	
C INCLUDE:
C COMMON BLOCKS:
C SIDE EFFECTS:
C RESTRICTIONS:
C PROCEDURE:
C	Latest 10/11/13 work
C	logical $CAM.
C MODIFICATION HISTORY:
C	October-2013, Bernard Jackson (UCSD)
C-
	program fft_fit2
C
	implicit real*8(a-g,o-z)
	parameter      (NdataM	=  1000,	! Max # data points
     &			MAM	=  4,		! Max # of coefficients
     &			MfitM	=  4,		! Max # of coefficients fit
     &			NspecM	=  4096)	! Max # of input spectral values
C
	real*8		X	(NdataM),	! Values of numbers input (spectral values)
     &			Xin	(NdataM),	! Input set of numbers (same as above but passed to mrqmin)
     &			Ydat	(NdataM),	! Data values at each value of Xin
     &			Ymodf	(NdataM),	! Model values at each value of Xin
     &			Y	(NdataM),	! The function values
     &			Sig	(NdataM),	! Individual standard deviations of the values fit
     &			A	(MAM),		! Coefficients to be fit
     &			Covar	(MfitM,MfitM),	! scratch array	
     &			Alpha	(MfitM,MfitM)	! scratch array
C
	real		hspecval(NspecM),	! The input spectral frequencies read in
     &			hpower	(NspecM)	! The input spectral power read in
	real		haspecval(NspecM),	! The input spectral frequencies read in
     &			hapower	(NspecM)	! The input spectral power read in
C
	integer		ListA	(MAM)		! Order of coefficients fit
	character 	src*8			! source name
	character 	srcl*8			! last source name
	character	cdate*18
	character	cWildspectrum*80,
     &			cVar(3)*100,
     &			cArg*100
	character	clist*80(1000)
	character	cfile*25		/'./data/fft_T20121016.data'/
	character	cfile1*26		/'./data/T20121016103900.dat'/
	character	cfilei*44
	character	cfilem*44
	character	cfileo*44
	character	cfilep*48
C
	common		elong,as,be, 		! extra parameters used by make_model
     &                  zs,ze,xs,xe,ccc,PI,au,au2,anu,alamb,freq,twopi,qxf,
     &			L1,L2,L3,L4,L5,L6,L7,L8,L9,A1,A2,A3,A4,A5,A6,A7,A8,A9,
     &			z,z0,R,p,Vx,Xnorm,XMnorm,ahpow,anoise,
     &			er1,er2,er3,er4,er5,er6,er7,er8,er9,er10,
     &			ijsing,iFs(500),iFt(500),iFd(500),aispecq(500),ispecq1,iq1(500),ispecq2,iq2(500)
C
	external	Funcs				! Put in to declare an external function by B. Jackson - needed here
C
C set up common block and model constant parameters
C
	c 	= 2.99793d8				! the speed of light in m/sec
	PI 	= 3.141592653589d0
	au 	= 149597870700.d0			! AU in m
	au2	= au*au					! AU^2
	hfreq	= 327.0					! array frequency in MHz
	as	= 1.0d-6				! lower limit on qy (in 1/m) ~ 1/(1000 km)
	be	= 500.0d-6				! upper limit on qy (in 1/m) ~ 1/(2 km)
	zs 	= -2.0d0				! beginning z
	ze	= 0.0d0					! ending z
        isrc 	= 1					! source counter
C
	Ndata 	= 100					! # of actual data points
	MA 	= MAM					! # of coefficients (4)
	Mfit	= MfitM					! # number of coeficients fit
	NCA	= Mfit					! a # greater or equal to Mfit
C
	src 	= '        '
	srcl	= '        '	
C
	Nspectry = NspecM
C
	iq = 0
	print *, ' '
	write(*,'(15X,A,I5)') 'Number of input spectral values requested ',Nspectry
C
C Guess the values of the coefficients.
C
C read in the fft power spectrum
C
	iVar = 3
C	call ForeignArg(' ',iVar,cVar,cArg)
C	call ForeignFile(iVar,cVar,'time_series',cWildspectrum)
	cWildspectrum = './data/T20121016103900.dat'
	cWildspectrum = './data/T??????????????.dat'
	call file_list(cWildspectrum,clist)
C
      do J=1, 1000
	write(cfile1(1:26),'(2A)') '.',clist(J)(49:73)
	call read_tsh(cfile1,src,helong,iyess)
C
	cfile = cfile1
	write(cfile(1:25),'(4A)') cfile(1:7),'fft_T',cfile(9:16),'.data'
	call read_fft(cfile,hspecval,hpower,Nspectry,Nspec,iyes)
	if(iyes.ne.0) go to 22
C
	if(src .ne. srcl) iq = 0
	srcl = src
C
	read(cfile1(9:12),'(I4)')  iyr
	read(cfile1(13:14),'(I2)') imo
	read(cfile1(15:16),'(I2)') ida
	read(cfile1(17:22),'(I4)') iut
	write(cdate(1:18),'(I4,A,I2,A,I2,A,I4,A)') iyr,'/',imo,'/',ida,' ',iut,' UT'
	write(cfilei(1:6),'(A)') './out/'
	write(cfilem(1:6),'(A)') './out/'
	write(cfileo(1:6),'(A)') './out/'
	write(cfilep(1:6),'(A)') './out/'
	write(cfilei(7:44),'(A,A8,A,I4.4,A,I2.2,A,I2.2,A,I4.4)')   'obs_spec___',src,'_',iyr,'_',imo,'_',ida,'_',iut
	write(cfilem(7:44),'(A,A8,A,I4.4,A,I2.2,A,I2.2,A,I4.4)')   'model_spec_',src,'_',iyr,'_',imo,'_',ida,'_',iut
	write(cfileo(7:44),'(A,A8,A,I4.4,A,I2.2,A,I2.2,A,I4.4)')   'output_____',src,'_',iyr,'_',imo,'_',ida,'_',iut
	write(cfilep(7:48),'(A,A8,A,I4.4,A,I2.2,A,I2.2,A,I4.4)')   'output_synopsis',src,' ',iyr,' ',imo,' ',ida,' ',iut
C
	Nsteps = 40
	print *, ' '
	if(Nspec .ne. 256) then
	  Nsteps = Nsteps*Nspec/256
	  write(*,'(A,I5,A,I4)') ' Unexpected number (',Nspec,')of spectral points. Nstep set to', Nsteps	  
	end if 
        if(iq.eq.0) call AskI4('Spectral values between 0.1 and 10.0 Hz to average?',Nsteps)
	call saverage(hspecval,hpower,Nspectry,Nsteps,haspecval,hapower,Ndat)
C	print *, Ndat
C
C  Determine the noise from the spectrum beyond 10.0 Hz assuming the noise level constant over frequency.
C  Also determine the level of sigma throughout the spectrum from the same assumption. 
C
	do i=1,Nspectry							! initalize spectral and power values to zero
	  hspecval(i) = 0.0
	  hpower(i)   = 0.0
	end do
C	iwrite = 0
	iwrite = 1
	if(iwrite.eq.1) write(*,'(/,A)') 'The averaged spectral values are:'
	hanoise = 0.0
	hnoise  = 0.0		
	istart  = 0
	hXnorm  = 0.0
	do i=1,Ndat
	  if(iwrite.eq.1) write(*,'(I4,2F15.7)') i,haspecval(i),hapower(i)
C
	  if(hXnorm .eq. 0.0 .and. haspecval(i) .gt. 0.29) hXnorm = haspecval(i) ! Presumed value of spectrum to normalize to ~.3
C
	  if(haspecval(i).ge. 10.0 .and. haspecval(i) .le. 18.0) then		! Begin to determine the noise level after 10.0 Hz
	    hanoise = hanoise + hapower(i)
	    hnoise  = hnoise + 1.0
	    if(istart.eq.0) istart = i
	  end if
	  hspecval(i) = haspecval(i)
	  hpower(i)   = hapower(i)
	end do
	hanoise = hanoise/hnoise
	anoise = dble(hanoise)							! Average value of the noise
	hsigv1 = 0.0
	print *, ' '
	print *, 'Variation of noise at the end of the spectrum'
	do i=istart,Ndat
	  hsigv1 = (hapower(i) - hanoise)**2 + hsigv1
	  print *, '              ',hapower(i) - hanoise
	end do
C	anoise = 0.0d0						! Put in to try
C	sigv = dble(sqrt(hsigv1/hnoise))*10.0d0			! now different
	sigv = anoise
	print *, ' '
	write(*,'(A,F9.6,A,F9.6))') ' Normalized noise =',anoise,' sigma value =',sigv
C
C  Now begin the parameter values to fit and how.
C
	Ndata = Nspec							! Number of spectral data values to fit
	print *, ' '
        if(iq.eq.0) call AskI4('Number of spectral values possible to fit?',Ndata)
	if(Ndata.gt.Nspec) Ndata = Nspec				! The number of values to fit must be less that the values present
C
C	hxs = hspecval(1)
	hxs = 0.29							! Lower limit on what is fit by Marquardt
        if(iq.eq.0) call AskR4('Beginning spectral value greater than?',hxs)
	xs = dble(hxs)
C
C	hxe = hspecval(Ndata)
C	hxe = 5.0							! Upper limit on what is fit by Marquardt
	hxe = 4.0							! This depends on the radio frequency observation, some
        if(iq.eq.0) call AskR4('Ending spectral value less than?',hxe)
	xe = dble(hxe)
C	ispec = 0
	ispec = 1
        if(iq.eq.0) call AskI4('Do you want to print out the spectral values, 1 = Yes?',ispec)
        if(iq.eq.0) call AskI4('Total number of coefficients?',MA)
C
	print *, ' '
        if(iq.eq.0) call AskR4('Radio source frequency for this radio system?',hfreq)
	freq = dble(hfreq)
C
	do i=1, Ndata
	  Xin(i)  = dble(hspecval(i))
	  Ydat(i) = dble(hpower(i))
C	  Sig(i)  = 0.0005d0
C	  Sig(i) = sigv					! weight evenly by noise throughout spectrum.
	  Sig(i) = Ydat(i) * 0.2d0			! weight by 20% of the observed power
C	  if(i+2.le.Ndata) then				! weight sig by 1/power observed normalized to third point
C	    Sig(i)  = sigv/Ydat(i+2) 
C	  else
C	    Sig(i)  = sigv/Ydat(Ndata)
C	  end if
	end do
C
C	initial model coefficient parameters
C
	hvel	= 500.0
	hdvel	= 0.001	
	halphav	= 3.5
	hdalphav= 0.001
	htheta	= 0.1
	hdtheta	= 0.001
	har	= 1.5
	hdar	= 0.001
C
	print *, ' '
	write(*,'(2A)'),      '       Source name   ', src
	write(*,'(A,F10.4)'), '       Elongation =', helong
	elong = dble(helong)
C
	Mfit = MA
	if(src .eq. ' 3C273  ')	then
	  Mfit 	 = 3
	  htheta = 0.001
	end if
	print *, ' '
        if(iq.eq.0) call AskI4('Number of coefficients to fit?',Mfit)
C
	print *, ' '
        if(iq.eq.0) call AskR4('Input value of coefficient 1 V?',hvel)
        if(iq.eq.0) call AskR4('Fractional derivative change of V?',hdvel)
	vel   = dble(hvel)*1.0d3						! velocity in m/sec
	er1   = dble(hdvel)
	print *, ' '
        if(iq.eq.0) call AskR4('Input value of coefficient 2 alpha?',halphav)
        if(iq.eq.0) call AskR4('Fractional derivative change of alpha?',hdalphav)
	alphav = dble(halphav)
	er2    = dble(hdalphav)
	print *, ' '
        if(iq.eq.0) call AskR4('Input value of coefficient 3 theta?',htheta)
        if(iq.eq.0) call AskR4('Fractional derivative change of theta?',hdtheta)
	theta = dble(htheta) * 4.84813681d-6					! theta input in radians (converted from arc seconds)	
	er3   = dble(hdtheta)
	print *, ' '
        if(iq.eq.0) call AskR4('Input value of coefficient 4 axial-ratio?',har)
        if(iq.eq.0) call AskR4('Fractional derivative change of axial-ratio?',hdar)
	ar = dble(har)
	er4 = dble(hdar)
C
	A(1) = vel
	A(2) = alphav
	A(3) = theta
	A(4) = ar
C
 35	continue
	print *, ' '
        if(iq.eq.0) call AskR4('Frequency at which to normalize the spectrum?',hXnorm)
	Xnorm = dble(hXnorm)
	hpow  = 0.0
	htiny = 1.0E-5
	hXnor = hXnorm
	do i=1,Ndat
	  if(hpow .eq. 0.0 .and. haspecval(i) .gt. (hXnorm - htiny)) then 		! Set the value of the power to normalize to
	    if(hapower(i) .lt. (1.0 - htiny) .or. hapower(i) .gt. (1.0 + tiny)) then
	      print *, ' The normalized power at this frequency is not 1.0'
	      call AskR4('What is the frequency you want to normalize power to?',hXnor)
	    end if
	    if(hXnor .eq. hXnorm) then	
	      hpow = hapower(i)	
	    else
	      go to 35
	    end if
	  end if	     
	end do
C
C Order to fit the parameters
C outside qsimp
 20	continue
	Lista(1) = 1
	Lista(2) = 2
	Lista(3) = 4
	Lista(4) = 3
	print *, ' '
        if(iq.eq.0) call AskI4('Order to fit coefficient 1 V?',Lista(1))
        if(iq.eq.0) call AskI4('Order to fit coefficient 2 alpha?',Lista(2))
        if(iq.eq.0) call AskI4('Order to fit coefficient 3 theta?',Lista(3))
        if(iq.eq.0) call AskI4('Order to fit coefficient 4 axial-ratio?',Lista(4))
C
	if((Lista(1)+Lista(2)+Lista(3)+Lista(4)).ne.10) then
	  print *, 'Parameter order values not correct, try again'
	  go to 20
	end if
C
 29	continue
C
	L1 = Lista(1)
	L2 = Lista(2)	
	L3 = Lista(3)
	L4 = Lista(4)
C
C Set up other additional model parameters
C
	anu 	= freq * 1.0d6				! array frequency in Hz
	alamb	= c / anu				! wavelength of radio waves in m
	print *, ' '
	print *, 'Wavelength of radio waves at this frequency in m', alamb
	twopi = 2.0d0 * PI
	elong	= elong * (PI/180.0d0)			! source elongation in radians					
C
	Alamda	 = -1.0d0				! first time through to set up Marquardt
	Alamdal	 = 0.0009d0
	Chisq 	 = 1.0d16
	Chisqold = Chisq
C
C Determine the initial normalization value
C
C Call the Minimization
C
	ijsing = 0			
	do jj=1,50					! initalize error messages for each iteration
	  iFs(jj)	= 0
	  iFt(jj)	= 0
	  iFd(jj)	= 0
	  iq1(jj)	= 0
	  iq2(jj)	= 0
	end do
C
	ahpow  = dble(hpow)
	XMnorm = -ahpow				! If - Normalize the model at the model value of hpow at frequency Xnorm
	call make_model(Xnorm,A,XMMnorm,MA)                
	XMnorm = XMMnorm
	print *, ' '
	write(*,'(A,E15.7,A,F11.7,A,F11.7)') 'XMMnorm =', XMnorm, ' At frequency',Xnorm,', to get',hpow
C
	inc = 0
C	go to 30
	Nint = 30
	do i=1,Nint	
C
C 	  print *, 'Outside mnqmin',A
C
	  print *, ' '
	  Write(*,'(A)') '*************************************************************************'
	  if(Chisq.eq. 1.0d16 .and. Chisqold.eq. 1.0d16)
     &	      write(*,'(A,I3,A)') 'Begin Loop', i,' Initalize values for Marquardt and begin first try.'
	  if(Chisq.ne. 1.0d16 .and. Chisqold.eq. 1.0d16)
     &	      write(*,'(A,I3,A,F15.3)') 'Begin Loop', i,' Chisq  =',Chisq
	  if(Chisq.ne. 1.0d16 .and. Chisqold.ne. 1.0d16)
     &	      write(*,'(A,I3,A,2F15.3)') 'Begin Loop', i,' Chisqold, Chisq  =',Chisqold, Chisq	    
	  Chisqold = Chisq
	  call mrqmin(Xin,Ydat,Sig,Ndata,A,MA,ListA,Mfit,Covar,Alpha,NCA,Chisq,Funcs,Alamda)	

	  if(Alamda .le. Alamdal) then
	    write(*,'(A,I3,A)') ' Marquardt #',i,' continues, alambda decreased'
	    iflag = 1
	  end if
	  if(Alamdal .lt. Alamda) then
     	    write(*,'(A,I3,A)') ' Marquardt #',i,' failed, alambda increased' 
	    iflag = 0
	  end if
	  if(iflag.eq.1) Chisqo = Chisqold
	  print *, 'Chisqold, Chisq  =',Chisqold, Chisq
	  Print *, 'After Mrqmin, new A =',A(1)/1.0d3
	  print *, 'Alamdal, Alamda  =',Alamdal,Alamda
C
	  ival	= i
	  if(((Chisqold-Chisq)/Chisq) .eq. 0.0d0 .and.ijsing.eq.1) then
		if(A(1).lt. 0.0d0) then
		  print *, 'Velocity < 0! use initial' 
		  A(1) = dble(hvel)*1.0d3
		  Lista(1) = 4
		  Lista(2) = 1
		  Lista(3) = 3
		  Lista(4) = 2
		end if
		if(A(2).lt. 0.0d0) then
		  print *, 'Alpha < 0! use initial' 
		  A(2) = dble(halphav)
		  Lista(2) = 4
		  Lista(1) = 1
		  Lista(3) = 3
		  Lista(4) = 2
		end if
		if(A(3).lt. 0.0d0) then
		  print *, 'Theta < 0! use initial' 
		  A(3) = dble(htheta)*4.84813681d-6
		  Lista(3) = 4
		  Lista(1) = 1
		  Lista(2) = 2
		  Lista(4) = 3
		end if
		if(A(4).lt. 0.0d0) then
		  print *, 'AR < 0! use initial' 
		  A(4) = dble(har)
		  Lista(4) = 4
		  Lista(1) = 1
		  Lista(2) = 2
		  Lista(3) = 3
		end if
		Mfit = Mfit - 1
                write(*,'(A,F10.4,F10.5,2F10.6)') ' singular matrix, new values of A',A(1)/1.0d3,A(2),A(3)/4.84813681d-6,A(4)
                write(*,'(A,I14)')                ' singular matrix, new Mfit       ', Mfit
	        go to 29
	  end if
	  if(((Chisqold-Chisq)/Chisq) .lt. 0.001d0 .and. Alamda .le. Alamdal ) go to 30
	  Alamdal = Alamda

	  if(i.eq.Nint) then
	    inc = 1
	    print *, '  '
	    print *, 'Program ended without converging completely'
	    Chisqold = Chisqo
	    print *, 'Chisqo',Chisqo
	    go to 30
	  end if
C 
	end do
 30	continue
	Alamdal = Alamda
	Chisql  = Chisq
	print *, '  '
	write(*,'(A,F15.7,A,F15.7)') ' Final value of Chisquare =',Chisq,', PERCENT from last =',100.d0*(Chisqold-Chisq)/Chisq
	print *, 'Outside last mrqmin'
	Alamda = 0.0d0
	call mrqmin(Xin,Ydat,Sig,Ndata,A,MA,ListA,Mfit,Covar,Alpha,NCA,Chisq,Funcs,Alamda) 
	print *, ' '
	if(inc.eq.0) print *, ' Iterations successfully completed'
	if(inc.eq.1) print *, ' Iterations completed - program ended without converging'
	print *, ' '
	write(*,'(4A)')' Source',src,' on ', cdate
	write(*,'(A)')' Final values Fit'
	write(*,'(A,F10.4,A,F10.4)') ' Velocity    =', A(1)/1.0d3,' +-', sqrt(Covar(1,1))/1.0d3
	write(*,'(A,F10.5,A,F10.5)') ' Alpha       =', A(2),' +-', sqrt(Covar(2,2))
	write(*,'(A,F10.6,A,F10.6)') ' Theta       =', A(3)/4.84813681d-6,' +-', sqrt(Covar(3,3))/4.84813681d-6
	write(*,'(A,F10.5,A,F10.5)') ' Axial Ratio =', A(4),' +-', sqrt(Covar(4,4))
	write(*,'(A,4E15.7)') 'Covarance values', Covar(1,1),Covar(2,2),Covar(3,3),Covar(4,4)
	write(*,'(A,4E15.7)') 'Alpha values', Alpha(1,1),Alpha(2,2),Alpha(3,3),Alpha(4,4)
C
	if(ispec.eq.1) then
	  call write_spec(cfilei,dble(hspecval),dble(hpower),Ndat)	! Write out to a file the input spectral data modeled
	  XMnorm = -1.0d0						! Normalize the model at frequency Xnorm
	  call make_model(Xnorm,A,XMMnorm,MA)                
	  XMnorm = XMMnorm
	  xe = dble(hspecval(Ndat))
	  print *, ' '
	  do i=1, Ndat
	    if(i.gt. 1 .and. (iq1(i-1).gt. 1000 .or. iq2(i-1) .gt. 1000 )) then
	      XMMnorm = 0.0d0
            else
	      call make_model(dble(hspecval(i)),A,XMMnorm,MA)
            end if
	    Ymodf(i) = XMMnorm
	    write(*,'(A,I3,F15.7,F15.7)') 'XMMnorm =', i, dble(hspecval(i)),XMMnorm
	  end do
	  call write_spec(cfilem,dble(hspecval),Ymodf,Ndat)		! Write out to a file the final model spectrum
	end if
	call write_val(cfileo,mfit,hvel,halphav,htheta,har,A,Covar,MA,Ndat,inc,iVal,Chisqold,Chisql,anoise,sigv,alamdal)
	call write_vals(cFilep,isrc,iUV,mfit,hvel,halphav,htheta,har,A,Covar,MA,Ndat,inc,iVal,Chisqold,Chisql,anoise,sigv,alamdal)
	isrc = isrc + 1
C
 22	continue
C
	iq = iq + 1
      end do
      close(iUV)
C
	print *, ' '
	print *, ' Program successfully completed'
	print *, ' '
C
      stop
      end
C************************************************************************************************
	subroutine Funcs(X,A,Xmod,Dxda,NA)
C
	implicit real*8(a-h,o-z)
	real*8	A	(NA),
     &		Dxda	(NA)
	real*8  er	(NA)

	common		elong,as,be, 		! extra parameters used by make_model
     &                  zs,ze,xs,xe,ccc,PI,au,au2,anu,alamb,freq,twopi,qxf,
     &			L1,L2,L3,L4,L5,L6,L7,L8,L9,A1,A2,A3,A4,A5,A6,A7,A8,A9,
     &			z,z0,R,p,Vx,Xnorm,XMnorm,ahpow,anoise,
     &			er1,er2,er3,er4,er5,er6,er7,er8,er9,er10,
     &			ijsing,iFs(500),iFt(500),iFd(500),aispecq(500),ispecq1,iq1(500),ispecq2,iq2(500)
C
C	print *, 'before make_model',A,er1,er2,er3,er4
C
	  call make_model(X,A,Xmod,NA)
C	  print *, ' '
C	  print *, 'after make_model',x,xmod
C	  print *, ' '
C
C	  call formula(X,A,Xmod,NA)
C	  Xmod = A(1)*X**2 + A(2)
C
C Determine A derivatives	subroutine make_model(X,A,Xmod,NA)
C
C	  er1, er2, etc., set to the order of A(i) in LISTA(i)
C
	  er(1) = er1				
	  er(2) = er2
	  er(3) = er3
	  er(4) = er4
C
C	  print *, 'Fractional E', er1,er2,er3,er4
C	  print *, 'A before der', A
C	  print *, 'New E       ', er
C
	  call derivatives(X,A,Xmod,Dxda,er,NA)
C
C	  write(*,'(A,E12.5,F9.3,F7.4,F12.6,F9.5,/,4E15.6)') 
C     &		'Inside Funcs at end', Xmod, A(1)/1.0d3, A(2), A(3)/4.84813681d-6, A(4), 
C     &		Dxda(1)*1.0d3, Dxda(2), Dxda(3)*4.84813681d-6, Dxda(4)
C
C end A derivative determination
C
	return
	end
C
C**************************************************************************************************
C
	subroutine make_model(X,A,Xmod,NA)
C
	implicit real*8(a-h,o-z)
	real*8		A	(NA)
	common		elong,as,be, 	! extra parameters used by the_model - zs and ze needed here, NA not doubly defined
     &                  zs,ze,xs,xe,ccc,PI,au,au2,anu,alamb,freq,twopi,qxf,
     &			L1,L2,L3,L4,L5,L6,L7,L8,L9,A1,A2,A3,A4,A5,A6,A7,A8,A9,
     &			z,z0,R,p,Vx,Xnorm,XMnorm,ahpow,anoise,
     &			er1,er2,er3,er4,er5,er6,er7,er8,er9,er10,
     &			ijsing,iFs(500),iFt(500),iFd(500),aispecq(500),ispecq1,iq1(500),ispecq2,iq2(500)
C
	external	ainteg1
	external	ainteg2
	NAA = NA
	A1 = A(1)			! Velocity
	if(A1 .gt. 2.0d6 ) then 
	  A1 = 2.0d6
	  ivul = 1
	end if
	if(A1 .lt. 2.0d5 ) then
	  A1 = 2.0d5
	  ivll = 1
	end if

	A2 = A(2)			! alpha
	if(A2 .gt. 5.0d0 ) then
	  A2 = 5.0d0
	  iaul = 1
	end if
	if(A2 .lt. 1.0d0 ) then
	  A2 = 1.0d0
	  iall = 1
	end if

	A3 = A(3)			! theta
	if(A3 .gt. 0.4d0*4.84813681d-6 ) then
	  A3 = 0.4d0*4.84813681d-6
	  itul = 1
	end if
	if(A3 .lt. 0.01d0*4.84813681d-6 ) then
	  A3 = 0.01d0*4.84813681d-6
	  itll = 1
	end if

	A4 = A(4)			! axial ratio
	if(A4 .gt. 5.0d0 ) then
	  A4 = 5.0d0
	  ixul = 1
	end if
	if(A4 .lt. 0.2d0 ) then
	  A4 = 0.2d0
	  ixll = 1
	end if
C
C The formula that fits the data
C
C	  Xmod = A(1)*X**2 + A(2)
C
	  qxf = (twopi)*X				! X in cycles/sec --> qxf

C	print *, 'outside qsimp for first integral', X,qxf,A,NA,zs,ze,xs,xe
C
	if(x.ge.xs.and.x.le.xe) then			! limit the beginning and end frequencies
C
	iold = 1
C	iold = 0
	if(iold.eq.1) go to 20
	  qs = as
	  qe = be
	  nzz = 100
	  zsau = zs*au					! z changed to meters here zs = -2 AU
	  zeau = ze*au
C	  dzz = (zeau-zsau)/float(nzz)			! Goes from -2 to 0
	  dzz = (zsau-zeau)/float(nzz)
C	  zz  = zsau					! Goes from -2 to 0
	  zz  = 0.0d0
	  xmod = 0.0d0
	  do i=1, nzz
	    if(i.eq.1) then
	      z  = zz
	      za = zz
	      call qsimp2(ainteg2,qs,qe,xmod1)
	    end if
	    zb  = zz + dzz
	    z   = zb
	    call qsimp2(ainteg2halphav,qs,qe,xmod2)
C	    xmod = 0.5*(zb-za)*(xmod1+xmod2) + xmod			! Goes from -2 to 0
	    xmod = 0.5*(za-zb)*(xmod1+xmod2) + xmod
C	    write(*,'(A,I5,2E10.3,E15.7/4E15.7)'), 'in make_model - i, qs,qe,xmod =',i,qs,qe,xmod,za,zb,xmod1,xmod2 
	    xmod1 = xmod2
	    za = zb
	    zz = zz + dzz
C	    print *, ' '
	  end do
 20	continue
	if(iold.eq.0) go to 21
	  qs = as
	  qe = be
	  nzz = 1
	  zsau = zs*au					! z changed to meters here zs = -2 AU
	  zeau = ze*au
	  dzz = (zsau-zeau)/float(nzz)			! Goes from 0 to -2
	  zz  = 0.0d0
	  xmod = 0.0d0
	  do i=1, nzz
	    if(i.eq.1) then
	      z  = zz
	      za = zz
	    end if
	    zb  = zz + dzz
	    z   = zb
	    call qsimp1(ainteg1,zb,za,xmod1)
	    xmod = xmod + xmod1
C	    write(*,'(A,I5,2E10.3,E15.7/4E15.7)'), 'in make_model - i, qs,qe,xmod =',i,qs,qe,xmod,za,zb,xmod1 
	    za = zb
	    zz = zz + dzz
C	    print *, ' '
	  end do
C	  call qsimp1(ainteg1,zsau,zeau,xmod)
 21	continue
C	  print *, ' '
C	  write(*,'(A,F12.7,E15.7)'), 'in make_model unnorm - X, xmod =', X, xmod
	  if(XMnorm .lt. 0.0d0) then		! normalize the model value to the normalized obs value of the input frequency X at XMMnorm
C           xmod = xmod + xmod*anoise
	    xmod = (xmod + xmod*anoise)/ahpow
C	    xmod = xmod/ahpow + xmod*anoise
	  else				
	    xmod = xmod/XMnorm + anoise
C	    xmod = (xmod + anoise)/XMnorm
	  end if
	end if

C	stop

	return
	end
C
C******************************************************
C
	function ainteg1(zz)
	implicit real*8(a-h,o-z)
C
	common		elong,as,be, 	! parameters used by the_model to calculate qyI
     &                  zs,ze,xs,xe,ccc,PI,au,au2,anu,alamb,freq,twopi,qxf,
     &			L1,L2,L3,L4,L5,L6,L7,L8,L9,A1,A2,A3,A4,A5,A6,A7,A8,A9,
     &			z,z0,R,p,Vx,Xnorm,XMnorm,ahpow,anoise,
     &			er1,er2,er3,er4,er5,er6,er7,er8,er9,er10,
     &			ijsing,iFs(500),iFt(500),iFd(500),aispecq(500),ispecq1,iq1(500),ispecq2,iq2(500)
C
	external ainteg2
C
	  z = zz
	  qs = as
	  qe = be
	  call qsimp2(ainteg2,qs,qe,xmod1)
	  ainteg1 = xmod1
	return
	end
C
C******************************************************
C
	function ainteg2(qy)
	implicit real*8(a-h,o-z)
C
	common		elong,as,be, 	! parameters used by the_model to calculate qyI
     &                  zs,ze,xs,xe,ccc,PI,au,au2,anu,alamb,freq,twopi,qxf,
     &			L1,L2,L3,L4,L5,L6,L7,L8,L9,A1,A2,A3,A4,A5,A6,A7,A8,A9,
     &			z,z0,R,p,Vx,Xnorm,XMnorm,ahpow,anoise,
     &			er1,er2,er3,er4,er5,er6,er7,er8,er9,er10,
     &			ijsing,iFs(500),iFt(500),iFd(500),aispecq(500),ispecq1,iq1(500),ispecq2,iq2(500)
C
		V 	= A1
		alpha	= A2
		theta	= A3
		ar	= A4
C
C	write(*,'(A,F9.3,F7.4,F12.6,F9.5)')' inside ainteg2',V/1.0d3,alpha,theta/4.84813681d-6,AR 
C	write(*,'(A,7E15.7)')' inside ainteg2',au,elong,qxf,alamb,PI,qy,z
C
C z come into the program in AU, z is 0 at the point p, (the z integral goes from z = -2 AU to 0 AU.)
C
		! z in m
C
		z0    = z + au*COS(elong)		! z0 in cm			

		p = au * TAN(elong)

		R2 = (z*z + p*p)
		R  = SQRT(R2)

		Vx = V * p / R				! V in cm/sec - Vx in cm/sec
C
		qx = qxf / Vx				! qxf in 1/sec --> qx in k if V in cm/sec (qy is integrated from 15 to 1500 km)
C
		qxSquared = qx*qx			! qx^2 in (1/km)^2
		qySquared = qy*qy			! qy^2 in (1/km)^2
		qSquared =  (qxSquared + qySquared)	! q^2 - (1/km)^2
C
C	print *, 'ainteg2 before Fd',z,R,V,Vx,qxf,qx,qy,qSquared,z0,alamb,PI
		Fd = SIN(qSquared * z0 * alamb / (4.0*PI))**2
C	if(Fd.eq.0.0) print *, 'Fd = 0.0, ainteg2 before Fd',z,R,V,Vx,qxf,qx,qy,qSquared,z0,alamb,PI
		if(Fd.eq.0.0d0)       iFd(ispecq1) = iFd(ispecq1) + 1
		aucost	= au * COS(elong)*theta/2.35
		aucost2	= aucost * aucost
		expval	= -qSquared * aucost2
C		expval	= -qSquared * (au * COS(elong)*theta/2.35)**2
		if(expval.lt.-600.0d0) then
C		  write(*,'(A,F8.1,A,E15.7)') 'Expval = ',expval,' limited to prevent underflow with small theta =', theta
		  expval = -600.0d0
		end if
C	if(Fd.eq.0.0.or.expval.eq.-40.0) print *, 'ainteg2 after Fd',Fd,qSquared,au,elong,COS(elong),theta,
C     &			-qSquared * (au * COS(elong)*theta/2.35)**2,expval
C	print *, 'ainteg2 after Fd',Fd,qSquared,au,elong,COS(elong),theta,
C     &			-qSquared * (au * COS(elong)*theta/2.35)**2,expval

		Fs = EXP(expval)

	if(Fs.eq.0.0) then 
C	  print *, 'ainteg2 after Fs',Fs,expval,R2,qxSquared,qySquared,AR,alpha,R2**(-2),
C     &			R2**(-2) * (qxSquared + qySquared / AR**2)**(-alpha/2.0d0)
	  iFs(ispecq1) = iFs(ispecq1) + 1
	  return
	end if

C		Ft = R2**(-2) * sqrt(qxSquared + qySquared/(AR**2))**(-alpha)
C		Ft = R2**(-2) * (qxSquared + qySquared / AR**2)**(-alpha/2.0d0)
		Ft = ((qxSquared + qySquared / (AR*AR))**(-alpha/2.0d0)) / (R2 * R2)
		if(Ft.eq.0.0d0)       iFt(ispecq1) = iFt(ispecq1) + 1
		ainteg2 = (Fd * Fs * Ft)/Vx
C
C	  print *, ' At end of ainteg2', V,alpha,theta,AR, z, qx,qy,Fd,Fs,Ft,ainteg2,aintegx
	if(Fd.eq.0.0.or.Fs.eq.0.0.or.Ft.eq.0.0.or.ainteg2.eq.0.0) then
C	  print *, ' At end of ainteg2', V,alpha,theta,AR, z, qx,qy,Fd,Fs,Ft,ainteg2,aintegx
C	  stop
	end if
C
	return
	end	
C
C******************************************************
C
	subroutine derivatives(X,A,Xmod,Dxda,er,NA)
	implicit real*8(a-h,o-z)
C
	real*8	A	(NA),
     &		Dxda	(NA)
	real*8  er	(NA),
     &		AA	(20)
C
C The derivatives
C
C	print *, 'inside derivatives', NA,A,er 
C
	do i = 1, NA
	  do j = 1, NA
	    AA(j) = A(j)
          end do
	  AA(i) = A(i) + A(i)*er(i)
	  call make_model(X,AA,XmodAP,NA)
C	  AA(i) = A(i) - A(i)*er(i)
C	  call make_model(X,AA,XmodAM,NA)
C	  Dxda(i) = (XmodAP-XmodAM)/(2.0d0*er(i)*A(i))
	  Dxda(i) = (XmodAP-Xmod)/(1.0d0*er(i)*A(i))
C	  if(i.eq.1) print *, 'X, DxDa Inside deriv',XmodAP, XmodAM, XmodAP-XmodAM, Dxda(i)
	end do
C
	return
	end