subroutine guess(acf,tau,npts,zero,amin,te,tr)
c
c     find zero crossing (zero), depth of minimum (amin), height of maximum
c
      real acf(npts),tau(npts)

      zero=0.0
      amin=1.0
      tmin=0.0
      jmin=0

      do i=npts,2,-1
         if(acf(i)*acf(i-1).lt.0.0) then
            zero=(tau(i-1)*acf(i)-tau(i)*acf(i-1))/(acf(i)-acf(i-1))
         end if
         if(acf(i).lt.amin) then
            amin=acf(i)
            jmin=i
         end if
      end do

      if(jmin.gt.0) then
         call parab1(tau(jmin-1),acf(jmin-1),a,b,c)
         tmin=-b/(2.0*a)
         amin=c+tmin*(b+tmin*a)
      end if
      
      tr=cdtr1(-amin)
      te=czte1(zero*1000.0,tr)
      return
      end
      
      subroutine parab1(x,y,a,b,c)
C-----
      dimension x(3),y(3)
      delta=x(1)-x(2)
      a=(y(1)-2.*y(2)+y(3))/(2.*delta*delta)
      b=(y(1)-y(2))/delta - a*(x(1)+x(2))
      c=y(1)-a*x(1)*x(1)-b*x(1)
      return
      end

      real function cdtr1(depth)
C-----convert depth to te/ti ratio
      dimension tr(4)
C     according to the curve published in farley et al 1967
c     modified for 2004 conditions on axis
c      data tr/7.31081,3.53286,5.92271,.174/
      data tr/9.5,4.0,8.5,.3/
      data nt/4/
      cdtr1=tr(1)
      do i=2,nt
        cdtr1=cdtr1*depth + tr(i)
      end do
      return
      end

      real function czte1(zlag,tr)
C-----convert zero crossing point to te
C     according to the curve published in farley et al 1967
c     modified for 2004 conditions on axis
      dimension dt(4)
c      data dt/0.00945025,-0.0774338,.203626,.812397/,nd/4/
      data dt/0.00945025,-0.0774338,.2,0.9/,nd/4/
      data t0/1000./
      tr1=min(abs(tr),5.)
      if(zlag .eq. 0)then
        czte1=1000000.
      else
        dt0=dt(1)
        do  i=2,nd
          dt0=dt0*tr1 + dt(i)
        end do
        czte1=t0*(dt0/zlag)**2
      end if
      return
      end

      subroutine fit(wl,taup,rhop,covar,cinv,sigma2p,paramp,ebp,
     &     bfldp,alphap,densp,alt,time,nl,ifitp,ist)
c
c     subroutine to fit measured ACF 
c     wavelength wl (m),lags taup (s), normalized acf rhop,
c     experimental variances sigma2p, covariances covar,
c     fit parameters params, error bars ebp, magnetic field bfldp (gauss),
c     alphap B-field angle (radians), density densp (gcs),
c     altitude alt (km), time (LT hours), nl lags
c     ifitp determines which parameters are fit (see below)
c
      real tol,pi,wl
      parameter(nlmax=100,npmax=10,lwa=2000,tol=1.0e-5)
      external fcn

      real covar(nl,nl),cinv(nl,nl)
      real ev(nlmax*nlmax),ap(nlmax*(nlmax+1)/2)
      real fv1(nlmax),fv2(nlmax),det(2),w(nlmax)

      real taup(nl),rhop(nl),sigma2p(nl),paramp(npmax)
      real p2(npmax),ebp(npmax),alphap,bfldp,densp
      real tau(nlmax),rho(nlmax),sigma2(nlmax),params(npmax)
      real fvec(nlmax),wa(lwa),wa2(lwa),eb(npmax)

      integer iwa(npmax),ifitp(npmax),ifit(npmax)

      include 'fitter.h'
      common /mode/imode

      common/fitter/tau,rho,sigma2,params,ifit
      common/trans/ev      
C-----

c     set ifit(1-5) to unity to fit:
c     1 normalization
c     2 Te 
c     3 Ti
c     4 H+
c     5 He+

      imode=2
c     
      pi=4.0*atan(1.0)
      ak=2.0*pi/wl

      wi(1)=16
      wi(2)=1
      wi(3)=4
      nion=3
      
c
c     invert covariances and find eigenvalues and eigenvectors
c
      l = 0
      do 20 j = 1, nl
         do 10 i = 1, j
            l = l + 1
            ap(l) = covar(i,j)
 10      continue
 20   continue      
      
      call sppfa(ap,nl,info)
      call sppdi(ap,nl,det,01)

      l = 0
      do 40 j = 1, nl
         do 30 i = 1, j
            l = l + 1
            cinv(i,j)=ap(l)
            cinv(j,i)=ap(l)
 30      continue
 40   continue   
      
c
c     transformation matrix is inverse (transpose) of eigenvectors
c
      
      matz=1
      call rs(nl,nl,cinv,w,matz,ev,fv1,fv2,ierr)      
      
      do i=1,nl
         sigma2(i)=1.0/w(i)
      end do

c      
c     extract desired parameters
c
      iparm=0
      do i=1,5
         eb(i)=0.0
         if(ifitp(i).eq.1) then 
            iparm=iparm+1
            p2(iparm)=paramp(i)            
         end if
         ifit(i)=ifitp(i)
         params(i)=paramp(i)
      end do
      np=iparm
      
      alpha=alphap
      dens=densp
      do i=1,nl
        tau(i)=taup(i)
        rho(i)=rhop(i)
        sigma2p(i)=sigma2(i)
      end do

      bfld=bfldp

c
c     no. equations is no. lags - do nlls fit
c
      
      call lmdif1(fcn,nl,np,p2,fvec,tol,info,iwa,wa,lwa)
      ist=info
      
c
c     generate error bars here
c
      call fdjac2(fcn,nl,np,p2,fvec,wa,nl,iflag,0.0e0,wa2)

      do i=1,np
         err=0.0
         do j=1,nl
            err=err+(wa(j+(i-1)*nl))**2
         end do
         if(err.gt.0.0) eb(i)=sqrt(err**-1)
      end do

c     reorder results
      iparm=0
      do i=1,5
         if(ifit(i).eq.1) then 
            iparm=iparm+1
            paramp(i)=p2(iparm)
            ebp(i)=eb(iparm)
         else            
            ebp(i)=0.0
            paramp(i)=params(i)
         end if
      end do

 9    continue
c      write(*,*) dens
c      write(*,*) te
      
      return
      end


      subroutine fcn(m,n,x,fvec,iflag)
c
c     provides m functions (fvec) in n variables (x) to minimize
c     in least-squares sense
c
      
      parameter(nlmax=100,npmax=10)
      integer m,n,iflag
      real fvec(m)
      
      integer ifit(npmax)
      real x(n),fv(nlmax),anorm,params(npmax),ev(nlmax*nlmax)
      real tau(nlmax),rho(nlmax),sigma2(nlmax)
      real chisq,acf

      include 'fitter.h'

      common/fitter/tau,rho,sigma2,params,ifit
      common /errs/ chisq
      common /trans/ev

c     1 0 0 0 0    fit zero lag (otherwise set to default)
c     0 1 0 0 0    fit Te (otherwise default)
c     0 0 1 0 0    fit Ti (otherwise set to Te)
c     0 0 0 1 0    fit H+ (otherwize default)
c     0 0 0 0 1    fit He+ (otherwise default)

 8    continue

c
c     ignore collisions
c
c      write(*,*) "starting fcn"
      ven=0.0
      do i=1,3
         vin(i)=0.0
      end do

      iparm=0

c
c     zero lag normalization constant
c
      if(ifit(1).eq.1) then
         iparm=iparm+1
         c=x(iparm)
      else
         c=params(1)
      end if
c
c     Te
c
      if(ifit(2).eq.1) then
         iparm=iparm+1
         te=x(iparm)
      else
         te=params(2)
      end if
c
c     Ti - default is Te rather than initial value
c
      if(ifit(3).eq.1) then
         iparm=iparm+1
         do i=1,3
            ti(i)=x(iparm)
         end do
      else
         do i=1,3
            ti(i)=te
         end do
         params(3)=te
      end if
  
c     three-ion plasma

      nion=3
c
c     composition H+ first
c
      fi(1)=1.0
      if(ifit(4).eq.1) then
         iparm=iparm+1
         fi(2)=(x(iparm))
         fi(1)=fi(1)-(x(iparm))
      else
         fi(2)=params(4)
         fi(1)=fi(1)-params(4)
      end if
c
c     He+
c
      if(ifit(5).eq.1) then
         iparm=iparm+1
         fi(3)=(x(iparm))
         fi(1)=fi(1)-(x(iparm))
      else
         fi(3)=params(5)
         fi(1)=fi(1)-params(5)
      end if

c      write(*,*) x,bfld,alpha,dens,ak,m
     
      call gaussq(0.0,anorm)
      anorm=anorm/c
 
      if(anorm.eq.0.0.or.m.eq.0)return
      chisq=0.0
      do i=1,m
        call gaussq(tau(i),acf)
        fv(i)=(acf/anorm-rho(i))
      end do

c
c     transform to space where s2 is diagonal
c     (are i and j transposed below?)
c      
      do i=1,m
         fvec(i)=0.0
         do j=1,m
            fvec(i)=fvec(i)+fv(j)*ev(j+(i-1)*m)/sqrt(sigma2(i))
         end do 
         chisq=chisq+fvec(i)**2
      end do

      chisq=chisq/float(m)

c      write(*,*) fvec
c      write(*,*) chisq
c      stop
      
      return
      end

c
c     not really fond of above code
c


      complex function cj_ion(theta,psi)
c     
      real theta,phi,psi,alpha
      complex z
      complex*16 zz
      z=zz(dcmplx(-theta,psi))
      cj_ion=z*cmplx(0.0,-1.0)
      return
      end

      complex function cj_electron(theta,phi,psi,alpha)
c
c     Theta, phi, and psi are the normalized frequency, gyrofrequency,
c     and collision frequency. Alpha is angle between wavevector and
c     magnetic field in radians.
c
      parameter(nterms=10)
      real theta,phi,psi,alpha
      real arg
      complex cj,cy(0:nterms)
      complex z
      complex*16 zz
      integer m,nz,ierr
      integer imode
      common/mode/imode
c
      if(imode.eq.3) then

         arg=0.5*(sin(alpha)/phi)**2

c     calculate modified Bessel functions using amos library
      
         call cbesi(cmplx(arg,0.0),0.0,2,nterms+1,cy,nz,ierr)

         cj=cmplx(0.0,0.0)
     
         do m=-nterms,nterms
            z=zz(dcmplx(-(theta-float(m)*phi)/cos(alpha),
     &           psi/cos(alpha)))
            cj=cj+z*cy(iabs(m))
         end do
         
         cj_electron=cj*cmplx(0.0,-1.0/cos(alpha))

      else

         z=zz(dcmplx(-theta/cos(alpha),psi/cos(alpha)))
         cj_electron=z*cmplx(0.0,-1.0/cos(alpha))
         cj_electron=cj_electron*(1.0-(sin(alpha)/phi)**2/2.0)

      end if

      return
      end

      complex function y_ion(theta,psi)
c
      real theta,phi,psi,alpha
      complex cj, cj_ion
      cj=cj_ion(theta,psi)
      y_ion=cj*cmplx(theta,-psi)+cmplx(0.0,1.0)
      y_ion=y_ion/(1-psi*cj)
      return
      end

      complex function y_electron(theta,phi,psi,alpha)
c
      real theta,phi,psi,alpha
      complex cj, cj_electron
      cj=cj_electron(theta,phi,psi,alpha)
      y_electron=cj*cmplx(theta,-psi)+cmplx(0.0,1.0)
      y_electron=y_electron/(1.0-psi*cj)
      return
      end

      real function spect1(omega)
c
c     Function to generate IS ion-line spectrum 
c     Provisions for nimax ions
c     No provisions for drifts
c
c     imode=1   Farley B-field treatment for electrons 
c           2   use Mike Sulzer's model for ye
c           3   calculate ye using sums of Bessel functions
c
      include 'fitter.h'

      real omega,thetae,thetai(nimax),psie,psii(nimax),phi,p
      real tr(nimax),vti(nimax),vte,omegae
      real bk,em,e,dlf,dl,pi
      real alpha2,densmks,freq
      complex ye,yed,yet,yi(nimax),sum1,sum2,sumdl
      complex y_ion, y_electron, y_esum
      integer i,j,k,imode
      common/mode/ imode
      data bk,em,e,dlf/1.38e-23,9.1e-31,1.6e-19,4772.9/
      
c     fi is ion fraction, wi is ion atomic weight 
c     dens is electron density (cgs), alpha is angle
c     between k and the magnetic field (radians)

      pi=4.0*atan(1.0)

      if(omega.eq.0.0) omega=1.0

      omegae=e*bfld*1.0e-4/em
      densmks=dens*1.0e6

      sum1=0.0
      sum2=0.0
      
      do i=1,nion
        tr(i)=te/ti(i)
        vti(i)=sqrt(bk*ti(i)/(1.67e-27*wi(i)))
        thetai(i)=(omega/ak)/(sqrt(2.0)*vti(i))
        psii(i)=(vin(i)/ak)/(sqrt(2.0)*vti(i))
        yi(i)=y_ion(thetai(i),psii(i))
        sum1=sum1+fi(i)*tr(i)*yi(i)
        sum2=sum2+fi(i)*yi(i)
      end do
      dl=ak**2*dlf*te/densmks
c      write(*,fmt='("Before YE")')
c      write(*,*) imode
c      call exit

      if(imode.eq.1.or.imode.eq.3) then
      
         vte=sqrt(bk*te/em)
         thetae=(omega/ak)/(sqrt(2.0)*vte)
         phi=(omegae/ak)/(sqrt(2.0)*vte)
         psie=(ven/ak)/(sqrt(2.0)*vte)
         ye=y_electron(thetae,phi,psie,alpha)
	 write(*,fmt='("AFTER YE")')
c         call exit
  
      else if(imode.eq.2) then
c
c     use Mike Sulzer's library here: alpha2 is angle off perp in degrees
c
         
         freq=omega/(2.0*pi)
         alpha2=abs(pi/2.0-alpha)*180.0/pi
c         write(*,*) "ye: ", ye
         call collision(densmks, te, freq, alpha2, ye)
c         write(*,fmt='("AFTER COLLISION")')
c         call exit
         ye=ye*omega+cmplx(0.0,1.0)        

      end if
      yed=ye+cmplx(0.0,dl)  
	
      p=(cabs(ye))**2*real(sum2)+cabs(sum1+cmplx(0.0,dl))**2*real(ye)
      p=p/(cabs(yed+sum1))**2
      spect1=p*2.0e0/(omega*pi)
      write(*,*) "spect1:",spect1
      return
      end

      subroutine acf2(wl, tau, te1, ti1, fi1, ven1, vin1, wi1, nion1,
     &     alpha1, dens1, bfld1, acf)
c
c     computes autocorrelation function for given plasma parameters
c     by integrating real spectrum
c     tau in sec., alpha in radians, density in cgs, bfield in cgs
c     scattering wavelength (wl) in meters
c
      include 'fitter.h'

      real wl,tau,te1,ti1(nion1),fi1(nion1),ven1,vin1(nion1),alpha1,
     &     dens1,bfld1,acf
      real pi
      integer nion1
      integer wi1(nion1)
      integer i,j,k,imode
      common /mode/imode
c
c      write(*,*) "INITIAL acf:",wl,tau,te1,ti1,fi1,ven1,vin1,wi1,alpha1
      write(*,*) "INITIAL acf:",dens1, bfld1, acf, nion1
c      write(*,fmt='("INIT")')
      pi=4.0*atan(1.0)
c
c     copy arguments to common block
c
      ak=2.0*pi/wl
      imode=2
      write(*,*) "imode:",imode

      nion=nion1
      alpha=alpha1
      te=te1
      ven=ven1
c      write(*,fmt='("INIT2")')
      do i=1,nion
c      write(*,fmt='("INIT2.5")')
        ti(i)=ti1(i)
        fi(i)=fi1(i)
        vin(i)=vin1(i)
        wi(i)=wi1(i)
      end do
c      write(*,fmt='("INIT3")')
      dens=dens1
      bfld=bfld1

c      write(*,*) wl,alpha1,bfld1,dens
c      call exit
c      write(*,fmt='("Before Gauss")')
      call gaussq(tau,acf)

      write(*,*) "FINAL acf:",acf

c      write(*,fmt='("After Gauss")')
      return
      end

      subroutine gaussq(tau,acf)
c
c     Computes cosine or sine transform of given real function
c     Uses quadpack, tau in sec.
c     
      real a,abserr,epsabs,acf,tau,work,chebmo
      integer ier,integr,iwork,last,leniw,lenw,limit,limlst,
     *     lst,maxp1,neval
      integer knorm,ksave,momcom,nrmom
      dimension iwork(300),work(1625),chebmo(61,25)
      external spect1
      include 'fitter.h'
c
c      write(*,*) "inside gaussq"
      a = 0.0
      b = 2.5e4*ak  ! upper integration limit open to debate (was 1.5, now 2.5)
      integr = 1
      epsabs = 1.0e-4
      limlst = 50
      limit = 100
      leniw = limit*2+limlst
      maxp1 = 61
      lenw = leniw*2+maxp1*25

c      write(*,*) leniw,lenw
      
      nrmom=0
      ksave=0
      momcom=0
      write(*,*) "Before qc25f:",acf,imode
c     much faster, more robust
c      write(*,*) "acf_in: ",acf
      call qc25f(spect1,a,b,tau,integr,nrmom,maxp1,ksave,acf,
     &     abserr,neval,resabs,resasc,momcom,chebmo)
c      write(*,*) "acf_out: ",acf
c      call qawf(spect1,a,tau,integr,epsabs,acf,abserr,neval,
c     &     ier,limlst,lst,leniw,maxp1,lenw,iwork,work)
      write(*,*) "After qc25f:",acf
c      call exit
      return
      end