      subroutine fitfld(itype,avetim,epoch,mono,nmax,nmaxt,nmaxtt 
     +     ,extfld)
c----------------------------------------------------------------
c contains the entry point to compon
c revised 5/31/1999 to neaten how leaving in Schmidt normalization
c jcain
c----------------------------------------------------------------

      integer extfld  
      common /datar/iskip,flatt,elong,alt,time,decl,declwt,dip,dipwt,hor
     +     ,horwt,x,xwt,y,ywt,z,zwt,b,bwt,xd,xdw,yd,ydw,zd,zdw,tfirst,tlast
c should find way to make array sizes variable maybe  
      common /coefs/g(61,61),h(61,61),gt(61,61),ht(61,61),gtt(61,61),
     +     htt(61,61),e(3),azero,maxn
      equivalence (e(1),e1),(e(2),e2),(e(3),e3) 
      dimension dxdh(61,61),dydh(61,61),dzdh(61,61),dxdg(61,61),dydg(61 
     +     ,61),dzdg(61,61),s(61,61)
      equivalence (dzdg(1,1),p(1,1)),(dzdh(1,1),dp(1,1))
       real*8 f
      common /fld/cd,ci,ch,cf,cz,cx,cy,cddeg,cideg,wt,tm, 
     +     cxd,cyd,czd,nor,nop,lat,lon,no,nonot,f(1) 
      dimension cp(61),sp(61),p(61,61),dp(61,61),const(61,61) 
      data a,flat/6378.16,298.25/, if/0./ 
      data p(1,1),cp(1),dp(1,2),dp(1,1),sp(1)/2*1.,3*0./
      save if,pi,pi2,rad,drad,a2,a4,b2,a2b2,a4b4,const,s

      if (if.eq.0) then
         if=1  
         pi=2.*asin(1.)  
         pi2=pi*2. 
         rad=180./pi 
         drad=pi/180.
         flat=1.-1./flat 
         a2=a**2 
         a4=a**4 
         b2=(a*flat)**2  
         a2b2=a2*(1.-flat**2)
         a4b4=a4*(1.-flat**4)

c-------------------------------------------------------------------
c     compute constants required for generating legendre polynomials
c------------------------------------------------------------------
c calculate factors here and review later what if mono=1 jc 5/27/99
         s(1,1)=-1.
         do n=2,maxn
           s(n,1) = s(n-1,1)*float(2*n-3)/float(n-1)
           s(1,n) = 0.
           j=2
           do m=2,n
             s(n,m) = s(n,m-1)*sqrt(float((n-m+1)*j)/float(n+m-2))
             s(m-1,n) = s(n,m)
             j=1
           enddo
         enddo
c don't see a neat way to combine do loops, later refinement
         do 1 n=2,maxn 
           do 1 m=1,n
 1          const(n,m)=float((n-2)**2-(m-1)**2)/float((2*n-3)*(2*n-5))

c-------------------------------------------------------------------
c     compute geocentric coordinates
c-------------------------------------------------------------------

      end if
 2    sinla=sin(flatt*drad) 
      sinla2=sinla**2 
      den2=a2-a2b2*sinla2 
      cosla2=1.-sinla2
      den=sqrt(den2)  
      fac=(((alt*den)+a2)/((alt*den)+b2))**2
      r=sqrt(alt*(alt+2.*den)+(a4-a4b4*sinla2)/den2)
      ct=sinla/sqrt(fac*cosla2+sinla2)
      aor=azero/r 
      st=sqrt(1.0-ct*ct)

c------------------------------------------------------------------
c     compute sine and cosine of difference between 
c     geodetic and geocentric latitudinal coordinates 
c------------------------------------------------------------------

      sind=sinla*st-sqrt(cosla2)*ct 
      cosd=sqrt(1.0-sind*sind)  

c------------------------------------------------------------------
c     longitude index 
c------------------------------------------------------------------

      if (elong.gt.180.) elong=elong-360.  
      lon=amin1(amax1(elong/10.0+19.0,1.0),36.0)

c------------------------------------------------------------------
c     latitude index  
c------------------------------------------------------------------

      lat=amin1(amax1(flatt*0.1+10.0,1.0),18.0) 
      rlong=elong*drad
      sp(2)=sin(rlong)
      cp(2)=cos(rlong)
      do 3 m=3,maxn 
c       sin(m-1)phi, eq.(5) 
        sp(m)=sp(2)*cp(m-1)+cp(2)*sp(m-1) 
c       cos(m-1)phi, eq.(6) 
        cp(m)=cp(2)*cp(m-1)-sp(2)*sp(m-1) 
 3    continue
c     generate associated legendre polynomials, eqs.(7),(8) 
      p(1,1)=1. 
      p(2,1)= ct
      dp(2,2)=ct
      p(2,2)=st 
      dp(2,1)=-st 
      do 5 n=3,maxn 
        do 5 m=1,n
          if (n.eq.m) then
             p(n,n)=st*p(n-1,n-1)  
             dp(n,n)=st*dp(n-1,n-1)+ct*p(n-1,n-1)
          else
 4           p(n,m)=ct*p(n-1,m)-const(n,m)*p(n-2,m)
             dp(n,m)=ct*dp(n-1,m)-st*p(n-1,m)-const(n,m)*dp(n-2,m) 
          end if
 5      continue

c-------------------------------------------------------------------
c     initialize to compute x,y,z 
c-------------------------------------------------------------------

      cx=0. 
      cy=0. 
      cz=0. 
      cxd=0.
      cyd=0.
      czd=0.
      ar=aor**mono
      tm=time-epoch 

c--------------------------------------------------------------------
c     compute x,y,z using best available  
c     parameters, eqs.(9),(10),(11) 
c--------------------------------------------------------------------

      do 6 n=mono,maxn
        fn=n
        ar=ar*aor 
        do 6 m=1,n
          fm=m-1
c because using Schmidt normalized, need s factor here
          p(n,m)=p(n,m)*ar*s(n,m)
          dp(n,m)=dp(n,m)*ar*s(n,m)
          temp=fn*p(n,m)*sind-dp(n,m)*cosd
          dxdg(n,m)=temp*cp(m)  
          dxdh(n,m)=temp*sp(m)  
          temp=fm*p(n,m)/st 
          dydg(n,m)=-temp*sp(m) 
          dydh(n,m)=temp*cp(m)  
          temp=fn*p(n,m)*cosd+dp(n,m)*sind
          dzdg(n,m)=temp*cp(m)  
          dzdh(n,m)=temp*sp(m)  
c-------------------------------------------------------------------
c         add time terms, eqs.(15) and (16) 
c-------------------------------------------------------------------
          gnm=(tm*gtt(n,m)+gt(n,m))*tm+g(n,m) 
          hnm=(tm*htt(n,m)+ht(n,m))*tm+h(n,m) 
          cx=cx+gnm*dxdg(n,m)+hnm*dxdh(n,m) 
          cy=cy+gnm*dydg(n,m)+hnm*dydh(n,m) 
          cz=cz+gnm*dzdg(n,m)+hnm*dzdh(n,m) 
          gnmd=gt(n,m)+2.*tm*gtt(n,m) 
          hnmd=ht(n,m)+2.*tm*htt(n,m) 
          cxd=cxd+gnmd*dxdg(n,m)+hnmd*dxdh(n,m) 
          cyd=cyd+gnmd*dydg(n,m)+hnmd*dydh(n,m) 
          czd=czd+gnmd*dzdg(n,m)+hnmd*dzdh(n,m) 
 6     continue
      if (extfld.ne.0) then  
         t1=e2*cp(2)+e3*sp(2)
         t2=e1*st-t1*ct  
         t1=e1*ct+t1*st  
         cx=cx-t2*cosd+t1*sind 
         cy=cy+e2*sp(2)-e3*cp(2) 
         cz=cz+t2*sind+t1*cosd 
      end if
c--------------------------------------------------------------------
c     compute horizontal, total field, dip, and 
c     declination, eqs.(17),(18),(19),(20)
c--------------------------------------------------------------------
      ch=sqrt(cx**2+cy**2)
      cf=sqrt(ch**2+cz**2)
      ci=asin(cz/cf)  
      cd=atan2(cy,cx) 
      cideg=ci*rad
      cddeg=cd*rad
      return

c--------------------------------------------------------------------
c     begin entry for compon
c     compute coefficients of g(n,m) and  
c     h(n,m) for the next approximation
c
c--------------------------------------------------------------------
 
      entry compon (itype,avetim,epoch,mono,nmax,nmaxt,nmaxtt,extfld) 

      i=1 
      go to (8,12,14,20,22,18,16,220,230,240),itype 
c     compute coefficients when declination is given
    8 fi=decl*drad-cd 
      if(fi.lt.-pi)fi=fi+pi2
      if(fi.gt.pi) fi=fi-pi2
      fi=fi*ch
      wt=1./declwt
      if(nmax.eq.0) go to 28
      t1=cx/ch
      t2=cy/ch
      do 9 n=mono,nmax
        do 9 m=1,n
c         ch*ddecl/dg(n,m)
          f(i)=(t1*dydg(n,m)-t2*dxdg(n,m))
          if(m.eq.1) go to 9
          i=i+1 
c         ch*ddecl/dh(n,m)
          f(i)=(t1*dydh(n,m)-t2*dxdh(n,m))
    9 i=i+1 
      go to 24
c     compute coefficients when dip is given
c     dip=arctan(z/h) 
   12 fi=(dip*drad-ci)*cf 
      wt=1./dipwt 
      if(nmax.eq.0)go to 28 
      t1=ch/cf
      temp=cz/(ch*cf) 
      t2=cx*temp
      t3=cy*temp
      do 13 n=mono,nmax 
        do 13 m=1,n 
c         f(i)=cf*ddip/dg(n,m)  
          f(i)=(t1*dzdg(n,m)-t2*dxdg(n,m)-t3*dydg(n,m)) 
          if(m.eq.1) go to 13 
          i=i+1 
c         f(i)=cf*ddip/dh(n,m)  
          f(i)=(t1*dzdh(n,m)-t2*dxdh(n,m)-t3*dydh(n,m)) 
   13 i=i+1 
      go to 24
c     compute coefficients when the horizontal
c     component (hor) is given, eqs,(25),(26) 
c     hor=sqrt(x*x+y*y) 
c     t1=dhor/dx
c     t2=dhor/dy
   14 fi=hor-ch 
      wt=1./horwt 
      if(nmax.eq.0) go to 28
      t1=cx/ch
      t2=cy/ch
      do 15 n=mono,nmax 
        do 15 m=1,n 
c         f(i)=dhor/dg(n,m) 
          f(i)=(t1*dxdg(n,m)+t2*dydg(n,m))
          if(m.eq.1) go to 15 
          i=i+1 
c         f(i)=dhor/dh(n,m) 
          f(i)=(t1*dxdh(n,m)+t2*dydh(n,m))
   15 i=i+1 
      go to 24
c     compute coefficients when total 
c     field b is given, eqs.(27),(28) 
c     b=sqrt(x*x+y*y+z*z) 
c     t1=db/dx
c     t2=db/dy
c     t3=db/dz
   16 fi=b-cf 
      wt=1./bwt 
      if(nmax.eq.0) go to 28
      t1=cx/cf
      t2=cy/cf
      t3=cz/cf
      do 17 n=mono,nmax 
        do 17 m=1,n 
c         f(i)=db/dg(n,m) 
          f(i)=(t1*dxdg(n,m)+t2*dydg(n,m)+t3*dzdg(n,m)) 
          if(m.eq.1) go to 17 
          i=i+1 
c         f(i)=db/dh(n,m) 
          f(i)=(t1*dxdh(n,m)+t2*dydh(n,m)+t3*dzdh(n,m)) 
   17 i=i+1 
      go to 24
c     compute coefficients when the z component is given
   18 fi=z-cz 
      wt=1./zwt 
      if(nmax.eq.0) go to 28
      do 19 n=mono,nmax 
        do 19 m=1,n 
          f(i)=dzdg(n,m)
          if (m.eq.1) goto19
          i=i+1 
          f(i)=dzdh(n,m)
   19 i=i+1 
      go to 24
c     compute coefficients when the x component is given
   20 fi=x-cx 
      wt=1./xwt 
      if(nmax.eq.0) go to 28
      do 21 n=mono,nmax 
        do 21 m=1,n 
          f(i)=dxdg(n,m)
          if (m.eq.1) goto21
          i=i+1 
          f(i)=dxdh(n,m)
   21 i=i+1 
      go to 24
c     compute coefficients when the y component is given
   22 fi=y-cy 
      wt=1./ywt 
      if(nmax.eq.0) go to 28
      do 23 n=mono,nmax 
        do 23 m=1,n 
          f(i)=dydg(n,m)
          if (m.eq.1) goto23
          i=i+1 
          f(i)=dydh(n,m)
   23 i=i+1 
      go to 24
c     compute partials when xdot is given 
  220 tfact=time-avetim 
      wt=1./xdw 
      if(no.eq.0) go to 222 
      do 221 n=1,no 
  221 f(n)=0. 
  222 i=no+1
      do 223 n=mono,nmaxt 
      do 223 m=1,n
      f(i)=dxdg(n,m)  
      if(m.eq.1) go to 223
      i=i+1 
      f(i)=dxdh(n,m)  
  223 i=i+1 
      fi=xd-cxd 
      if(nmaxtt.eq.0) go to 28  
      ttf=2.*tfact
      do 224 n=mono,nmaxtt
      do 224 m=1,n
      f(i)=ttf*dxdg(n,m)
      if(m.eq.1) go to 224
      i=i+1 
      f(i)=ttf*dxdh(n,m)
  224 i=i+1 
      go to 28
c     compute partials when ydot is given 
  230 tfact=time-avetim 
      wt=1./ydw 
      if(no.eq.0) go to 232 
      do 231 n=1,no 
  231 f(n)=0. 
  232 i=no+1
      do 233 n=mono,nmaxt 
      do 233 m=1,n
      f(i)=dydg(n,m)  
      if(m.eq.1) go to 233
      i=i+1 
      f(i)=dydh(n,m)  
  233 i=i+1 
      fi=yd-cyd 
      if(nmaxtt.eq.0) go to 28  
      ttf=2.*tfact
      do 234 n=mono,nmaxtt
      do 234 m=1,n
      f(i)=ttf*dydg(n,m)
      if(m.eq.1) go to 234
      i=i+1 
      f(i)=ttf*dydh(n,m)
  234 i=i+1 
      go to 28
c     compute partials if zdot is given 
  240 tfact=time-avetim 
      wt=1./zdw 
      if(no.eq.0) go to 242 
      do 241 n=1,no 
  241 f(n)=0. 
  242 i=no+1
      do 243 n=mono,nmaxt 
      do 243 m=1,n
      f(i)=dzdg(n,m)  
      if(m.eq.1) go to 243
      i=i+1 
      f(i)=dzdh(n,m)  
  243 i=i+1 
      fi=zd-czd 
      if(nmaxtt.eq.0) go to 28  
      ttf=2.*tfact
      do 244 n=mono,nmaxtt
      do 244 m=1,n
      f(i)=ttf*dzdg(n,m)
      if(m.eq.1) go to 244
      i=i+1 
      f(i)=ttf*dzdh(n,m)
  244 i=i+1 
      go to 28
c     add time**1 terms 
   24 tfact=time-avetim 
      if(nmaxt.eq.0) go to 26 
      j=1 
      do 25 n=mono,nmaxt
        do 25 m=1,n 
          f(i)=f(j)*tfact 
          i=i+1 
          j=j+1 
          if (m.eq.1) goto25
          f(i)=f(j)*tfact 
          i=i+1 
          j=j+1 
   25 continue
c     add time**2 terms 
   26 if(nmaxtt.eq.0) go to 28  
      j=1 
      tfact=tfact*tfact 
      do 27 n=mono,nmaxtt 
        do 27 m=1,n 
          f(i)=f(j)*tfact 
          i=i+1 
          j=j+1 
          if (m.eq.1) goto27
          f(i)=f(j)*tfact 
          i=i+1 
          j=j+1 
   27 continue
   28 if(extfld.eq.0) go to 36  
      if(itype.le.7) go to 280  
      f(i)=0. 
      f(i+1)=0. 
      f(i+2)=0. 
      i=i+3 
      go to 36
c     add external field terms  
  280 dxde1=ct*sind-st*cosd 
      dzde1=st*sind+ct*cosd 
      dxde2=cp(2)*dzde1 
      dxde3=sp(2)*dzde1 
      dzde2=-cp(2)*dxde1
      dzde3=-sp(2)*dxde1
      go to (29,30,31,34,35,33,32),itype  
c     coefficients when declination (d) is given
c     f(i)=dd/de1 
   29 f(i)=-t2*dxde1  
      i=i+1 
c     f(i)=dd/de2 
      f(i)=t1*sp(2)-t2*dxde2
      i=i+1 
c     f(i)=dd/de3 
      f(i)=-t1*cp(2)-t2*dxde3 
      i=i+1 
      go to 36
c     coefficients when dip (i) is given  
c     f(i)=di/de1 
   30 f(i)=t1*dzde1-t2*dxde1
      i=i+1 
c     f(i)=di/de2 
      f(i)=t1*dzde2-t2*dxde2-t3*sp(2) 
      i=i+1 
c     f(i)=di/de3 
      f(i)=t1*dzde3-t2*dxde3+t3*cp(2) 
      i=i+1 
      go to 36
c     coefficients when the horizontal
c     component (hor) is given  
c     f(i)=dhor/de1 
   31 f(i)=t1*dxde1 
      i=i+1 
c     f(i)=dhor/de2 
      f(i)=t1*dxde2+t2*sp(2)
      i=i+1 
c     f(i)=dhor/de3 
      f(i)=t1*dxde3-t2*cp(2)
      i=i+1 
      go to 36
c     coefficients when total field is given
c     f(i)=df/de1 
   32 f(i)=t1*dxde1+t3*dzde1
      i=i+1 
c     f(i)=df/de2 
      f(i)=t1*dxde2+t2*sp(2)+t3*dzde2 
      i=i+1 
c     f(i)=df/de3 
      f(i)=t1*dxde3-t2*cp(2)+t3*dzde3 
      i=i+1 
      go to 36
c     coefficients when the z component is given
   33 f(i)=dzde1
      i=i+1 
      f(i)=dzde2
      i=i+1 
      f(i)=dzde3
      i=i+1 
      go to 36
c     coefficients when the x component is given
   34 f(i)=dxde1
      i=i+1 
      f(i)=dxde2
      i=i+1 
      f(i)=dxde3
      i=i+1 
      go to 36
c     coefficients when the y component is given
   35 f(i)=0.0
      i=i+1 
      f(i)=sp(2)
      i=i+1 
      f(i)=-cp(2) 
      i=i+1 
c     add observation term
   36 f(i)=fi 
      return
      end 
