sub delaz(plon1 as float, plat1 as float, plon2 as float, plat2 as float, delta as float, azimuth as float)
' get distance and azimuth from  (plat1, plon1) to (plat2, plon2)


dim f,g, tplat2,tplon2,glat,dca,dcb,dcc,gplat1,epdca,epdcb,epdcc as float
dim par1, t1,t2,t3,par2,t4,apar2,par3 as float

  f=3.1415926/180.0 '''1.745329e-02 
  g=1.0/f
' --- distaz
        tplat2 = plat2*f
        tplon2 = plon2*f
  glat=atn(.993277*tan(tplat2))
        dca = cos(glat)*cos(tplon2)
        dcb = cos(glat)*sin(tplon2)
        dcc = sin(glat)
        plat2 = tplat2/f
  gplat1 = atn(.993277*tan(plat1*f))
  epdca = cos(gplat1)*cos(plon1*f)
  epdcb = cos(gplat1)*sin(plon1*f)
  epdcc = sin(gplat1)

 '---- calculate distance and azimuth from earthquake to station

  par1=(1-.5*((dca -epdca)^2 + (dcb -epdcb)^2+(dcc  - epdcc)^2))
        DELTA  = g*atn(sqr((1/(par1*par1)) - 1))
        if (par1<0) then 
            DELTA  = 180 - DELTA
        elseif(delta=0) then 
            DELTA = .001
        end if

        t1 = dca  - sin(plon1*f)
        t2 = dcb  + cos(plon1*f)
        t3 = sin(delta *f)
        par2 = (0.5*(t1*t1 + t2*t2 + dcc *dcc ) - 1.0 )/t3
        t1 = dca  - sin(gplat1)*cos(plon1*f)
        t2 = dcb  - sin(gplat1)*sin(plon1*f)
        t4 = dcc  + cos(gplat1)
        apar2 = abs(par2)
        AZIMUTH  = 90.0

        if ( apar2 < 1.0) then 
         if par2=0.0 then 
            par2=.00000001
         end if
         if par2=1.0 then 
            par2=.99999999
         end if
         AZIMUTH  = g*atn(1/sqr((1/par2^2)-1))
        end if

        par3 = (.5*(t1*t1+t2*t2+t4*t4) - 1)/t3

      if (par2<0) then 

        if (par3<0) then 
           AZIMUTH  = 180 + AZIMUTH 
        else
           AZIMUTH  = 360 - AZIMUTH 
        end if

      else
        if(par3<0) then 
           AZIMUTH  = 180 - AZIMUTH
        end if
      end if

end sub