在日常工作和科研中我们经常需要知道台站相对某个地震的震中距、方位角和反方位角。

目前比较流行的 distaz 程序采用了椭球模型进行计算,精度还算可以。

以下给出 distaz 的 C 版本和 Matlab 版本。同时附上编译好的 distaz.exe

distaz 程序的 C 版本

源码地址: http://www.seis.sc.edu/software/distaz/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/*
c
c   compile on Linux use :
c   $ cc distaz.c -o distaz -lm
c   compile on Windows use :
c   > gcc distaz.c -o distaz.exe -lm
c
*/

int main(int argc, char *argv[])
{
   double stalat, stalon, evtlat, evtlon;
   double delta, az, baz;
   double scolat, slon, ecolat, elon;
   double a,b,c,d,e,aa,bb,cc,dd,ee,g,gg,h,hh,k,kk;
   double rhs1,rhs2,sph,rad,del,daz,dbaz,pi,piby2;

   if (argc != 5) {
      printf("Usage: distaz sta_lat sta_lon evt_lat evt_lon\n");
      printf("       Returns:  Delta Baz Az\n");
      exit(1);
   } else {
      stalat = atof(argv[1]);
      stalon = atof(argv[2]);
      evtlat = atof(argv[3]);
      evtlon = atof(argv[4]);
   }
   if ((stalat == evtlat)&&(stalon == evtlon)) {
      printf("%6.3f %6.3f %6.3f\n", 0.0, 0.0, 0.0);
      exit(0);
   }

   pi=3.141592654;
   piby2=pi/2.0;
   rad=2.*pi/360.0;

   sph=1.0/298.257;

   scolat=piby2 - atan((1.-sph)*(1.-sph)*tan(stalat*rad));
   ecolat=piby2 - atan((1.-sph)*(1.-sph)*tan(evtlat*rad));
   slon=stalon*rad;
   elon=evtlon*rad;

   a=sin(scolat)*cos(slon);
   b=sin(scolat)*sin(slon);
   c=cos(scolat);
   d=sin(slon);
   e=-cos(slon);
   g=-c*e;
   h=c*d;
   k=-sin(scolat);

   aa=sin(ecolat)*cos(elon);
   bb=sin(ecolat)*sin(elon);
   cc=cos(ecolat);
   dd=sin(elon);
   ee=-cos(elon);
   gg=-cc*ee;
   hh=cc*dd;
   kk=-sin(ecolat);

   del=acos(a*aa + b*bb + c*cc);
   delta=del/rad;

   rhs1=(aa-d)*(aa-d)+(bb-e)*(bb-e)+cc*cc - 2.;
   rhs2=(aa-g)*(aa-g)+(bb-h)*(bb-h)+(cc-k)*(cc-k) - 2.;
   dbaz=atan2(rhs1,rhs2);
   if (dbaz<0.0) {
      dbaz=dbaz+2*pi;
   }
   baz=dbaz/rad;

   rhs1=(a-dd)*(a-dd)+(b-ee)*(b-ee)+c*c - 2.;
   rhs2=(a-gg)*(a-gg)+(b-hh)*(b-hh)+(c-kk)*(c-kk) - 2.;
   daz=atan2(rhs1,rhs2);
   if(daz<0.0) {
      daz=daz+2*pi;
   }
   az=daz/rad;
/*
c
c   Make sure 0.0 is always 0.0, not 360.
c
*/
   if(abs(baz-360.) < .00001) baz=0.0;
   if(abs(az-360.) < .00001) az=0.0;

   printf("%6.3f %6.3f %6.3f\n", delta, baz, az);
   exit(0);
}

distaz 程序的 Matlab 版本

源码地址: https://github.com/geodynamics/specfem3d_globe/blob/master/utils/Visualization/VTK_ParaView/matlab/distaz.m

function [dk,dd,daze,dazs]=distaz (sta,sto,epa,epo)
        rad=pi/180.0d0;

        sa  = atan(.993270*tan(sta*rad));
        ea  = atan(.993270*tan(epa*rad));
        ssa = sin(sa);
        csa = cos(sa);
        so  = sto*rad;
        eo  = epo*rad;
        sea = sin(ea);
        cea = cos(ea);
        ces = cos(eo-so);
        ses = sin(eo-so);
        
        if  (sa==ea)
           if (sto==epo)
              	dk =0.00;
					dd =0.00;
	  				daze=0.0;
               dazs=0.0;
                 return
  end
end

if sta==90.
   if epa==90.0
     dk =0.00;
	  dd =0.00;
	  daze=0.00;
     dazs=0.00;
     return
  end
end

      
      
if sta==-90.0
   if epa==-90.0
	  dk =0.00;
	  dd =0.00;
	  daze=0.00;
	  dazs=0.00;
	  return
	end
end
        dd = ssa*sea+csa*cea*ces;
        if dd~=0. , dd=atan(sqrt(1.0-dd*dd)/dd); end
        if dd==0. , dd=pi/2.0; end
        if dd<0.0 , dd=dd+pi; end
        dd = dd/rad;
        dk = dd*111.19;

        dazs = atan2(-ses,(ssa/csa*cea-sea*ces));
        daze = atan2(ses,(sea*csa/cea-ssa*ces));
        dazs = dazs/rad;
        daze = daze/rad;
        if (dazs<0.00), dazs=dazs+360.0; end
        if (daze<0.00), daze=daze+360.0; end
return

参考文献