C ------------------------------------------------------------------------
C       PROGRAM Field: T(x) version      
C       DEDUCE CONTINUOUS THERMODYNAMIC FIELD FROM PARTICLE DATA USING 
C       MAXIMUM LIKELIHOOD INFERENCE, ASSUMING QUASI-EQUILIBRIUM PHASE 
C       SPACE DISTRIBUTION.
C
C                                       JAN.30.1997
C                                       DEVELOPED BY JU LI AT MIT
C ------------------------------------------------------------------------
            
      PROGRAM FIELD_TX
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MAX_NP=50000,MAX_ORDER=30,MESH=500,NDIM=3)
      DIMENSION VV(MAX_NP*NDIM+1),X(MAX_NP),V2(MAX_NP),
     A          COEF_TX(MAX_ORDER+1)
      DATA LP_TEMP,LP_PAR /25,26/
      INTEGER ORDER_TX
      CHARACTER *70 BUF
      
c     =====================================================================
c     Read instructions from input file "con":
      
      READ (*,'(A70)') buf
C     "Input control parameter 1:"
      read *, par1

      READ (*,'(/,A70)') buf
c     "Input control parameter 2:"
      read *, par2
      
      READ (*,'(/,A70)') buf
c     "Input number of particles:"
      read *, np

      READ (*,'(/,A70)') buf
c     "Highest order of Chebyshev expansion:"
      read *, ORDER_TX

      READ (*,'(/,A70)') buf
c     "Give the random number generator a kick:"
      read *, KICK
      do 35 I=1,KICK
 35   a = randsafe()

      READ (*,'(/,A70)') buf
c     "Fit particle data from ** to **:"
      read *, bound_low, bound_high
c     ======================================================================
      
      xmin = 0.
      xmax = 1.
      PI = dacos(-1.d0)

c     generate random number series NDIM*NP long
      ip = 1 
 78   s1=dsqrt(-2.d0*log(randsafe()))
      s2=2*PI*ran1(KICK)
      VV(ip)=s1*cos(s2)
      VV(ip+1)=s1*sin(s2)
      ip = ip+2
      if (ip.le.NDIM*NP) goto 78
      
c     convert it to V2 series NP long
      ip = 0
      do 34 i=1,NP
      v2(i) = 0.
      do 56 j=1,NDIM
 56   v2(i) = v2(i)+VV(ip+j)**2
      ip = ip+NDIM
 34   enddo
      
      do i=1,NP
c     generate the particle coordinates,
      x(i) = bound_low + (bound_high-bound_low)*randsafe()
c     and then we scale the particle velocities according 
c     to prescribed field
      v2(i) = v2(i) * temp(x(i),par1,par2,xmin,xmax)
      enddo

c     write the particle data into file "par.out"
      OPEN (UNIT=LP_PAR, STATUS='UNKNOWN', FORM='FORMATTED',
     a      FILE='par.out')      
      DO 64 I=1,NP
 64   WRITE (LP_PAR,'(2(1X,F10.5))'),x(i),v2(i)/NDIM
      CLOSE (LP_PAR)
      
      call fit_TX (np,x,v2,ORDER_TX,COEF_TX,NDIM,xmin,xmax)     
c     get the temperature at two ends and compare fitting with original
      T_xmin = get_TX (xmin,ORDER_TX,COEF_TX,xmin,xmax)
      T_xmax = get_TX (xmax,ORDER_TX,COEF_TX,xmin,xmax)
      print *,'T_Xmin: Original = ', temp(xmin,par1,par2,xmin,xmax),
     a        'Fitted = ', T_xmin
      print *,'T_Xmax: Original = ', temp(xmax,par1,par2,xmin,xmax),
     a        'Fitted = ', T_xmax

      
      OPEN (UNIT=LP_TEMP, STATUS='UNKNOWN', FORM='FORMATTED',
     a      FILE='temp.out')
c     give the average drift and fluctuation of the temperature field
      XDEL = (XMAX-XMIN)/(MESH-1)
      xx = xmin
      T_sum = 0.
      T_var = 0.
      do i = 1,mesh
      T_orig = temp(xx,par1,par2,xmin,xmax)
      T_fit  = get_TX(xx,ORDER_TX,COEF_TX,xmin,xmax)
      WRITE (LP_TEMP,'(3(1X,F10.5))') xx,T_orig,T_fit
      diff = T_fit - T_orig
      T_sum = T_sum + diff
      T_var = T_var + diff*diff
      xx = xx+XDEL
      enddo
      CLOSE (LP_TEMP)
      
      print *,'The average drift in T is',T_sum/mesh
      print *,'The average variation in T is',sqrt(T_var/mesh)   
      
      STOP
      END
      
      
      FUNCTION temp (x,par1,par2,xmin,xmax)           
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)      
      xv = (x-xmin)/(xmax-xmin)
c     <Linear>
      temp = par1+(par2-par1)*xv
c     <Quadratic>
c     temp=par1+(par2-par1)*xv+par1*(x-xmax)*(x-xmin)
c     <Cubic>
c     temp=par1-par2*(x-xmax)*(x-xmin)*x+(par2-par1)*xv
c     <Quadruple>
c     temp=par1+par2*x*x*x*x+(par2-par1)*xv
c     <Sinusoidal> --- one period
c     temp=par1+par2*sin(2*dacos(-1.d0)*xv)
c     <Sinusoidal> --- two period
c     temp=par1+par2*sin(4*dacos(-1.d0)*xv)
c     <Step function>
c     if (x.lt.xmin+(xmax-xmin)*par2/(par1+par2)) then 
c     temp = par1
c     else
c     temp = par2
c     endif
c     <Bilinear>
c     if (xv.lt.0.5) then 
c     temp = par1+(par2-par1)*xv*2
c     else
c     temp = par1+(par2-par1)*(1-xv)*2
c     endif
c     <BiCubic>
c     if (xv.gt.0.5) xv = 1-xv
c     temp = par1-par2*xv*(xv-0.5)*(xv-0.3)
c     <Trilinear>
c     xv = mod(xv+0.25,1)
c     if (xv.lt.0.5) then 
c     temp = par1+(par2-par1)*xv*2
c     else
c     temp = par1+(par2-par1)*(1-xv)*2
c     endif
c     <TriCubic>
c     xv = mod(xv+0.25,1)
c     if (xv.gt.0.5) xv = 1-xv
c     temp = par1-par2*xv*(xv-0.5)*(xv-0.3)

      RETURN
      END

      
      SUBROUTINE fit_TX (NN,X,V,ORDER_TX,COEF,NDD,xmin,xmax)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)    
      PARAMETER (MAX_NP=50000,MAX_ORDER=30,MESH=500,NDIM=3)
      PARAMETER (SAMP_MIN = -0.95, SAMP_MAX = 0.95)
      COMMON /SHARE1/ TX(MAX_ORDER+1,MAX_NP),BETA(MAX_NP),V2(MAX_NP)
      COMMON /SHARE2/ NP,ORDER,ND
      DIMENSION COEF(MAX_ORDER+1),X(MAX_NP),V(MAX_NP)
      INTEGER ORDER_TX,ORDER
      DOUBLE PRECISION MU(MAX_NP),WEIGHT(MAX_NP),SUM(MAX_ORDER+1)
      
c     pass the SHARE information to common blocks 
      NP = NN
      ND = NDD
      ORDER = ORDER_TX+1
      DO I=1,NP
      V2(I) = V(I)
      ENDDO
      
      DO 532 I=1,NP
c     the range of Chebyshev polynomials are from -1 to 1
      MU(I) = -1 + 2*(X(I)-XMIN)/(XMAX-XMIN)
      WEIGHT(I) = 1/SQRT(1-MU(I)*MU(I))
      TX(1,I) = 1
      TX(2,I) = MU(I)
      DO 532 L=3,ORDER
      TX(L,I) = 2*MU(I)*TX(L-1,I)-TX(L-2,I)
 532  CONTINUE
      
c     Assign initial values to the expansion coefficients:
c     use the orthogonality formula with probalistic interpretation. 
      DO 653 L=1,ORDER
 653  SUM(L) = 0.
c     sum is the slate of inner products 
      NSUM = 0
c     nsum is the number of particles within sampling range
      DO 564 I=1,NP
      IF (MU(I).GT.SAMP_MIN.AND.MU(I).LT.SAMP_MAX) THEN 
      NSUM = NSUM+1
      DO 678 L=1,ORDER
 678  SUM(L) = SUM(L)+ND/V2(I)*TX(L,I)*WEIGHT(I)
      ENDIF
 564  CONTINUE
      PI = DACOS(-1.D0)
      COEF(1) = SUM(1)/NSUM*2/PI
      DO 435 L=2,ORDER
 435  COEF(L) = SUM(L)/NSUM/PI

c     invoke conjugate gradient minimization procedure
      CALL FRPRMN(COEF,ORDER,1D-8,ITER,FMIN)
      print *,' '
      print *,'MINIMIZATION CONVERGED AFTER',ITER,' ITERATIONS.'
      print *,' '
      
      RETURN
      END
      
      
      FUNCTION VALUE(COEF)
c     ***************************************
c     The objective function of minimization
c     ***************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
      PARAMETER (MAX_NP=50000,MAX_ORDER=30,MESH=500,NDIM=3)
      PARAMETER (TINY=1D-20,HUGE=1D+20)
      COMMON /SHARE1/ TX(MAX_ORDER+1,MAX_NP),BETA(MAX_NP),V2(MAX_NP)
      COMMON /SHARE2/ NP,ORDER,ND
      DIMENSION COEF(MAX_ORDER+1)
      INTEGER ORDER
      SAVE VALUE_LAST
      
      SUM = 0.
      DO 721 I=1,NP
      BETA(I) = 0.
      DO 722 L=1,ORDER
      BETA(I) = BETA(I)+COEF(L)*TX(L,I)
 722  CONTINUE
      SUM = SUM+BETA(I)*V2(I)
c     To prevent negative temperature
      IF (BETA(I).LE.0) THEN
      VALUE = VALUE_LAST*(1+SIGN(BETA(I)**2,VALUE_LAST))
      GOTO 564
      ENDIF
 721  CONTINUE

      PRODUCT = 1.
      DO I=1,NP
      PRODUCT = BETA(I)*PRODUCT
      IF (PRODUCT.le.TINY.or.PRODUCT.gt.HUGE) then 
      SUM = SUM - ND*LOG(PRODUCT)
      PRODUCT = 1.
      ENDIF
      ENDDO
      SUM = SUM - ND*LOG(PRODUCT)
      VALUE = SUM/NP
564   VALUE_LAST = VALUE
      RETURN
      END
      
      
      SUBROUTINE FORCE(COEF,F)
c     **************************************
c     force field of the objective function
c     **************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
      PARAMETER (MAX_NP=50000,MAX_ORDER=30,MESH=500,NDIM=3)
      COMMON /SHARE1/ TX(MAX_ORDER+1,MAX_NP),BETA(MAX_NP),V2(MAX_NP)
      COMMON /SHARE2/ NP,ORDER,ND
      DIMENSION COEF(MAX_ORDER+1),F(MAX_ORDER+1),AA(MAX_NP)
      INTEGER ORDER
      
c     Before calling force(), it should have called value(), so
c     there is no need to re-evaluate BETA.
      print *,'evaluating force'
      
      DO I=1,NP
      AA(I) = (V2(I)-ND/BETA(I))/NP
      ENDDO
      
      DO 439 L=1,ORDER
      F(L) = 0.
      DO 439 I=1,NP
      F(L) = F(L)-TX(L,I)*AA(I)
 439  CONTINUE
      
      RETURN
      END
      
      
      FUNCTION GET_TX (X,ORDER_TX,COEF_TX,xmin,xmax)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)      
      PARAMETER (MAX_NP=50000,MAX_ORDER=30,MESH=500,NDIM=3)
      DIMENSION COEF_TX(MAX_ORDER+1)
      INTEGER ORDER_TX
      DOUBLE PRECISION MU
      
      MU = -1 + 2*(X-XMIN)/(XMAX-XMIN)
c     remember what's being expanded is beta(x), not T(x)
      GET_TX = 1/CHEBY(MU,ORDER_TX,COEF_TX)

      RETURN
      END
      
      
      FUNCTION cheby (x,order,coef)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)      
      PARAMETER (MAX_NP=50000,MAX_ORDER=30,MESH=500,NDIM=3)
      DIMENSION T(MAX_ORDER+1),COEF(MAX_ORDER+1)
      INTEGER ORDER
      
      T(1) = 1
      T(2) = x
      cheby = 0.
      do L=1,order+1
      if (L.gt.2) T(L)=2*x*T(L-1)-T(L-2)
      cheby = cheby + coef(L)*T(L)
      enddo
      RETURN
      END

      
      FUNCTION randsafe()
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      save iseed 
      data iseed /1/
      
 54   randsafe=ran1(iseed)
      if(randsafe.le.0) goto 54
      if(randsafe.ge.1) goto 54
      RETURN
      END

      SUBROUTINE mnbrak(ax,bx,cx,fa,fb,fc,func)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION ax,bx,cx,fa,fb,fc,func,GOLD,GLIMIT,TINY
      EXTERNAL func
      PARAMETER (GOLD=1.618034, GLIMIT=10.d0, TINY=1.d-20)
      DOUBLE PRECISION dum,fu,q,r,u,ulim
      fa=func(ax)
      fb=func(bx)
      if(fb.gt.fa)then
        dum=ax
        ax=bx
        bx=dum
        dum=fb
        fb=fa
        fa=dum
      endif
      cx=bx+GOLD*(bx-ax)
      fc=func(cx)
1     if(fb.ge.fc)then
        r=(bx-ax)*(fb-fc)
        q=(bx-cx)*(fb-fa)
        u=bx-((bx-cx)*q-(bx-ax)*r)/(2.d0*sign(max(dabs(q-r),TINY),q-r))
        ulim=bx+GLIMIT*(cx-bx)
        if((bx-u)*(u-cx).gt.0.d0)then
          fu=func(u)
          if(fu.lt.fc)then
            ax=bx
            fa=fb
            bx=u
            fb=fu
            return
          else if(fu.gt.fb)then
            cx=u
            fc=fu
            return
          endif
          u=cx+GOLD*(cx-bx)
          fu=func(u)
        else if((cx-u)*(u-ulim).gt.0.d0)then
          fu=func(u)
          if(fu.lt.fc)then
            bx=cx
            cx=u
            u=cx+GOLD*(cx-bx)
            fb=fc
            fc=fu
            fu=func(u)
          endif
        else if((u-ulim)*(ulim-cx).ge.0.d0)then
          u=ulim
          fu=func(u)
        else
          u=cx+GOLD*(cx-bx)
          fu=func(u)
        endif
        ax=bx
        bx=cx
        cx=u
        fa=fb
        fb=fc
        fc=fu
        goto 1
      endif
      return
      END

      DOUBLE PRECISION FUNCTION brent(ax,bx,cx,f,tol,xmin)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER ITMAX
      DOUBLE PRECISION ax,bx,cx,tol,xmin,f,CGOLD,ZEPS
      EXTERNAL f
      PARAMETER (ITMAX=1000,CGOLD=.3819660,ZEPS=1.0d-10)
      INTEGER iterr
      DOUBLE PRECISION a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,
     a                 tol1,tol2,u,v,w,x,xm
      a=min(ax,cx)
      b=max(ax,cx)
      v=bx
      w=v
      x=v
      e=0.d0
      fx=f(x)
      fv=fx
      fw=fx
      do 11 iterr=1,ITMAX
        xm=0.5*(a+b)
        tol1=tol*dabs(x)+ZEPS
        tol2=2.d0*tol1
        if(dabs(x-xm).le.(tol2-.5*(b-a))) goto 3
        if(dabs(e).gt.tol1) then
          r=(x-w)*(fx-fv)
          q=(x-v)*(fx-fw)
          p=(x-v)*q-(x-w)*r
          q=2.*(q-r)
          if(q.gt.0.) p=-p
          q=dabs(q)
          etemp=e
          e=d
         if(dabs(p).ge.dabs(.5*q*etemp).or.p.le.q*(a-x).or.p.ge.q*(b-x)) 
     *goto 1
          d=p/q
          u=x+d
          if(u-a.lt.tol2 .or. b-u.lt.tol2) d=sign(tol1,xm-x)
          goto 2
        endif
1       if(x.ge.xm) then
          e=a-x
        else
          e=b-x
        endif
        d=CGOLD*e
2       if(dabs(d).ge.tol1) then
          u=x+d
        else
          u=x+sign(tol1,d)
        endif
        fu=f(u)
        if(fu.le.fx) then
          if(u.ge.x) then
            a=x
          else
            b=x
          endif
          v=w
          fv=fw
          w=x
          fw=fx
          x=u
          fx=fu
        else
          if(u.lt.x) then
            a=u
          else
            b=u
          endif
          if(fu.le.fw .or. w.eq.x) then
            v=w
            fv=fw
            w=u
            fw=fu
          else if(fu.le.fv .or. v.eq.x .or. v.eq.w) then
            v=u
            fv=fu
          endif
        endif
11    continue
      pause 'brent exceed maximum iterations'
3     xmin=x
      brent=fx
      return
      END

      DOUBLE PRECISION FUNCTION f1dim(x)
      PARAMETER (MAX_NP=50000,MAX_ORDER=30,NDIM=3,MOP=MAX_ORDER+1)  
      DOUBLE PRECISION VALUE,x
CU    USES VALUE
      INTEGER j,ncom
      DOUBLE PRECISION pcom(MOP),xicom(MOP),xt(MOP)
      COMMON /f1com/ pcom,xicom,ncom
      do 11 j=1,ncom
        xt(j)=pcom(j)+x*xicom(j)
11    continue
      f1dim=VALUE(xt)
      return
      END

      
      SUBROUTINE linmin(p,xi,n,fret)
      PARAMETER (MAX_NP=50000,MAX_ORDER=30,NDIM=3,MOP=MAX_ORDER+1)  
      INTEGER n
      DOUBLE PRECISION fret,p(MOP),xi(MOP),TOL
      PARAMETER (TOL=1.d-2)
CU    USES brent,f1dim,mnbrak
      INTEGER j,ncom
      DOUBLE PRECISION ax,bx,fa,fb,fx,xmin,xx,pcom(MOP),xicom(MOP),brent
      COMMON /f1com/ pcom,xicom,ncom
      EXTERNAL f1dim
      ncom=n
      do 11 j=1,n
        pcom(j)=p(j)
        xicom(j)=xi(j)
11    continue
      ax=0.
      xx=1.
      call mnbrak(ax,xx,bx,fa,fx,fb,f1dim)
      fret=brent(ax,xx,bx,f1dim,TOL,xmin)
      do 12 j=1,n
        xi(j)=xmin*xi(j)
        p(j)=p(j)+xi(j)
12    continue
      return
      END

      
      SUBROUTINE frprmn(p,n,ftol,iterr,fret)
      PARAMETER (MAX_NP=50000,MAX_ORDER=30,NDIM=3,MOP=MAX_ORDER+1)  
      INTEGER iterr,n,ITMAX
      DOUBLE PRECISION fret,ftol,p(MOP),EPS,VALUE
      EXTERNAL VALUE
      PARAMETER (ITMAX=100000,EPS=1.d-7)
CU    USES FORCE,VALUE,linmin
      INTEGER its,j
      DOUBLE PRECISION dgg,fp,gam,gg,g(MOP),h(MOP),xi(MOP)

      fp=VALUE(p)
      call FORCE(p,xi)
      do 11 j=1,n
        g(j)=-xi(j)
        h(j)=g(j)
        xi(j)=h(j)
11    continue
      do 14 its=1,ITMAX
        iterr=its
        call linmin(p,xi,n,fret)
        if(2.*dabs(fret-fp).le.ftol*(dabs(fret)+dabs(fp)+EPS)) return
        fp=VALUE(p)
        call FORCE(p,xi)
        gg=0.
        dgg=0.
        do 12 j=1,n
          gg=gg+g(j)**2
C         dgg=dgg+xi(j)**2
          dgg=dgg+(xi(j)+g(j))*xi(j)
12      continue
        if(gg.eq.0.)return
        gam=dgg/gg
        do 13 j=1,n
          g(j)=-xi(j)
          h(j)=g(j)+gam*h(j)
          xi(j)=h(j)
13      continue
14    continue
      pause 'frprmn maximum iterations exceeded'
      return
      END

      
      FUNCTION ran1(idum)
      INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
      DOUBLE PRECISION ran1,AM,EPS,RNMX
      PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
     *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      INTEGER j,k,iv(NTAB),iy
      SAVE iv,iy
      DATA iv /NTAB*0/, iy /0/
      if (idum.le.0.or.iy.eq.0) then
        idum=max(-idum,1)
        do 11 j=NTAB+8,1,-1
          k=idum/IQ
           idum=IA*(idum-k*IQ)-IR*k
          if (idum.lt.0) idum=idum+IM
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum=idum+IM
      j=1+iy/NDIV
      iy=iv(j)
      iv(j)=idum

      ran1=min(AM*iy,RNMX)
      return
      END