Skip to content

Commit

Permalink
fixes utm_geo errors in formula for utm/geo conversions (a 20-year ol…
Browse files Browse the repository at this point in the history 
…d mistake, finally fixed...)
  • Loading branch information
@danielpeter
danielpeter committed Apr 30, 2021
1 parent c7590b6 commit b994d55
Show file tree
Hide file tree
Showing 10 changed files with 222 additions and 149 deletions.
4 changes: 2 additions & 2 deletions EXAMPLES/Mount_StHelens/convert_lonlat2utm.pl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -179,7 +179,7 @@ sub geo2utm {
$f1 = (1. - $e2/4. - 3.*$e4/64. - 5.*$e6/256)*$rlat;
$f2 = 3.*$e2/8. + 3.*$e4/32. + 45.*$e6/1024.;
$f2 = $f2*sin(2.*$rlat);
$f3 = 15.*$e4/256.*45.*$e6/1024.;
$f3 = 15.*$e4/256. + 45.*$e6/1024.;
$f3 = $f3*sin(4.*$rlat);
$f4 = 35.*$e6/3072.;
$f4 = $f4*sin(6.*$rlat);
Expand Down Expand Up @@ -244,7 +244,7 @@ sub geo2utm {

$f1 = $rn1*tan($rlat1)/$r1;
$f2 = $d**2/2.;
$f3 = 5.*3.*$t1 + 10.*$c1 - 4.*$c1**2 - 9.*$ep2;
$f3 = 5. + 3.*$t1 + 10.*$c1 - 4.*$c1**2 - 9.*$ep2;
$f3 = $f3*$d**2*$d**2/24.;
$f4 = 61. + 90.*$t1 + 298.*$c1 + 45.*$t1**2. - 252.*$ep2 - 3.*$c1**2;
$f4 = $f4*($d**2)**3./720.;
Expand Down
171 changes: 86 additions & 85 deletions src/inverse_problem_for_source/CMT3D/cmt3d/utm_geo.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,8 @@ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway,SUPPRESS_UTM_PROJECT
logical SUPPRESS_UTM_PROJECTION

double precision, parameter :: degrad=PI/180.d0, raddeg=180.d0/PI
double precision, parameter :: semimaj=6378206.4d0, semimin=6356583.8d0
!double precision, parameter :: semimaj=6378206.4d0, semimin=6356583.8d0 ! Clarke 1866 ellipsoid
double precision, parameter :: semimaj=6378137.0d0, semimin=6356752.314245d0 ! WGS-84 ellipsoid
double precision, parameter :: scfa=0.9996d0
double precision, parameter :: north=0.d0, east=500000.d0

Expand Down Expand Up @@ -90,96 +91,96 @@ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway,SUPPRESS_UTM_PROJECT
!
if (iway == ILONGLAT2UTM) then

rlon = degrad*dlon
rlat = degrad*dlat

delam = dlon - cm
if (delam < -180.) delam = delam + 360.
if (delam > 180.) delam = delam - 360.
delam = delam*degrad

f1 = (1. - e2/4. - 3.*e4/64. - 5.*e6/256)*rlat
f2 = 3.*e2/8. + 3.*e4/32. + 45.*e6/1024.
f2 = f2*sin(2.*rlat)
f3 = 15.*e4/256.*45.*e6/1024.
f3 = f3*sin(4.*rlat)
f4 = 35.*e6/3072.
f4 = f4*sin(6.*rlat)
rm = semimaj*(f1 - f2 + f3 - f4)
if (dlat == 90. .or. dlat == -90.) then
xx = 0.
yy = scfa*rm
else
rn = semimaj/sqrt(1. - e2*sin(rlat)**2)
t = tan(rlat)**2
c = ep2*cos(rlat)**2
a = cos(rlat)*delam

f1 = (1. - t + c)*a**3/6.
f2 = 5. - 18.*t + t**2 + 72.*c - 58.*ep2
f2 = f2*a**5/120.
xx = scfa*rn*(a + f1 + f2)
f1 = a**2/2.
f2 = 5. - t + 9.*c + 4.*c**2
f2 = f2*a**4/24.
f3 = 61. - 58.*t + t**2 + 600.*c - 330.*ep2
f3 = f3*a**6/720.
yy = scfa*(rm + rn*tan(rlat)*(f1 + f2 + f3))
endif
xx = xx + east
yy = yy + north
rlon = degrad*dlon
rlat = degrad*dlat

delam = dlon - cm
if (delam < -180.) delam = delam + 360.
if (delam > 180.) delam = delam - 360.
delam = delam*degrad

f1 = (1. - e2/4. - 3.*e4/64. - 5.*e6/256)*rlat
f2 = 3.*e2/8. + 3.*e4/32. + 45.*e6/1024.
f2 = f2*sin(2.*rlat)
f3 = 15.*e4/256. + 45.*e6/1024.
f3 = f3*sin(4.*rlat)
f4 = 35.*e6/3072.
f4 = f4*sin(6.*rlat)
rm = semimaj*(f1 - f2 + f3 - f4)
if (dlat == 90. .or. dlat == -90.) then
xx = 0.
yy = scfa*rm
else
rn = semimaj/sqrt(1. - e2*sin(rlat)**2)
t = tan(rlat)**2
c = ep2*cos(rlat)**2
a = cos(rlat)*delam

f1 = (1. - t + c)*a**3/6.
f2 = 5. - 18.*t + t**2 + 72.*c - 58.*ep2
f2 = f2*a**5/120.
xx = scfa*rn*(a + f1 + f2)
f1 = a**2/2.
f2 = 5. - t + 9.*c + 4.*c**2
f2 = f2*a**4/24.
f3 = 61. - 58.*t + t**2 + 600.*c - 330.*ep2
f3 = f3*a**6/720.
yy = scfa*(rm + rn*tan(rlat)*(f1 + f2 + f3))
endif
xx = xx + east
yy = yy + north

!
!---- UTM to Lat/Lon conversion
!
else

xx = xx - east
yy = yy - north
e1 = sqrt(1. - e2)
e1 = (1. - e1)/(1. + e1)
rm = yy/scfa
u = 1. - e2/4. - 3.*e4/64. - 5.*e6/256.
u = rm/(semimaj*u)

f1 = 3.*e1/2. - 27.*e1**3./32.
f1 = f1*sin(2.*u)
f2 = 21.*e1**2/16. - 55.*e1**4/32.
f2 = f2*sin(4.*u)
f3 = 151.*e1**3./96.
f3 = f3*sin(6.*u)
rlat1 = u + f1 + f2 + f3
dlat1 = rlat1*raddeg
if (dlat1 >= 90. .or. dlat1 <= -90.) then
dlat1 = dmin1(dlat1,dble(90.) )
dlat1 = dmax1(dlat1,dble(-90.) )
dlon = cm
else
c1 = ep2*cos(rlat1)**2.
t1 = tan(rlat1)**2.
f1 = 1. - e2*sin(rlat1)**2.
rn1 = semimaj/sqrt(f1)
r1 = semimaj*(1. - e2)/sqrt(f1**3)
d = xx/(rn1*scfa)

f1 = rn1*tan(rlat1)/r1
f2 = d**2/2.
f3 = 5.*3.*t1 + 10.*c1 - 4.*c1**2 - 9.*ep2
f3 = f3*d**2*d**2/24.
f4 = 61. + 90.*t1 + 298.*c1 + 45.*t1**2. - 252.*ep2 - 3.*c1**2
f4 = f4*(d**2)**3./720.
rlat = rlat1 - f1*(f2 - f3 + f4)
dlat = rlat*raddeg

f1 = 1. + 2.*t1 + c1
f1 = f1*d**2*d/6.
f2 = 5. - 2.*c1 + 28.*t1 - 3.*c1**2 + 8.*ep2 + 24.*t1**2.
f2 = f2*(d**2)**2*d/120.
rlon = cmr + (d - f1 + f2)/cos(rlat1)
dlon = rlon*raddeg
if (dlon < -180.) dlon = dlon + 360.
if (dlon > 180.) dlon = dlon - 360.
endif
xx = xx - east
yy = yy - north
e1 = sqrt(1. - e2)
e1 = (1. - e1)/(1. + e1)
rm = yy/scfa
u = 1. - e2/4. - 3.*e4/64. - 5.*e6/256.
u = rm/(semimaj*u)

f1 = 3.*e1/2. - 27.*e1**3./32.
f1 = f1*sin(2.*u)
f2 = 21.*e1**2/16. - 55.*e1**4/32.
f2 = f2*sin(4.*u)
f3 = 151.*e1**3./96.
f3 = f3*sin(6.*u)
rlat1 = u + f1 + f2 + f3
dlat1 = rlat1*raddeg
if (dlat1 >= 90. .or. dlat1 <= -90.) then
dlat1 = dmin1(dlat1,dble(90.) )
dlat1 = dmax1(dlat1,dble(-90.) )
dlon = cm
else
c1 = ep2*cos(rlat1)**2.
t1 = tan(rlat1)**2.
f1 = 1. - e2*sin(rlat1)**2.
rn1 = semimaj/sqrt(f1)
r1 = semimaj*(1. - e2)/sqrt(f1**3)
d = xx/(rn1*scfa)

f1 = rn1*tan(rlat1)/r1
f2 = d**2/2.
f3 = 5. + 3.*t1 + 10.*c1 - 4.*c1**2 - 9.*ep2
f3 = f3*d**2*d**2/24.
f4 = 61. + 90.*t1 + 298.*c1 + 45.*t1**2. - 252.*ep2 - 3.*c1**2
f4 = f4*(d**2)**3./720.
rlat = rlat1 - f1*(f2 - f3 + f4)
dlat = rlat*raddeg

f1 = 1. + 2.*t1 + c1
f1 = f1*d**2*d/6.
f2 = 5. - 2.*c1 + 28.*t1 - 3.*c1**2 + 8.*ep2 + 24.*t1**2.
f2 = f2*(d**2)**2*d/120.
rlon = cmr + (d - f1 + f2)/cos(rlat1)
dlon = rlon*raddeg
if (dlon < -180.) dlon = dlon + 360.
if (dlon > 180.) dlon = dlon - 360.
endif
endif

if (iway == IUTM2LONGLAT) then
Expand Down
Loading

0 comments on commit b994d55

Please sign in to comment.