c
c     $Id: full_profile.f,v 1.3 2013/09/12 19:11:48 daveh Exp daveh $ DLH 
c      

      subroutine profile(acf_sum,acf_err,power,en,alag,thb2,bfm2,
     & ote,ete,oti,eti,oph,eph,ophe,ephe,range2,ut,
     &     NHTS,NACF,IBITS,acf_avg_real,status)
c
      complex acf_sum(4,NHTS,IBITS)
      real acf_avg_real(NHTS,IBITS)
      real acf_err(NHTS,IBITS),range2(NHTS)
      real power(NHTS),en(NHTS),alag(IBITS),thb2(NHTS),bfm2(NHTS)
      real ote(NHTS),ete(NHTS),oti(NHTS),eti(NHTS),oph(NHTS),eph(NHTS)
      real ophe(NHTS),ephe(NHTS)
      integer NHTS,NACF,IBITS
      real status(1)
      common /chisq/chi2
c
c     arguments above, internal variables below: have to make these match
c
      include 'fpa.h'
      include 'bfield.h'
c
      common /mode/imode
      common /utime/uttime
c
      include 'spline.h'
      include 'fitter.h'
c
c     functions to minimize are all lagged profiles plus penalties
c     variables to minimize them with are all spline values
c

      parameter(nfun=nrange*nl+nstate*(nspline-2)+npen)
      parameter(nvar=nspline*nstate)
      parameter(lwa=nfun*nvar+nstate*nvar+nfun)
      real fvec(nfun),wa(lwa),xe(nvar)
      integer iwa(nvar)
c      integer lwa
c
      external fcn_lpreg
c
      pi=4.0*atan(1.0)
      uttime=ut
c      write(*,*) "ut: ",ut
c      write(*,*) "NHTS: ",NHTS
c      write(*,*) "NACF: ",NACF
c      write(*,*) "IBITS: ",IBITS
c      call exit
c      write(*,*) "acf_sum: ", acf_sum(1,111,11)      
c      write(*,*) "Starting Profile Process"
c     set up basic physical parameters
c      dens=639747.44
c      dens=509664.781
c      te=3971.98926
c      write(*,*) "acf_sum: ", acf_sum(1,31,2)
c      call    exit

      wl=3.0
      ak=2.0*pi/wl

c     analysis starts at gate 15

      ir0=15

      r0=range2(ir0)
      dr=range2(ir0+1)-range2(ir0)

      do i=1,nfield ! nrange+10   memory hole
         bfld_prof(i)=bfm2(i+ir0-1)
         alpha_prof(i)=thb2(i+ir0-1)*pi/180.0
      end do

      do i=1,nrange+nl ! careful ... check nfield compared to nrange+nl
         densp(i)=en(i+ir0-1)
         do j=1,nl
c            write(*,*) acf_sum(1,i+ir0-1,j)
            plag(j,i)=acf_sum(1,i+ir0-1,j)
            
            plag_errors(j,i)=acf_err(i+ir0-1,j)
c            write(*,*) j,i,plag_errors(j,i)
         end do
      end do      
c      call exit
c     set some fitting conventions
c      write(*,*) "plag: ",plag(2,2)
      imode=2
      nion=3

      ven=0.0
      do i=1,nion
         vin(i)=0.0
      end do
      wi(1)=16
      wi(2)=1
      wi(3)=4

c     system constant

      sconst=1.0

c     set up splines: DLH need to modify for different nu. splines

      do i=1,nspline+norder
         ta(i)=r0+float(i-1-2)*delr ! 2-knot offset for 4th order spline
      end do

c     initial guess - set density, He+ profiles here
         
      do i=1,nspline ! DLH modify for different no. splines 
         bcoef(i,1)=log10(abs(densp(1+int(float(i-1)*delr/dr))))
                                ! need good guess for density profile 
         tmp=0.03
         bcoef(i,5)=atanh(2.0*tmp-1.0) ! not so critical
      end do        

c      write(*,*) bcoef
c      call exit

c     grid search for other params to help set initial guess
c     could include Ne, He+ in grid search too
c      write(*,fmt='("before GRID")')
      call grid()
c      write(*,fmt='("After GRID")')
      write(*,*) "chi2 after grid:",chi2
c      write(*,fmt='("Just After chi")')
c      call exit
      if(chi2.le.3.5) then

c      write(*,fmt='("After chi")')
c      call exit
      tol=1.0e-9

c      call lmdif1(fcn_lpreg,nfun-nl*nl,nvar,bcoef,fvec,tol,
c     &     info,iwa,wa,lwa)

      info=0
      m=nfun-nl*nl
      n=nvar
      zero=0.0
      factor=1.0e-2   ! rather small initial step size

      maxfev = 200*(n + 1)
      ftol = tol
      xtol = tol
      gtol = zero
      epsfcn = zero
      mode = 1
      nprint = 0
      mp5n = m + 5*n
c      write(*,*) wa(5*n+1)
c      call exit
      write(*,*) "before densp:", densp(1),bcoef(1,2)
      call lmdif(fcn_lpreg,m,n,bcoef,fvec,ftol,xtol,gtol,
     &     maxfev,epsfcn,wa(1),
     *           mode,factor,nprint,info,nfev,wa(mp5n+1),m,iwa,
     *           wa(n+1),wa(2*n+1),wa(3*n+1),wa(4*n+1),wa(5*n+1))

c      write(*,fmt='("After lmdif")')
c      write(*,*) info
      write(*,*) "densp_after:", densp(1),bcoef(1,2)
c      call exit

c
c     error handling; evalate solution point first, then Jacobian
c
      iflag=0
c      write(*,fmt='("before fcn_lpreg")')
      call fcn_lpreg(m,n,bcoef,fvec,iflag) ! common sys constant throughout
c      write(*,fmt='("After fcn_lpreg")')
      nm1=n*m+1
      call fdjac2(fcn_lpreg,m,n,bcoef,fvec,wa,m,iflag,epsfcn,wa(nm1))

      do i=1,n
         err=0.0
         xe(i)=0.0
         do j=1,m-npen
c         do j=2+(nspline-2)*nstate,m-npen ! skip equations unrelated to data
            err=err+MIN(wa(j+(i-1)*m)**2,10.0)
         end do
         if(err.gt.0.0) xe(i)=sqrt(err**(-1))
      end do
    
      call propagate(xe)

      do i=1,nrange
         ete(i+ir0)=etep(i)
         eti(i+ir0)=etip(i)
         eph(i+ir0)=ehfp(i)
         ephe(i+ir0)=ehefp(i)
      end do
      
c
c     output parameters
c

      do i=1,nrange
         altp(i)=r0+dr*float(i-1)
         call get_spline(altp(i),densp(i),tep(i),tip(i),
     &        hfp(i),hefp(i))
         ote(i+ir0)=MAX(tep(i),tip(i))
         oti(i+ir0)=tip(i)
         oph(i+ir0)=hfp(i)
         ophe(i+ir0)=hefp(i)
      end do
c
c     theoretical lag products and other quantities for plotting
c

      call lagp(plag,wl,r0,dr,nl,nrange)
      write(*,*) "plag_plag",REAL(plag(12,58))
      write(*,*) "plag_plag",REAL(plag(12,72))
c      write(*,*) "plag_plag",AIMAG(plag(12,72))
      do i=1,nrange                        !
         en(i+ir0-1)=densp(i)
         do j=1,nl
            acf_avg_real(i+ir0-1,j)=REAL(plag(j,i))
c            acf_avg_imag(i+ir0-1,j)=AIMAG(acf_sum(1,i+ir0-1,j))
         end do
      end do      
      write(*,*) "acf",acf_avg_real(72,12)
c      write(*,*) "acf",acf_avg_imag(72,12)
      end if
      status(1)=chi2
      
c      write(*,*) "status: ",status

 5    return
      end


      subroutine grid()
c
c     minimize chi-squared statistic through grid search
c     just optmize temperature and H+ profiles
c
      include 'spline.h'
      include 'fpa.h'
      include 'fitter.h'
c
      real plag2(nl,nrange)
      real fvec(nrange*nl)
      real chmin,uttime
      integer lmin,mmin,nmin
      integer iloop,iskip,l1,l2,m1,m2,n1,n2
c
      common /chisq/chi2
      common /utime/uttime
c
      chmin=1.0e10

      do i=1,nrange
         altp(i)=r0+dr*float(i-1)
      end do

      goto 3
      call newplot(.true.)
      call graphinit()
      call erase

      call setorigin(6.5,3.5,90.0)
      call linearaxis("L",altp(1),altp(nrange),1.5,5.0,2,0.0,50.0)
      call linearaxis("B",0.0,5000.0,1.5,5.0,2,0.0,500.0)
      call drawaxis("T",2)
      call drawaxis("R",2)
      call labelaxis("B","(f5.0)",.1,2,0.0,2000.0)
      call labelaxis("L","(f5.0)",.1,2,0.0,100.0)
      call ctitle("L","Altitude (km)",.125,2)
      call ctitle("B","T\\se\\u (K)",.125,2)

      call setorigin(6.5,7.5,90.0)
      call linearaxis("L",altp(1),altp(nrange),1.5,5.0,2,0.0,50.0)
      call linearaxis("B",0.0,1.0,1.5,5.0,2,0.0,0.25)
      call drawaxis("T",2)
      call drawaxis("R",2)
      call labelaxis("B","(f5.1)",.1,2,0.0,0.5)
      call labelaxis("L","(f5.0)",.1,2,0.0,100.0)
      call ctitle("B","H\\S+\\U frac.",.125,2)

 3    continue

      do iloop=1,2
         if(iloop.eq.1) then
            l1=2
            l2=26
            m1=2
            m2=10
            n1=2
            n2=16
            iskip=2
         else
            l1=lmin-2
            l2=lmin+2
            m1=mmin-2
            m2=mmin+2
            n1=nmin-2
            n2=nmin+2
            iskip=1
         end if

      do l=l1,l2,iskip ! temperature search
         do k=1,nspline

            tmp=-2.5+float(k-1)*float(l-1)*1.3e-2 
c            if(uttime.ge.2.0.and.uttime.lt.10.0) then
c               tmp=-1.5+float(k-1)*float(l-1)*8.0e-3 
c            else
c               tmp = -0.75 -0.2*float(l) + 3.5e-1*float(k-1) ! was 3.5
c               tmp=-1.5 + 3.0e-2*float(l-1) 
c     &              *(float(k-1)**3.0/float(nspline-1)**2.0)
c            end if

            bcoef(k,2)=tmp
            if(uttime.ge.10.0.and.uttime.le.14.0) then ! was 10.5
               bcoef(k,3)=tmp-0.1 
            else
               bcoef(k,3)=tmp-0.02
            end if
         end do

         do m=m1,m2,iskip ! H+ scale height search DLH 2015
            
            do n=n1,n2,iskip ! H+ crossing  height search
	      
               do k=1,nspline
                  tmp=float(k-nspline/4-n)*(float(m-1)*5.0e-2+0.1)
                  bcoef(k,4)=tmp
               end do
               
               do i=1,nrange
                  call get_spline(altp(i),densp(i),tep(i),tip(i),
     &                 hfp(i),hefp(i))

               end do
	       
               goto 4
               call setorigin(6.5,3.5,90.0)
               call setscale("L",altp(1),altp(nrange),1)
               call setscale("B",0.0,5000.0,1)
               call setclip("L",250.0,1500.0)
               call plotline("B,L",tep,altp,nrange,2,0)
              
               if(l.eq.1) then
                  call setorigin(6.5,7.5,90.0)
                  call setscale("L",altp(1),altp(nrange),1)
                  call setscale("B",0.0,1.0,1)
                  call setclip("L",250.0,1500.0)
                  call plotline("B,L",hfp,altp,nrange,2,0)
               end if
 4             continue
	       
c     calculate chi-2
c               write(*,*) "wl: ",wl, "r0: ",r0
c               write(*,*) "dr: ",dr,"nl: ",n,"nrange", nrange
c	       write(*,*) "before lagp inside grid"
               call lagp(plag2,wl,r0,dr,nl,nrange)
c               write(*,*) "plag2: ",plag2(12,72)
c               write(*,*) "plag2: ",plag2
c               call exit	       
               call get_scale(plag2)
c               write(*,*) "plag: ",plag(13,70)
c               write(*,*) "plag2: ",plag2(13,70)
c	       write(*,*) "plag_errors: ",plag_errors(13,70)
c	       write(*,*) "sconst: ",sconst              
c               call exit
c               write(*,*) "cal chi2.1: ", chi2
c               call exit
               chi2=0.0
               ipt=0
c               write(*,*) "cal chi2.1: ", chi2
               do j=nl+1,nrange
                  do i=1,nl
                     ipt=ipt+1
                     fvec(ipt)=(plag(i,j)-plag2(i,j)*sconst)/
     &                    plag_errors(i,j)
c                     write(*,*) j,i,plag(i,j),plag2(i,j),
c     &                    plag_errors(i,j),sconst
c		     call exit
                     chi2=chi2+fvec(ipt)**2
c                     write(*,*) "cal chi2.1: ", chi2
                  end do
c                  call exit
               end do
c               write(*,*) "cal chi2.2: ", chi2
c               call exit
               chi2=sqrt(chi2/float(nrange*nl))!+float(l)*0.005 ! DLH??
               if(chi2.lt.chmin) then
                  lmin=l
                  mmin=m
                  nmin=n
                  chmin=chi2
               end if

               write(*,*) "A",l,m,n,chi2,chmin,"B" ! best
c               call exit
c	        write(*,fmt='("intermediate")')
            end do
         end do
      end do

      end do
c      write(*,fmt='("befire flush")')
c      call where()
      write(22,*) lmin,mmin,nmin,chmin
      call flush(22)
c      write(*,fmt='("after flush")')
c     best result here

      do k=1,nspline
         tmp=-2.5+float(k-1)*float(lmin-1)*1.3e-2 
c         if(uttime.ge.2.0.and.uttime.lt.10.0) then
c         tmp=-1.5+float(k-1)*float(lmin-1)*8.0e-3 
c         else
c            tmp = -0.75 -0.2*float(l) + 3.5e-1*float(k-1) ! was 3.5
c         tmp=-1.5 + 3.0e-2*float(lmin-1) ! new family of temps 
c     &     *(float(k-1)**3.0/float(nspline-1)**2.0)
c         end if
         bcoef(k,2)=tmp
         if(uttime.ge.10.0.and.uttime.le.14.0) then  ! was 10.5
            bcoef(k,3)=tmp-0.1
         else
            bcoef(k,3)=tmp-0.02
         end if
         tmp=float(k-nspline/4-nmin)*(float(mmin-1)*5.0e-2+0.1)
         bcoef(k,4)=tmp
      end do

c     reset system constant

      do i=1,nrange
         call get_spline(altp(i),densp(i),tep(i),tip(i),
     &        hfp(i),hefp(i))
      end do
c      write(*,fmt='("------------before lagp------------")')
c      write(*,*) chi2
c      call exit
      call lagp(plag2,wl,r0,dr,nl,nrange)   
c      write(*,fmt='("after lagp")')      
      call get_scale(plag2)
c      write(*,*) chi2
      chi2=chmin
      return
      end

      subroutine propagate(xe)
c
c     propagate errors through splines
c
      include 'fpa.h'
      include 'spline.h'
      parameter(nvar=nspline*nstate)
      real xe(nvar)


      do i=1,nrange

         write(*,*) "densp:", densp(i),bcoef(1,2)
c         call exit
c                  
         edensp(i)=bvalue(ta,xe(1+0*nspline),nspline,norder,altp(i),0)
     &        *densp(i)*2.30
         x=bvalue(ta,bcoef(1,2),nspline,norder,altp(i),0)
         etep(i)=bvalue(ta,xe(1+1*nspline),nspline,norder,altp(i),0)
     &        *(t1/2.0)/cosh(x)**2 ! DLH 10/14 was 3500

         write(*,*) i,altp(i),x,
     &        bvalue(ta,xe(1+1*nspline),nspline,norder,altp(i),0),
     &        (3500.0/2.0)/cosh(x)**2,xe(1+i/4+1*nspline)

         x=bvalue(ta,bcoef(1,3),nspline,norder,altp(i),0)
         etip(i)=bvalue(ta,xe(1+2*nspline),nspline,norder,altp(i),0)
     &        *(t1/2.0)/cosh(x)**2  ! DLH 10/14 was 3500
         x=bvalue(ta,bcoef(1,4),nspline,norder,altp(i),0)
         ehfp(i)=bvalue(ta,xe(1+3*nspline),nspline,norder,altp(i),0)
     &        *(1.0/2.0)/cosh(x)**2
         x=bvalue(ta,bcoef(1,5),nspline,norder,altp(i),0)
         ehefp(i)=bvalue(ta,xe(1+4*nspline),nspline,norder,altp(i),0)
     &        *(1.0/2.0)/cosh(x)**2

         write(*,*) altp(i),edensp(i),etep(i),etip(i),ehfp(i),ehefp(i)
c         call exit
      end do

c      stop

      return
      end

      subroutine fcn_lpreg(m,n,x,fvec,iflag)
c
      include 'spline.h'
      include 'fpa.h'
      include 'fitter.h'
c
      integer m,n,iflag,iflaglast
      real x(n),fvec(m),pen(nstate+npen)
      real plag2(nl,nrange)
      real regu(nstate+npen)
      data regu/1.0e1,5.0e1,10.0e1,1.0e2,5.0e1,
     &     1.0e1,2.0e0,1.0e0,5.0e-3,1.0e0/ 

c     2019 
      common /iflaglast/iflaglast
      common /chisq/chi2
      common /utime/uttime
c
c      write(*,*)uttime," 3 *"

c     zero out, copy over some arrays

c      write(*,*) "Starting fcn_lpreg"
      write(*,*) "x:",x(2)

      do j=1,nstate+npen
         pen(j)=0.0
      end do

      ipt=0
      do j=1,nstate
         do i=1,nspline
            ipt=ipt+1
            bcoef(i,j)=x(ipt)
c            write(*,*)j,i,bcoef(i,j)
         end do
      end do

c     populate profiles from splines
      write(*,*) "densp1:",densp(1)
      do i=1,nrange
         altp(i)=r0+dr*float(i-1)
c         write(*,*) "densp1:",densp(i)
         call get_spline(altp(i),densp(i),tep(i),tip(i),
     &        hfp(i),hefp(i))
c         write(*,*) "densp2:",densp(i)
c         call exit
      end do
      write(*,*) "densp2:",densp(1)
c     regularize spline coefficients, establish boundary conditions
c      write(*,*) "After spline_fcn_lpreg"
      ipt=0
      do i=2,nspline-1
         do j=1,nstate

            fac=1.2  ! DLH 7/14 could be bendy (or else get spurious He+)
            if((uttime.lt.10.0.or.uttime.gt.15.0).and.j.eq.2) fac=4.0
            
            ipt=ipt+1
            tmp=(bcoef(i-1,j)-2.*bcoef(i,j)+bcoef(i+1,j))
            fvec(ipt)=tmp*regu(j)*fac
            pen(j)=pen(j)+fvec(ipt)**2
         end do
      end do

c     rescale for clearer diagnostic (just for printing)

      do i=1,nstate
         pen(i)=sqrt(pen(i)/float(nspline)) 
      end do
       
c      write(*,*) x

c     compute model prediction error to minimize in LS sense
c      write(*,*) "before lagp_fcn_lpreg"
      call lagp(plag2,wl,r0,dr,nl,nrange)
c      write(*,*) "before lagp_fcn_lpreg"
c     adjust system constant
      write(*,*) "flgas:",iflag,iflaglast
      if(iflag.eq.1.and.iflaglast.eq.2) then
         sclast=sconst
         call get_scale(plag2)
         sconst=0.1*sconst+0.9*sclast
      end if
      iflaglast=iflag
         
      chi2=0.0
      do j=nl+1,nrange
         do i=1,nl
            ipt=ipt+1
            fvec(ipt)=(plag(i,j)-plag2(i,j)*sconst)/plag_errors(i,j)
            chi2=chi2+fvec(ipt)**2
         end do
      end do
      chi2=sqrt(chi2/float(nrange*nl))

c     add penalties for Tr, composition, He+, lower boundary values, etc.
c     have to adjust penalties for different conditions (e.g. sunrise, night)

      do i=1,nrange
         pen(nstate+1)=pen(nstate+1)+hefp(i)**2
         
         if(tep(i).gt.tip(i)) then ! 3/2015 500 -> 1000 ?
            if(uttime.ge.10.0.and.uttime.le.15.0) then ! was 10.5, 14->15
               scale=500.0  ! dlh 7/14 was 1500, but this defeated He+, so 500
            else           
               scale=50.0  
            end if
         else
            scale=10.0           
         end if 

         pen(nstate+2)=pen(nstate+2)+((tip(i)-tep(i))/scale)**2
         pen(nstate+3)=pen(nstate+3)+exp((hfp(i)+hefp(i)-1.0)/2.0)**2

      end do

      pen(nstate+4)=(tip(1)-700.0)**2*float(nrange) ! or 800?
      pen(nstate+5)=(hfp(nrange)-1.0)**2*float(nrange)

      do i=1,npen
         tmp=sqrt(pen(nstate+i))*regu(nstate+i)
         ipt=ipt+1
         fvec(ipt)=tmp
         pen(nstate+i)=tmp/sqrt(float(nrange))
      end do
      
c      write(*,*) "chi2:",chi2," sys const:",sconst,
c     &     " flag:",iflag," penalties:"
c      write(*,*)(pen(i),i=1,5)
c      write(*,*)(pen(i),i=6,10)
c      write(*,*) "ut:",uttime
c      write(*,*) fvec(222)

      return
      end

      subroutine get_scale(plag2)     
c
      include 'fpa.h'
      real plag2(nl,nrange)
c      
      a=0.
      b=0.
c      do j=3*nrange/4,nrange
      do j=nl+1,nrange 
c         write(*,*) j,plag(1,j),plag2(1,j)
         a=a+plag(1,j)*plag2(1,j)/plag_errors(1,j)**2
         b=b+plag(1,j)*plag(1,j)/plag_errors(1,j)**2
      end do

      sconst=b/a
      
      write(*,*)"sconst_inside_get_scape:",sconst,a,b
c      stop
      
      return
      end