subroutine mkfact(tm,h,bfm,thb,bki,nhts)
c
c      create table of magnetic field components
c
      real h(nhts),bfm(nhts),thb(nhts),bki(nhts)
      data pio180/0.017453292/,pi/3.14159265/
C     rate of rotation produced by n=1 cm-3, b=1 gauss
      data rate/1.8978873e-6/
C     Jicamarca position
      data gdlat/-11.951/,elon/-76.8667/,gdalt/0./,re/6280.0/
C     3.0, 4.5, 6.0 and 3.38N degree positions 
c      data decd/-12.88/,ham/-4.617/
c      data decd/-13.1/,ham/-5.0/  ! made up
c     data decd/-14.56/,ham/-5.33/
c      data decd/-16.24/,ham/-6.03/
      data decd/-9.52/,ham/-3.20/
c=======
c
c     convert from geodetic to geocentric coordinates
c
      call convrt(1,gdlat,gdalt,gclat,re)
c
c     convert from degrees to radians
c
      gclat=gclat*pio180
      elong=elon*pio180
c
c     create components of vector from center of earth
c     to antenna center

      ca=cos(elong)
      cd=cos(gclat)
      sa=sin(elong)
      sd=sin(gclat)
      rx=re*ca*cd
      ry=re*sa*cd
      rz=re*sd
c
c     create antenna normal vector (unit vector in direction of k)
c
      dec=decd*pio180
      ha=ham*pi/720.+elong
      ch=cos(ha)
      sh=sin(ha)
      cdec=cos(dec)
      sdec=sin(dec)
      ax=ch*cdec
      ay=sh*cdec
      az=sdec
      do i=1,nhts
c
c       create vector from center of earth to scattering volume
c
        x=rx+ax*h(i)
        y=ry+ay*h(i)
        z=rz+az*h(i)
        x2y2=x*x+y*y
        r=sqrt(x2y2+z*z)
        x2y2=sqrt(x2y2)
c
c       convert to spherical coordinates
c
        st=x2y2/r
        ct=z/r
        sp=y/x2y2
        cp=x/x2y2
c       call milmag(tm,r,st,ct,sp,cp,br,bt,bp,b)
        theta=atan2(st,ct)
        phi=atan2(sp,cp)
        call geobfield(tm,r,theta,phi,br,bt,bp,b)

c
c       convert to cartesian coordinates
c
        bx=br*st*cp+bt*ct*cp-bp*sp
        by=br*st*sp+bt*ct*sp+bp*cp
        bz=br*ct-bt*st
        bk=(ax*bx+ay*by+az*bz)
        thb(i)=acos(abs(bk/b))*57.29577951
        bfm(i)=b
        bki(i)=1.0/bk

      end do
      return
      end

      subroutine gausfit(drift,hist,vz,m)
c
c     fit procedure for estimating velocities
c
      real drift(m),hist(m),vz
      real x(1000),y(1000)
      parameter(n=6)
      real a(n),tol,wa(m*n+5*n+m),wa2(m*n+5*n+m),fvec(m)
      integer info,iwa(n),lwa
      external gauss
      common/gf/x,y
      data vlast/0.0/

      lwa=m*n+5*n+m
      tol=1.0e-6
      info=0
      amax=0.0

      do i=1,m
         x(i)=drift(i)
         y(i)=hist(i)
         amax=max(amax,hist(i))
      end do

      a(1)=amax
      a(2)=-20.0
      a(3)=20.0

      a(4)=amax/4.0
      a(5)=20.0
      a(6)=20.0

c      write(*,*) m,n,a,tol,info

      call lmdif1(gauss,m,n,a,fvec,tol,info,iwa,wa,lwa)
     
      call fdjac2(gauss,m,n,a,fvec,wa,m,iflag,0.0e0,wa2)
      if(abs(a(1)).gt.abs(a(4))) then
         vz=a(2)
         i=2
      else
         vz=a(5)
         i=5
      end if

      err=0.0
      do j=1,m
         err=err+(wa(j+(i-1)*m))**2
      end do
      if(err.gt.0.0) err=sqrt(err**-1)
      
c
c     try fitting single distribution
c
      a(1)=amax
      a(2)=-20.0
      a(3)=20.0

      call lmdif1(gauss,m,n/2,a,fvec,tol,info,iwa,wa,lwa)
      call fdjac2(gauss,m,n/2,a,fvec,wa,m,iflag,0.0e0,wa2)
      vz1=a(2)
      i=2

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


      if(err1.lt.err.or.abs(vz-vz1).lt.10.0) vz=vz1
      vz=vz1
      vlast=vz

      return
      end

      subroutine gauss(m,n,a,fvec,iflag)
c
      integer m,n,iflac
      real fvec(m),a(n),b
      real x(1000),y(1000)
      common/gf/x,y
c
      
      do i=1,m
         sig=sqrt(max(1.0,y(i)))
         b=abs(a(1))*exp(-(a(2)-x(i))**2/(2.0*a(3)**2))
         if(n.gt.3) then
            c=abs(a(4))*exp(-(a(5)-x(i))**2/(2.0*a(6)**2))
         else
            c=0.0
         end if
         fvec(i)=(b+c-y(i))**2/sig**2
      end do
      return
      end