source: trunk/src/libsrc/util/geo_to_km.c @ 7513

Revision 7513, 17.9 KB checked in by baker, 6 months ago (diff)

cleanup warning: variable set but not used [-Wunused-but-set-variable]

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1#include <stdio.h>
2#include <math.h>
3
4#include "geo.h"
5
6/*
7Calculations are based upon the reference spheroid of 1968 and
8are defined by the major radius (RAD) and the flattening (FL).
9*/
10#define semi_major 6378.160
11#define semi_minor 6356.775   
12
13int geo_to_km(double lat1,double lon1,double lat2,double lon2,double* dist,double* azm) 
14{
15
16/*
17=====================================================================
18 PURPOSE:  To compute the distance and azimuth between locations.
19=====================================================================
20 INPUT ARGUMENTS:
21    lat1:     Event latitude in decimal degrees, North positive. [r]
22    lon1:     Event longitude, East positive. [r]
23    lat2:     station latitude. [r]
24    lon2:     station longitude. [r]
25=====================================================================
26 OUTPUT ARGUMENTS:
27    DIST:    epicentral distance in km. [r]
28    AZM:      azimuth in degrees. [r]
29=====================================================================
30*/
31/*
32Calculations are based upon the reference spheroid of 1968 and
33are defined by the major radius (RAD) and the flattening (FL).
34*/
35
36      const double a = semi_major;
37      const double b = semi_minor;   
38      double torad, todeg;
39      double aa, bb, cc, dd, top, bottom, lambda12, az, temp;
40      double v1, v2;
41      double fl, e2, eps, eps0;
42      double b0, x2, y2, z2, z1, u1p, u2p, xdist;
43      double lat1rad, lat2rad, lon1rad, lon2rad;
44      double coslon1, sinlon1, coslon2, sinlon2;
45      double coslat1, sinlat1, coslat2, sinlat2;
46      double tanlat1, tanlat2, cosazm, sinazm;
47
48      double c0, c2, c4, c6;
49
50      double c00=1.0, c01=0.25, c02=-0.046875, c03=0.01953125;
51      double c21=-0.125, c22=0.03125, c23=-0.014648438;
52      double c42=-0.00390625, c43=0.0029296875;
53      double c63=-0.0003255208;
54
55/*  Check for special conditions */
56     if( lat1 == lat2 && lon1 == lon2 ) {
57         *azm = 0.0;
58         *dist= 0.0;
59         return(1);
60     }
61/* - Initialize.             */
62
63      torad = PI / 180.0;
64      todeg = 1.0 / torad;
65      fl = ( a - b ) / a;
66      e2 = 2.0*fl - fl*fl;
67      eps = e2 / ( 1.0 - e2);
68/*
69* - Convert event location to radians.
70*   (Equations are unstable for latidudes of exactly 0 degrees.)
71*/
72      temp=lat1;
73      if(temp == 0.) temp=1.0e-08;
74      lat1rad=torad*temp;
75      lon1rad=torad*lon1;
76
77      temp=lat2;
78      if(temp == 0.) temp=1.0e-08;
79      lat2rad=torad*temp;
80      lon2rad=torad*lon2;
81
82/*
83      Compute some of the easier and often used terms.
84*/
85      coslon1 = cos(lon1rad);
86      sinlon1 = sin(lon1rad);
87      coslon2 = cos(lon2rad);
88      sinlon2 = sin(lon2rad);
89      tanlat1 = tan(lat1rad);
90      tanlat2 = tan(lat2rad);
91      sinlat1 = sin(lat1rad);
92      coslat1 = cos(lat1rad);
93      sinlat2 = sin(lat2rad);
94      coslat2 = cos(lat2rad);
95/*
96    The radii of curvature are compute from an equation defined in
97    GEODESY by Bomford, Appendix A (page 647).
98    v = semi_major/sqrt(1-e*e*sin(lat)*sin(lat))
99*/
100      v1 = a / sqrt( 1.0 - e2*sinlat1*sinlat1 );  /* radii of curvature */
101      v2 = a / sqrt( 1.0 - e2*sinlat2*sinlat2 );  /* radii of curvature */
102      aa = tanlat2 / ((1.0+eps)*tanlat1);
103      bb = e2*(v1*coslat1)/(v2*coslat2);
104      lambda12 = aa + bb;
105      top = sinlon2*coslon1 - coslon2*sinlon1;
106      bottom = lambda12*sinlat1-coslon2*coslon1*sinlat1-sinlon2*sinlon1*sinlat1;
107      az = atan2(top,bottom)*todeg;
108      if( az < 0.0 ) az = 360 + az;
109      *azm = az;
110      az = az * torad;
111      cosazm = cos(az);
112      sinazm = sin(az);
113
114/*
115   Now compute the distance using the equations on page 121 in GEODESY by
116   Bomford (2.15 Reverse formulae).  There is some numerical problem with
117   the following formulae.
118   If the station is in the southern hemisphere and the event is in the
119   northern, these equations give the longer, not the shorter distance between
120   the two locations.  Since the equations are messy, the simplist solution
121   is to reverse the order of the lat,lon pairs.  This means that the azimuth
122   must also be recomputed to get the correct distance.
123*/
124      if( lat2rad < 0.0 ) {
125          temp = lat1rad;
126          lat1rad = lat2rad;
127          lat2rad = temp;
128          temp = lon1rad;
129          lon1rad = lon2rad;
130          lon2rad = temp;
131
132          coslon1 = cos(lon1rad);
133          sinlon1 = sin(lon1rad);
134          coslon2 = cos(lon2rad);
135          sinlon2 = sin(lon2rad);
136          tanlat1 = tan(lat1rad);
137          tanlat2 = tan(lat2rad);
138          sinlat1 = sin(lat1rad);
139          coslat1 = cos(lat1rad);
140          sinlat2 = sin(lat2rad);
141          coslat2 = cos(lat2rad);
142
143          v1 = a / sqrt( 1.0 - e2*sinlat1*sinlat1 ); 
144          v2 = a / sqrt( 1.0 - e2*sinlat2*sinlat2 ); 
145
146          aa = tanlat2 / ((1.0+eps)*tanlat1);
147          bb = e2*(v1*coslat1)/(v2*coslat2);
148          lambda12 = aa + bb;
149
150          top = sinlon2*coslon1 - coslon2*sinlon1;
151          bottom =lambda12*sinlat1-coslon2*coslon1*sinlat1-
152                  sinlon2*sinlon1*sinlat1;
153          az = atan2(top,bottom);
154          cosazm = cos(az);
155          sinazm = sin(az);
156           
157       }
158
159       eps0 = eps * ( coslat1*coslat1*cosazm*cosazm + sinlat1*sinlat1 );
160       b0 = (v1/(1.0+eps0)) * sqrt(1.0+eps*coslat1*coslat1*cosazm*cosazm);
161     
162       x2 = v2*coslat2*(coslon2*coslon1+sinlon2*sinlon1);
163       y2 = v2*coslat2*(sinlon2*coslon1-coslon2*sinlon1);
164       z2 = v2*(1.0-e2)*sinlat2;
165       z1 = v1*(1.0-e2)*sinlat1;
166
167       c0 = c00 + c01*eps0 + c02*eps0*eps0 + c03*eps0*eps0*eps0;
168       c2 =       c21*eps0 + c22*eps0*eps0 + c23*eps0*eps0*eps0;
169       c4 =                  c42*eps0*eps0 + c43*eps0*eps0*eps0;
170       c6 =                                  c63*eps0*eps0*eps0;
171
172       bottom = cosazm*sqrt(1.0+eps0);
173       u1p = atan2(tanlat1,bottom);
174         
175       top = v1*sinlat1+(1.0+eps0)*(z2-z1);
176       bottom = (x2*cosazm-y2*sinlat1*sinazm)*sqrt(1.0+eps0);
177       u2p = atan2(top,bottom);
178
179       aa = c0*(u2p-u1p);
180       bb = c2*(sin(2.0*u2p)-sin(2.0*u1p));
181       cc = c4*(sin(4.0*u2p)-sin(4.0*u1p));
182       dd = c6*(sin(6.0*u2p)-sin(6.0*u1p));
183
184       xdist = fabs(b0*(aa+bb+cc+dd));
185       *dist = xdist;
186           return(1);
187}
188
189
190
191int     geo_to_km_deg (double lat1, double lon1, double lat2, double lon2,
192                                        double *dist, double *xdeg, double *azm)
193{
194
195/*
196=====================================================================
197 PURPOSE:  To compute the distance and azimuth between locations.
198=====================================================================
199 INPUT ARGUMENTS:
200    lat1:     Event latitude in decimal degrees, North positive. [r]
201    lon1:     Event longitude, East positive. [r]
202    lat2:     Array of station latitudes. [r]
203    lon2:     Array of station longitudes. [r]
204=====================================================================
205 OUTPUT ARGUMENTS:
206    DIST:    epicentral distance in km. [r]
207    XDEG:    epicentral distance in degrees. [r]
208    AXM:     azimuth in degrees. [r]
209=====================================================================
210*/
211/*
212Calculations are based upon the reference spheroid of 1968 and
213are defined by the major radius (RAD) and the flattening (FL).
214*/
215
216      const double a = semi_major;
217      const double b = semi_minor;
218      double torad, todeg;
219      double aa, bb, cc, dd, top, bottom, lambda12, az, temp;
220      double v1, v2;
221      double fl, e2, eps, eps0;
222      double b0, x2, y2, z2, z1, u1p, u2p, xdist;
223      double lat1rad, lat2rad, lon1rad, lon2rad;
224      double coslon1, sinlon1, coslon2, sinlon2;
225      double coslat1, sinlat1, coslat2, sinlat2;
226      double tanlat1, tanlat2, cosazm, sinazm;
227          double onemec2;
228          double aa2, bb2, cc2, dd2, ee2, ff2, aa3, bb3, cc3, dd3, ee3, ff3;
229          double aminus, bminus, cminus, aplus, bplus, cplus;
230          double thg, sd, sc, deg;
231
232      double c0, c2, c4, c6;
233
234      double c00=1.0, c01=0.25, c02=-0.046875, c03=0.01953125;
235      double c21=-0.125, c22=0.03125, c23=-0.014648438;
236      double c42=-0.00390625, c43=0.0029296875;
237      double c63=-0.0003255208;
238
239/*  Check for special conditions                                  */
240
241     if( lat1 == lat2 && lon1 == lon2 ) {
242         *azm = 0.0;
243         *dist= 0.0;
244         return(0);
245     }
246/* - Initialize.             */
247
248      torad = PI / 180.0;
249      todeg = 1.0 / torad;
250      fl = ( a - b ) / a;
251      e2 = 2.0*fl - fl*fl;
252      eps = e2 / ( 1.0 - e2);
253      onemec2 = 1.0 - e2;
254/*
255* - Convert event location to radians.
256*   (Equations are unstable for latidudes of exactly 0 degrees.)
257*/
258
259      temp=lat1;
260      if(temp == 0.) temp=1.0e-08;
261      lat1rad=torad*temp;
262      lon1rad=torad*lon1;
263
264      temp=lat2;
265      if(temp == 0.) temp=1.0e-08;
266      lat2rad=torad*temp;
267      lon2rad=torad*lon2;
268
269/*
270      Compute some of the easier and often used terms.
271*/
272      coslon1 = cos(lon1rad);
273      sinlon1 = sin(lon1rad);
274      coslon2 = cos(lon2rad);
275      sinlon2 = sin(lon2rad);
276      tanlat1 = tan(lat1rad);
277      tanlat2 = tan(lat2rad);
278      sinlat1 = sin(lat1rad);
279      coslat1 = cos(lat1rad);
280      sinlat2 = sin(lat2rad);
281      coslat2 = cos(lat2rad);
282
283/*
284    The radii of curvature are compute from an equation defined in
285    GEODESY by Bomford, Appendix A (page 647).
286    v = semi_major/sqrt(1-e*e*sin(lat)*sin(lat))
287*/
288      v1 = a / sqrt( 1.0 - e2*sinlat1*sinlat1 );  /* radii of curvature */
289      v2 = a / sqrt( 1.0 - e2*sinlat2*sinlat2 );  /* radii of curvature */
290
291      aa = tanlat2 / ((1.0+eps)*tanlat1);
292      bb = e2*(v1*coslat1)/(v2*coslat2);
293      lambda12 = aa + bb;
294
295      top = sinlon2*coslon1 - coslon2*sinlon1;
296      bottom = lambda12*sinlat1-coslon2*coslon1*sinlat1-sinlon2*sinlon1*sinlat1;
297      az = atan2(top,bottom)*todeg;
298      if( az < 0.0 ) az = 360 + az;
299      *azm = az;
300      az = az * torad;
301      cosazm = cos(az);
302      sinazm = sin(az);
303
304/*
305   Now compute the distance using the equations on page 121 in GEODESY by
306   Bomford (2.15 Reverse formulae).  There is some numerical problem with
307   the following formulae.
308   If the station is in the southern hemisphere and the event is in the
309   northern, these equations give the longer, not the shorter distance between
310   the two locations.  Since the equations are messy, the simplist solution
311   is to reverse the order of the lat,lon pairs.  This means that the azimuth
312   must also be recomputed to get the correct distance.
313*/
314
315      if( lat2rad < 0.0 ) {
316          temp = lat1rad;
317          lat1rad = lat2rad;
318          lat2rad = temp;
319          temp = lon1rad;
320          lon1rad = lon2rad;
321          lon2rad = temp;
322
323          coslon1 = cos(lon1rad);
324          sinlon1 = sin(lon1rad);
325          coslon2 = cos(lon2rad);
326          sinlon2 = sin(lon2rad);
327          tanlat1 = tan(lat1rad);
328          tanlat2 = tan(lat2rad);
329          sinlat1 = sin(lat1rad);
330          coslat1 = cos(lat1rad);
331          sinlat2 = sin(lat2rad);
332          coslat2 = cos(lat2rad);
333
334          v1 = a / sqrt( 1.0 - e2*sinlat1*sinlat1 );
335          v2 = a / sqrt( 1.0 - e2*sinlat2*sinlat2 );
336
337          aa = tanlat2 / ((1.0+eps)*tanlat1);
338          bb = e2*(v1*coslat1)/(v2*coslat2);
339          lambda12 = aa + bb;
340
341          top = sinlon2*coslon1 - coslon2*sinlon1;
342          bottom =lambda12*sinlat1-coslon2*coslon1*sinlat1-
343                  sinlon2*sinlon1*sinlat1;
344          az = atan2(top,bottom);
345          cosazm = cos(az);
346          sinazm = sin(az);
347
348       }
349
350       eps0 = eps * ( coslat1*coslat1*cosazm*cosazm + sinlat1*sinlat1 );
351       b0 = (v1/(1.0+eps0)) * sqrt(1.0+eps*coslat1*coslat1*cosazm*cosazm);
352
353       x2 = v2*coslat2*(coslon2*coslon1+sinlon2*sinlon1);
354       y2 = v2*coslat2*(sinlon2*coslon1-coslon2*sinlon1);
355       z2 = v2*(1.0-e2)*sinlat2;
356       z1 = v1*(1.0-e2)*sinlat1;
357
358       c0 = c00 + c01*eps0 + c02*eps0*eps0 + c03*eps0*eps0*eps0;
359       c2 =       c21*eps0 + c22*eps0*eps0 + c23*eps0*eps0*eps0;
360       c4 =                  c42*eps0*eps0 + c43*eps0*eps0*eps0;
361       c6 =                                  c63*eps0*eps0*eps0;
362
363       bottom = cosazm*sqrt(1.0+eps0);
364       u1p = atan2(tanlat1,bottom);
365
366       top = v1*sinlat1+(1.0+eps0)*(z2-z1);
367       bottom = (x2*cosazm-y2*sinlat1*sinazm)*sqrt(1.0+eps0);
368       u2p = atan2(top,bottom);
369
370       aa = c0*(u2p-u1p);
371       bb = c2*(sin(2.0*u2p)-sin(2.0*u1p));
372       cc = c4*(sin(4.0*u2p)-sin(4.0*u1p));
373       dd = c6*(sin(6.0*u2p)-sin(6.0*u1p));
374
375       xdist = fabs(b0*(aa+bb+cc+dd));
376       *dist = xdist;
377
378/* compute the distance in degrees                                */
379                thg=atan(onemec2*tan(lat1rad));
380                dd2=sin(lon1rad);
381                ee2=-cos(lon1rad);
382                ff2=-cos(thg);
383                cc2=sin(thg);
384                aa2= ff2*ee2;
385                bb2=-ff2*dd2;
386
387/* -- Calculate some trig constants.          */
388
389                thg=atan(onemec2*tan(lat2rad));
390                dd3=sin(lon2rad);
391                ee3=-cos(lon2rad);
392                ff3=-cos(thg);
393                cc3=sin(thg);
394                aa3=ff3*ee3;
395                bb3=-ff3*dd3;
396                sc=aa2*aa3+bb2*bb3+cc2*cc3;
397
398/* - Spherical trig relationships used to compute angles.       */
399
400        aminus = aa2 - aa3;
401                bminus = bb2 - bb3;
402                cminus = cc2 - cc3;
403                aplus  = aa2 + aa3;
404                bplus  = bb2 + bb3;
405                cplus  = cc2 + cc3;
406
407                sd=0.5*sqrt((aminus*aminus + bminus*bminus+cminus*cminus)*(aplus*aplus
408                        + bplus*bplus + cplus*cplus));
409                deg=atan2(sd,sc)*todeg;
410                *xdeg = deg;
411
412           return(1);
413}
414
415void distaz_geo(olat,olon,dist,azimuth,lat,lon)
416double olat, olon, azimuth, dist, *lat, *lon;
417{
418
419/*
420=====================================================================
421 PURPOSE:  To compute latitude and longitude given a distance and
422            azimuth from a geographic point (olat,olon).
423=====================================================================
424 INPUT ARGUMENTS:
425    olat:    latitude in decimal degrees, North positive.
426    olon:    longitude in decimal degrees, East positive.
427    dist:    distance from point (olat,olon) in kms.
428    azimuth: azimuath between point (olat,olon) and (lat,lon) in degrees.
429=====================================================================
430 OUTPUT ARGUMENTS:
431    lat:    latitude in decimal degrees, North positive.
432    lon:    longitude in decimal degrees, East positive.
433=====================================================================
434*/
435
436      double C, C_sqr;
437      double torad, todeg;
438      double aa, bb, a0, b0;
439      double g0, g2, g4, g6;
440      double latrad, lonrad, azmrad;
441      double ec2, eps, eps0,       fl;
442      double tanmu;
443      double tanlat, coslat, sinlat, tanazm, sinazm, cosazm;
444      double v1, u1p, u2p, u1pbot, x2, y2, z2, ztop, zbot;
445      double temp,        mu, sd, sig1      ;
446      /* double sig2, eps2, delta; */
447
448/*
449Calculations are based upon the reference spheroid of 1967.
450See GEODESY by Bomford for definitions and reference values.
451Reference spheriod is found in GEODESY by Bomford (page 426).  Definitions
452for flattening, eccentricity and second eccentricity are found in
453Appendix A (page 646).
454*/
455
456      const double a = semi_major;
457      const double b = semi_minor;
458      double a1=0.25, b1=-0.125, c1=-0.0078125;
459      double g00=1.0, g01=-0.25, g02=0.109375, g03=-0.058139535;
460      double g21=0.125, g22=-0.06250, g23=0.034667969;
461      double g42=0.01953125, g43=-0.01953125;
462      double g63=0.004720052;
463
464/* - Initialize.             */
465
466      fl = ( a - b ) / a;         /* earth flattening                 */
467      ec2 = 2.0*fl-fl*fl;         /* square of eccentricity           */
468      eps = ec2/(1.0-ec2);        /* second eccentricity e'*e' = eps  */
469      torad = PI / 180.0;
470      todeg = 1.0 / torad;
471
472/* - Convert event location to radians.                               */
473
474      temp=olat;
475      if(temp == 0.) temp=1.0e-08;
476      latrad=torad*temp;
477      lonrad=torad*olon;
478      azmrad=azimuth*torad;
479
480/*  Compute some of the easier terms               */
481
482      coslat = cos(latrad);
483      sinlat = sin(latrad);
484      cosazm = cos(azmrad);
485      sinazm = sin(azmrad);
486      tanazm = tan(azmrad);
487      tanlat = tan(latrad);
488
489      C_sqr = coslat*coslat*cosazm*cosazm+ sinlat*sinlat;
490      C = sqrt(C_sqr);
491      eps0 = C_sqr * eps;
492      v1 = a / sqrt( 1.0 - ec2*sinlat*sinlat ); /* Radii of curvature */
493      b0 = v1*sqrt(1.0+eps*coslat*coslat*cosazm*cosazm)/(1.0+eps0);
494
495      g0 = g00 + g01*eps0 + g02*eps0*eps0 + g03*eps0*eps0*eps0;
496      g2 =       g21*eps0 + g22*eps0*eps0 + g23*eps0*eps0*eps0;
497      g4 =                  g42*eps0*eps0 + g43*eps0*eps0*eps0;
498      g6 =                                  g63*eps0*eps0*eps0;
499
500
501      u1pbot = cosazm * sqrt(1.0+eps0);
502      u1p=atan2(tanlat,u1pbot);
503      sig1 = 0.5*( 2.0*u1p-(a1*eps0+b1*eps0*eps0)*sin(2.0*u1p)+
504                 c1*eps0*eps0*sin(4.0*u1p));
505/*
506      aa = ( sig2 - sig1 ) page 117, GEODESY
507      bb = ( sig1 + sig2 ) page 117, GEODESY
508*/
509      aa =  (dist*g0)/b0;
510      bb = 2.0 * sig1 + aa;
511      u2p=u1p+aa+2.0*g2*sin(aa)*cos(bb)+2.0*g4*sin(2.*aa)*cos(2.*bb)+
512          2.*g6*sin(3.*aa)*cos(3.*bb);
513
514/*  This calculation of latitude is based on Rudoe's formulation, which
515    has not been tested  */
516
517/*    sinu1=tanlat/sqrt(1.0+eps+tanlat*tanlat);                       */
518/*    sinu2 = ((b0*C)/b)*sin(u2p)-((eps-eps0)/(1+eps0))*sinu1;        */
519/*    u2 = asin(sinu2);                                               */
520/*    arg = sinu2/sqrt(1.0-eps*eps*cos(u2)*cos(u2));                  */
521/*    *lat = asin(arg)*todeg;                                         */
522
523      a0 = b0*sqrt(1.0+eps0);
524      tanmu = sinlat*tanazm;
525      mu = atan(tanmu);
526
527/*  This calculation of longitude is based on Ruloe's formulation, which
528    has not been tested  */
529
530/*    arg = (a0*cos(u2p))/(a*cos(u2));                                */
531/*    delta = acos(arg)-mu;                                           */
532/*    *lon = (lonrad + delta)*todeg;                                  */
533
534/*
535   This calculation of latitude and longitude is an alternative to Rudoe's
536   formulation.  See GEODESY by Bomford (page 118).
537*/
538      sd=(ec2*v1*sinlat*coslat*sinazm)/(1.-ec2*coslat*coslat*sinazm*sinazm);
539      x2=a0*cos(u2p)*cos(mu)+b0*sin(u2p)*sin(mu)*coslat*sinazm+
540         sd*sinlat*sinazm;
541      y2=-a0*cos(u2p)*sin(mu)+b0*sin(u2p)*cos(mu)*coslat*sinazm+
542         sd*cosazm;
543      z2=b0*C*sin(u2p)-(sd*coslat*sinazm)/(1.0+eps);
544      ztop = (1.0+eps)*z2;
545      zbot = sqrt(x2*x2+y2*y2);
546      *lat = atan2(ztop,zbot)*todeg;
547      *lon = (atan(y2/x2)+lonrad)*todeg;
548      if((atan(y2/x2)+lonrad)<0.0) *lon=((atan(y2/x2)+lonrad)*todeg);
549}
Note: See TracBrowser for help on using the repository browser.