184 subroutine direct(a, f, lat1, lon1, azi1, s12a12, flags,
185 + lat2, lon2, azi2, omask, a12s12, m12, mm12, mm21, ss12)
187 double precision a, f, lat1, lon1, azi1, s12a12
190 double precision lat2, lon2, azi2
192 double precision a12s12, m12, MM12, MM21, SS12
194 integer ord, nC1, nC1p, nC2, nA3, nA3x, nC3, nC3x, nC4, nC4x
195 parameter(ord = 6, nc1 = ord, nc1p = ord,
196 + nc2 = ord, na3 = ord, na3x = na3,
197 + nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2,
198 + nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
199 double precision A3x(0:na3x-1), C3x(0:nc3x-1), C4x(0:nc4x-1),
200 + c1a(nc1), c1pa(nc1p), c2a(nc2), c3a(nc3-1), c4a(0:nc4-1)
202 double precision csmgt, atanhx, hypotx,
203 + angnm, angnm2, angrnd, trgsum, a1m1f, a2m1f, a3f
204 logical arcmod, unroll, arcp, redlp, scalp, areap
205 double precision e2, f1, ep2, n, b, c2,
206 + lon1x, azi1x, phi, alp1, salp0, calp0, k2, eps,
207 + salp1, calp1, ssig1, csig1, cbet1, sbet1, dn1, somg1, comg1,
208 + salp2, calp2, ssig2, csig2, sbet2, cbet2, dn2, somg2, comg2,
209 + ssig12, csig12, salp12, calp12, omg12, lam12, lon12,
210 + sig12, stau1, ctau1, tau12, s12a, t, s, c, serr, e,
211 + a1m1, a2m1, a3c, a4, ab1, ab2,
212 + b11, b12, b21, b22, b31, b41, b42, j12
214 double precision dblmin, dbleps, pi, degree, tiny,
215 + tol0, tol1, tol2, tolb, xthrsh
216 integer digits, maxit1, maxit2
218 common /geocom/ dblmin, dbleps, pi, degree, tiny,
219 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
221 if (.not.init)
call geoini
230 arcmod = mod(flags/1, 2) == 1
231 unroll = mod(flags/2, 2) == 1
233 arcp = mod(omask/1, 2) == 1
234 redlp = mod(omask/2, 2) == 1
235 scalp = mod(omask/4, 2) == 1
236 areap = mod(omask/8, 2) == 1
241 else if (e2 .gt. 0)
then 242 c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
244 c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
250 if (areap)
call c4cof(n, c4x)
253 azi1x = angrnd(angnm(azi1))
257 alp1 = azi1x * degree
260 salp1 = csmgt(0d0, sin(alp1), azi1x .eq. -180)
261 calp1 = csmgt(0d0, cos(alp1), abs(azi1x) .eq. 90)
265 sbet1 = f1 * sin(phi)
266 cbet1 = csmgt(tiny, cos(phi), abs(lat1) .eq. 90)
267 call norm2(sbet1, cbet1)
268 dn1 = sqrt(1 + ep2 * sbet1**2)
272 salp0 = salp1 * cbet1
275 calp0 = hypotx(calp1, salp1 * sbet1)
286 somg1 = salp0 * sbet1
287 csig1 = csmgt(cbet1 * calp1, 1d0, sbet1 .ne. 0 .or. calp1 .ne. 0)
290 call norm2(ssig1, csig1)
294 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
298 b11 = trgsum(.true., ssig1, csig1, c1a, nc1)
302 stau1 = ssig1 * c + csig1 * s
303 ctau1 = csig1 * c - ssig1 * s
307 if (.not. arcmod)
call c1pf(eps, c1pa)
309 if (redlp .or. scalp)
then 312 b21 = trgsum(.true., ssig1, csig1, c2a, nc2)
319 call c3f(eps, c3x, c3a)
320 a3c = -f * salp0 * a3f(eps, a3x)
321 b31 = trgsum(.true., ssig1, csig1, c3a, nc3-1)
324 call c4f(eps, c4x, c4a)
326 a4 = a**2 * calp0 * salp0 * e2
327 b41 = trgsum(.false., ssig1, csig1, c4a, nc4)
336 sig12 = s12a12 * degree
338 s12a = s12a - 180 * aint(s12a / 180)
339 ssig12 = csmgt(0d0, sin(sig12), s12a .eq. 0)
340 csig12 = csmgt(0d0, cos(sig12), s12a .eq. 90)
345 tau12 = s12a12 / (b * (1 + a1m1))
349 b12 = - trgsum(.true.,
350 + stau1 * c + ctau1 * s, ctau1 * c - stau1 * s, c1pa, nc1p)
351 sig12 = tau12 - (b12 - b11)
354 if (abs(f) .gt. 0.01d0)
then 376 ssig2 = ssig1 * csig12 + csig1 * ssig12
377 csig2 = csig1 * csig12 - ssig1 * ssig12
378 b12 = trgsum(.true., ssig2, csig2, c1a, nc1)
379 serr = (1 + a1m1) * (sig12 + (b12 - b11)) - s12a12 / b
380 sig12 = sig12 - serr / sqrt(1 + k2 * ssig2**2)
388 ssig2 = ssig1 * csig12 + csig1 * ssig12
389 csig2 = csig1 * csig12 - ssig1 * ssig12
390 dn2 = sqrt(1 + k2 * ssig2**2)
391 if (arcmod .or. abs(f) .gt. 0.01d0)
392 + b12 = trgsum(.true., ssig2, csig2, c1a, nc1)
393 ab1 = (1 + a1m1) * (b12 - b11)
396 sbet2 = calp0 * ssig2
398 cbet2 = hypotx(salp0, calp0 * csig2)
399 if (cbet2 .eq. 0)
then 406 somg2 = salp0 * ssig2
411 calp2 = calp0 * csig2
415 omg12 = csmgt(e * (sig12
416 + - (atan2( ssig2, csig2) - atan2( ssig1, csig1))
417 + + (atan2(e * somg2, comg2) - atan2(e * somg1, comg1))),
418 + atan2(somg2 * comg1 - comg2 * somg1,
419 + comg2 * comg1 + somg2 * somg1),
422 lam12 = omg12 + a3c *
423 + ( sig12 + (trgsum(.true., ssig2, csig2, c3a, nc3-1)
425 lon12 = lam12 / degree
428 lon2 = csmgt(lon1 + lon12, angnm(lon1x + angnm2(lon12)), unroll)
429 lat2 = atan2(sbet2, f1 * cbet2) / degree
431 azi2 = 0 - atan2(-salp2, calp2) / degree
433 if (redlp .or. scalp)
then 434 b22 = trgsum(.true., ssig2, csig2, c2a, nc2)
435 ab2 = (1 + a2m1) * (b22 - b21)
436 j12 = (a1m1 - a2m1) * sig12 + (ab1 - ab2)
440 if (redlp) m12 = b * ((dn2 * (csig1 * ssig2) -
441 + dn1 * (ssig1 * csig2)) - csig1 * csig2 * j12)
443 t = k2 * (ssig2 - ssig1) * (ssig2 + ssig1) / (dn1 + dn2)
444 mm12 = csig12 + (t * ssig2 - csig2 * j12) * ssig1 / dn1
445 mm21 = csig12 - (t * ssig1 - csig1 * j12) * ssig2 / dn2
449 b42 = trgsum(.false., ssig2, csig2, c4a, nc4)
450 if (calp0 .eq. 0 .or. salp0 .eq. 0)
then 452 salp12 = salp2 * calp1 - calp2 * salp1
453 calp12 = calp2 * calp1 + salp2 * salp1
457 if (salp12 .eq. 0 .and. calp12 .lt. 0)
then 458 salp12 = tiny * calp1
470 salp12 = calp0 * salp0 *
471 + csmgt(csig1 * (1 - csig12) + ssig12 * ssig1,
472 + ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1),
474 calp12 = salp0**2 + calp0**2 * csig1 * csig2
476 ss12 = c2 * atan2(salp12, calp12) + a4 * (b42 - b41)
479 if (arcp) a12s12 = csmgt(b * ((1 + a1m1) * sig12 + ab1),
480 + sig12 / degree, arcmod)
528 subroutine invers(a, f, lat1, lon1, lat2, lon2,
529 + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
531 double precision a, f, lat1, lon1, lat2, lon2
534 double precision s12, azi1, azi2
536 double precision a12, m12, MM12, MM21, SS12
538 integer ord, nC1, nC2, nA3, nA3x, nC3, nC3x, nC4, nC4x
539 parameter(ord = 6, nc1 = ord, nc2 = ord, na3 = ord, na3x = na3,
540 + nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2,
541 + nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
542 double precision A3x(0:na3x-1), C3x(0:nc3x-1), C4x(0:nc4x-1),
543 + c1a(nc1), c2a(nc2), c3a(nc3-1), c4a(0:nc4-1)
545 double precision csmgt, atanhx, hypotx,
546 + angnm, angdif, angrnd, trgsum, lam12f, invsta
547 integer latsgn, lonsgn, swapp, numit
548 logical arcp, redlp, scalp, areap, merid, tripn, tripb
550 double precision e2, f1, ep2, n, b, c2,
551 + lat1x, lat2x, phi, salp0, calp0, k2, eps,
552 + salp1, calp1, ssig1, csig1, cbet1, sbet1, dbet1, dn1,
553 + salp2, calp2, ssig2, csig2, sbet2, cbet2, dbet2, dn2,
554 + slam12, clam12, salp12, calp12, omg12, lam12, lon12,
555 + salp1a, calp1a, salp1b, calp1b,
556 + dalp1, sdalp1, cdalp1, nsalp1, alp12, somg12, domg12,
557 + sig12, v, dv, dnm, dummy,
558 + a4, b41, b42, s12x, m12x, a12x
560 double precision dblmin, dbleps, pi, degree, tiny,
561 + tol0, tol1, tol2, tolb, xthrsh
562 integer digits, maxit1, maxit2
564 common /geocom/ dblmin, dbleps, pi, degree, tiny,
565 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
567 if (.not.init)
call geoini
576 arcp = mod(omask/1, 2) == 1
577 redlp = mod(omask/2, 2) == 1
578 scalp = mod(omask/4, 2) == 1
579 areap = mod(omask/8, 2) == 1
584 else if (e2 .gt. 0)
then 585 c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
587 c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
593 if (areap)
call c4cof(n, c4x)
598 lon12 = angdif(angnm(lon1), angnm(lon2))
600 lon12 = angrnd(lon12)
602 if (lon12 .ge. 0)
then 607 lon12 = lon12 * lonsgn
612 if (abs(lat1x) .ge. abs(lat2x))
then 617 if (swapp .lt. 0)
then 619 call swap(lat1x, lat2x)
622 if (lat1x .lt. 0)
then 627 lat1x = lat1x * latsgn
628 lat2x = lat2x * latsgn
643 sbet1 = f1 * sin(phi)
644 cbet1 = csmgt(tiny, cos(phi), lat1x .eq. -90)
645 call norm2(sbet1, cbet1)
649 sbet2 = f1 * sin(phi)
650 cbet2 = csmgt(tiny, cos(phi), abs(lat2x) .eq. 90)
651 call norm2(sbet2, cbet2)
662 if (cbet1 .lt. -sbet1)
then 663 if (cbet2 .eq. cbet1) sbet2 = sign(sbet1, sbet2)
665 if (abs(sbet2) .eq. -sbet1) cbet2 = cbet1
668 dn1 = sqrt(1 + ep2 * sbet1**2)
669 dn2 = sqrt(1 + ep2 * sbet2**2)
671 lam12 = lon12 * degree
673 if (lon12 .eq. 180) slam12 = 0
679 merid = lat1x .eq. -90 .or. slam12 == 0
695 csig1 = calp1 * cbet1
697 csig2 = calp2 * cbet2
700 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, 0d0),
701 + csig1 * csig2 + ssig1 * ssig2)
702 call lengs(n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
703 + cbet1, cbet2, s12x, m12x, dummy,
704 + scalp, mm12, mm21, ep2, c1a, c2a)
713 if (sig12 .lt. 1 .or. m12x .ge. 0)
then 716 a12x = sig12 / degree
724 if (.not. merid .and. sbet1 .eq. 0 .and.
725 + (f .le. 0 .or. lam12 .le. pi - f * pi))
then 735 m12x = b * sin(sig12)
741 else if (.not. merid)
then 746 sig12 = invsta(sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12,
747 + f, a3x, salp1, calp1, salp2, calp2, dnm, c1a, c2a)
749 if (sig12 .ge. 0)
then 751 s12x = sig12 * b * dnm
752 m12x = dnm**2 * b * sin(sig12 / dnm)
754 mm12 = cos(sig12 / dnm)
757 a12x = sig12 / degree
758 omg12 = lam12 / (f1 * dnm)
780 do 10 numit = 0, maxit2-1
783 v = lam12f(sbet1, cbet1, dn1, sbet2, cbet2, dn2,
784 + salp1, calp1, f, a3x, c3x, salp2, calp2, sig12,
785 + ssig1, csig1, ssig2, csig2,
786 + eps, omg12, numit .lt. maxit1, dv,
787 + c1a, c2a, c3a) - lam12
791 + .not. (abs(v) .ge. csmgt(8d0, 2d0, tripn) * tol0))
794 if (v .gt. 0 .and. (numit .gt. maxit1 .or.
795 + calp1/salp1 .gt. calp1b/salp1b))
then 798 else if (v .lt. 0 .and. (numit .gt. maxit1 .or.
799 + calp1/salp1 .lt. calp1a/salp1a))
then 803 if (numit .lt. maxit1 .and. dv .gt. 0)
then 807 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1
808 if (nsalp1 .gt. 0 .and. abs(dalp1) .lt. pi)
then 809 calp1 = calp1 * cdalp1 - salp1 * sdalp1
811 call norm2(salp1, calp1)
815 tripn = abs(v) .le. 16 * tol0
827 salp1 = (salp1a + salp1b)/2
828 calp1 = (calp1a + calp1b)/2
829 call norm2(salp1, calp1)
831 tripb = abs(salp1a - salp1) + (calp1a - calp1) .lt. tolb
832 + .or. abs(salp1 - salp1b) + (calp1 - calp1b) .lt. tolb
835 call lengs(eps, sig12, ssig1, csig1, dn1,
836 + ssig2, csig2, dn2, cbet1, cbet2, s12x, m12x, dummy,
837 + scalp, mm12, mm21, ep2, c1a, c2a)
840 a12x = sig12 / degree
841 omg12 = lam12 - omg12
847 if (redlp) m12 = 0 + m12x
851 salp0 = salp1 * cbet1
852 calp0 = hypotx(calp1, salp1 * sbet1)
853 if (calp0 .ne. 0 .and. salp0 .ne. 0)
then 856 csig1 = calp1 * cbet1
858 csig2 = calp2 * cbet2
860 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
862 a4 = a**2 * calp0 * salp0 * e2
863 call norm2(ssig1, csig1)
864 call norm2(ssig2, csig2)
865 call c4f(eps, c4x, c4a)
866 b41 = trgsum(.false., ssig1, csig1, c4a, nc4)
867 b42 = trgsum(.false., ssig2, csig2, c4a, nc4)
868 ss12 = a4 * (b42 - b41)
874 if (.not. merid .and. omg12 .lt. 0.75d0 * pi
875 + .and. sbet2 - sbet1 .lt. 1.75d0)
then 880 domg12 = 1 + cos(omg12)
883 alp12 = 2 * atan2(somg12 * (sbet1 * dbet2 + sbet2 * dbet1),
884 + domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) )
887 salp12 = salp2 * calp1 - calp2 * salp1
888 calp12 = calp2 * calp1 + salp2 * salp1
893 if (salp12 .eq. 0 .and. calp12 .lt. 0)
then 894 salp12 = tiny * calp1
897 alp12 = atan2(salp12, calp12)
899 ss12 = ss12 + c2 * alp12
900 ss12 = ss12 * swapp * lonsgn * latsgn
906 if (swapp .lt. 0)
then 907 call swap(salp1, salp2)
908 call swap(calp1, calp2)
909 if (scalp)
call swap(mm12, mm21)
912 salp1 = salp1 * swapp * lonsgn
913 calp1 = calp1 * swapp * latsgn
914 salp2 = salp2 * swapp * lonsgn
915 calp2 = calp2 * swapp * latsgn
918 azi1 = 0 - atan2(-salp1, calp1) / degree
919 azi2 = 0 - atan2(-salp2, calp2) / degree
945 subroutine area(a, f, lats, lons, n, AA, PP)
948 double precision a, f, lats(n), lons(n)
950 double precision AA, PP
952 integer i, omask, cross, trnsit
953 double precision s12, azi1, azi2, dummy, SS12, b, e2, c2, area0,
954 + atanhx, aacc(2), pacc(2)
956 double precision dblmin, dbleps, pi, degree, tiny,
957 + tol0, tol1, tol2, tolb, xthrsh
958 integer digits, maxit1, maxit2
960 common /geocom/ dblmin, dbleps, pi, degree, tiny,
961 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
968 call invers(a, f, lats(i+1), lons(i+1),
969 + lats(mod(i+1,n)+1), lons(mod(i+1,n)+1),
970 + s12, azi1, azi2, omask, dummy, dummy, dummy, dummy, ss12)
971 call accadd(pacc, s12)
972 call accadd(aacc, -ss12)
973 cross = cross + trnsit(lons(i+1), lons(mod(i+1,n)+1))
980 else if (e2 .gt. 0)
then 981 c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
983 c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
986 if (mod(abs(cross), 2) .eq. 1)
then 987 if (aacc(1) .lt. 0)
then 988 call accadd(aacc, +area0/2)
990 call accadd(aacc, -area0/2)
993 if (aacc(1) .gt. area0/2)
then 994 call accadd(aacc, -area0)
995 else if (aacc(1) .le. -area0/2)
then 996 call accadd(aacc, +area0)
1006 double precision dblmin, dbleps, pi, degree, tiny,
1007 + tol0, tol1, tol2, tolb, xthrsh
1008 integer digits, maxit1, maxit2
1011 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1012 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1016 double precision dblmin, dbleps, pi, degree, tiny,
1017 + tol0, tol1, tol2, tolb, xthrsh
1018 integer digits, maxit1, maxit2
1020 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1021 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1024 dblmin = 0.5d0**1022
1025 dbleps = 0.5d0**(digits-1)
1027 pi = atan2(0d0, -1d0)
1038 xthrsh = 1000 * tol2
1040 maxit2 = maxit1 + digits + 10
1047 subroutine lengs(eps, sig12,
1048 + ssig1, csig1, dn1, ssig2, csig2, dn2,
1049 + cbet1, cbet2, s12b, m12b, m0,
1050 + scalp, mm12, mm21, ep2, c1a, c2a)
1052 double precision eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1056 double precision s12b, m12b, m0
1058 double precision MM12, MM21
1060 double precision C1a(*), C2a(*)
1062 integer ord, nC1, nC2
1063 parameter(ord = 6, nc1 = ord, nc2 = ord)
1065 double precision A1m1f, A2m1f, TrgSum
1066 double precision A1m1, AB1, A2m1, AB2, J12, csig12, t
1074 ab1 = (1 + a1m1) * (trgsum(.true., ssig2, csig2, c1a, nc1) -
1075 + trgsum(.true., ssig1, csig1, c1a, nc1))
1077 ab2 = (1 + a2m1) * (trgsum(.true., ssig2, csig2, c2a, nc2) -
1078 + trgsum(.true., ssig1, csig1, c2a, nc2))
1080 j12 = m0 * sig12 + (ab1 - ab2)
1084 m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
1085 + csig1 * csig2 * j12
1087 s12b = (1 + a1m1) * sig12 + ab1
1089 csig12 = csig1 * csig2 + ssig1 * ssig2
1090 t = ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2)
1091 mm12 = csig12 + (t * ssig2 - csig2 * j12) * ssig1 / dn1
1092 mm21 = csig12 - (t * ssig1 - csig1 * j12) * ssig2 / dn2
1098 double precision function astrd(x, y)
1102 double precision x, y
1104 double precision cbrt, csmgt
1105 double precision k, p, q, r, S, r2, r3, disc, u,
1106 + t3, t, ang, v, uv, w
1111 if ( .not. (q .eq. 0 .and. r .lt. 0) )
then 1120 disc = s * (s + 2 * r3)
1122 if (disc .ge. 0)
then 1128 t3 = t3 + csmgt(-sqrt(disc), sqrt(disc), t3 .lt. 0)
1133 if (t .ne. 0) u = u + t + r2 / t
1136 ang = atan2(sqrt(-disc), -(s + r3))
1139 u = u + 2 * r * cos(ang / 3)
1145 uv = csmgt(q / (v - u), u + v, u .lt. 0)
1147 w = (uv - q) / (2 * v)
1151 k = uv / (sqrt(uv + w**2) + w)
1163 double precision function invsta(sbet1, cbet1, dn1,
1164 + sbet2, cbet2, dn2, lam12, f, a3x,
1165 + salp1, calp1, salp2, calp2, dnm,
1171 double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12,
1174 double precision salp1, calp1, salp2, calp2, dnm
1176 double precision C1a(*), C2a(*)
1178 double precision csmgt, hypotx, A3f, Astrd
1180 double precision f1, e2, ep2, n, etol2, k2, eps, sig12,
1181 + sbet12, cbet12, sbt12a, omg12, somg12, comg12, ssig12, csig12,
1182 + x, y, lamscl, betscl, cbt12a, bt12a, m12b, m0, dummy,
1185 double precision dblmin, dbleps, pi, degree, tiny,
1186 + tol0, tol1, tol2, tolb, xthrsh
1187 integer digits, maxit1, maxit2
1189 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1190 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1206 etol2 = 0.1d0 * tol2 /
1207 + sqrt( max(0.001d0, abs(f)) * min(1d0, 1 - f/2) / 2 )
1212 sbet12 = sbet2 * cbet1 - cbet2 * sbet1
1213 cbet12 = cbet2 * cbet1 + sbet2 * sbet1
1214 sbt12a = sbet2 * cbet1 + cbet2 * sbet1
1216 shortp = cbet12 .ge. 0 .and. sbet12 .lt. 0.5d0 .and.
1217 + cbet2 * lam12 .lt. 0.5d0
1221 sbetm2 = (sbet1 + sbet2)**2
1224 sbetm2 = sbetm2 / (sbetm2 + (cbet1 + cbet2)**2)
1225 dnm = sqrt(1 + ep2 * sbetm2)
1226 omg12 = omg12 / (f1 * dnm)
1231 salp1 = cbet2 * somg12
1232 calp1 = csmgt(sbet12 + cbet2 * sbet1 * somg12**2 / (1 + comg12),
1233 + sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12),
1236 ssig12 = hypotx(salp1, calp1)
1237 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12
1239 if (shortp .and. ssig12 .lt. etol2)
then 1241 salp2 = cbet1 * somg12
1242 calp2 = sbet12 - cbet1 * sbet2 *
1243 + csmgt(somg12**2 / (1 + comg12), 1 - comg12, comg12 .ge. 0)
1244 call norm2(salp2, calp2)
1246 sig12 = atan2(ssig12, csig12)
1247 else if (abs(n) .gt. 0.1d0 .or. csig12 .ge. 0 .or.
1248 + ssig12 .ge. 6 * abs(n) * pi * cbet1**2)
then 1257 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
1258 lamscl = f * cbet1 * a3f(eps, a3x) * pi
1259 betscl = lamscl * cbet1
1260 x = (lam12 - pi) / lamscl
1264 cbt12a = cbet2 * cbet1 - sbet2 * sbet1
1265 bt12a = atan2(sbt12a, cbt12a)
1268 call lengs(n, pi + bt12a,
1269 + sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
1270 + cbet1, cbet2, dummy, m12b, m0, .false.,
1271 + dummy, dummy, ep2, c1a, c2a)
1272 x = -1 + m12b / (cbet1 * cbet2 * m0 * pi)
1273 betscl = csmgt(sbt12a / x, -f * cbet1**2 * pi,
1275 lamscl = betscl / cbet1
1276 y = (lam12 - pi) / lamscl
1279 if (y .gt. -tol1 .and. x .gt. -1 - xthrsh)
then 1282 salp1 = min(1d0, -x)
1283 calp1 = - sqrt(1 - salp1**2)
1285 calp1 = max(csmgt(0d0, 1d0, x .gt. -tol1), x)
1286 salp1 = sqrt(1 - calp1**2)
1325 + csmgt(-x * k/(1 + k), -y * (1 + k)/k, f .ge. 0)
1326 somg12 = sin(omg12a)
1327 comg12 = -cos(omg12a)
1329 salp1 = cbet2 * somg12
1330 calp1 = sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12)
1334 if (.not. (salp1 .le. 0))
then 1335 call norm2(salp1, calp1)
1345 double precision function lam12f(sbet1, cbet1, dn1,
1346 + sbet2, cbet2, dn2, salp1, calp1, f, a3x, c3x, salp2, calp2,
1347 + sig12, ssig1, csig1, ssig2, csig2, eps, domg12, diffp, dlam12,
1350 double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2,
1351 + salp1, calp1, f, a3x(*), c3x(*)
1354 double precision salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
1357 double precision dlam12
1359 double precision C1a(*), C2a(*), C3a(*)
1362 parameter(ord = 6, nc3 = ord)
1364 double precision csmgt, hypotx, A3f, TrgSum
1366 double precision f1, e2, ep2, salp0, calp0,
1367 + somg1, comg1, somg2, comg2, omg12, lam12, b312, h0, k2, dummy
1369 double precision dblmin, dbleps, pi, degree, tiny,
1370 + tol0, tol1, tol2, tolb, xthrsh
1371 integer digits, maxit1, maxit2
1373 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1374 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1381 if (sbet1 .eq. 0 .and. calp1 .eq. 0) calp1 = -tiny
1384 salp0 = salp1 * cbet1
1386 calp0 = hypotx(calp1, salp1 * sbet1)
1391 somg1 = salp0 * sbet1
1392 csig1 = calp1 * cbet1
1394 call norm2(ssig1, csig1)
1401 salp2 = csmgt(salp0 / cbet2, salp1, cbet2 .ne. cbet1)
1406 calp2 = csmgt(sqrt((calp1 * cbet1)**2 +
1407 + csmgt((cbet2 - cbet1) * (cbet1 + cbet2),
1408 + (sbet1 - sbet2) * (sbet1 + sbet2),
1409 + cbet1 .lt. -sbet1)) / cbet2,
1410 + abs(calp1), cbet2 .ne. cbet1 .or. abs(sbet2) .ne. -sbet1)
1414 somg2 = salp0 * sbet2
1415 csig2 = calp2 * cbet2
1417 call norm2(ssig2, csig2)
1421 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, 0d0),
1422 + csig1 * csig2 + ssig1 * ssig2)
1425 omg12 = atan2(max(comg1 * somg2 - somg1 * comg2, 0d0),
1426 + comg1 * comg2 + somg1 * somg2)
1428 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
1429 call c3f(eps, c3x, c3a)
1430 b312 = (trgsum(.true., ssig2, csig2, c3a, nc3-1) -
1431 + trgsum(.true., ssig1, csig1, c3a, nc3-1))
1432 h0 = -f * a3f(eps, a3x)
1433 domg12 = salp0 * h0 * (sig12 + b312)
1434 lam12 = omg12 + domg12
1437 if (calp2 .eq. 0)
then 1438 dlam12 = - 2 * f1 * dn1 / sbet1
1440 call lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1441 + cbet1, cbet2, dummy, dlam12, dummy,
1442 + .false., dummy, dummy, ep2, c1a, c2a)
1443 dlam12 = dlam12 * f1 / (calp2 * cbet2)
1451 double precision function a3f(eps, A3x)
1453 integer ord, nA3, nA3x
1454 parameter(ord = 6, na3 = ord, na3x = na3)
1457 double precision eps
1459 double precision A3x(0: na3x-1)
1461 double precision polval
1462 a3f = polval(na3 - 1, a3x, eps)
1467 subroutine c3f(eps, C3x, c)
1470 integer ord, nC3, nC3x
1471 parameter(ord = 6, nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2)
1474 double precision eps, C3x(0:nc3x-1)
1476 double precision c(nc3-1)
1479 double precision mult, polval
1483 do 10 l = 1, nc3 - 1
1486 c(l) = mult * polval(m, c3x(o), eps)
1493 subroutine c4f(eps, C4x, c)
1496 integer ord, nC4, nC4x
1497 parameter(ord = 6, nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
1500 double precision eps, C4x(0:nc4x-1)
1502 double precision c(0:nc4-1)
1505 double precision mult, polval
1509 do 10 l = 0, nc4 - 1
1511 c(l) = mult * polval(m, c4x(o), eps)
1519 double precision function a1m1f(eps)
1522 double precision eps
1525 integer ord, nA1, o, m
1526 parameter(ord = 6, na1 = ord)
1527 double precision polval, coeff(na1/2 + 2)
1528 data coeff /1, 4, 64, 0, 256/
1532 t = polval(m, coeff(o), eps**2) / coeff(o + m + 1)
1533 a1m1f = (t + eps) / (1 - eps)
1538 subroutine c1f(eps, c)
1541 parameter(ord = 6, nc1 = ord)
1544 double precision eps
1546 double precision c(nc1)
1548 double precision eps2, d
1550 double precision polval, coeff((nc1**2 + 7*nc1 - 2*(nc1/2))/4)
1553 + -9, 64, -128, 2048,
1564 c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
1572 subroutine c1pf(eps, c)
1575 parameter(ord = 6, nc1p = ord)
1578 double precision eps
1580 double precision c(nc1p)
1582 double precision eps2, d
1584 double precision polval, coeff((nc1p**2 + 7*nc1p - 2*(nc1p/2))/4)
1586 + 205, -432, 768, 1536,
1587 + 4005, -4736, 3840, 12288,
1589 + -7173, 2695, 7680,
1598 c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
1607 double precision function a2m1f(eps)
1609 double precision eps
1612 integer ord, nA2, o, m
1613 parameter(ord = 6, na2 = ord)
1614 double precision polval, coeff(na2/2 + 2)
1615 data coeff /25, 36, 64, 0, 256/
1619 t = polval(m, coeff(o), eps**2) / coeff(o + m + 1)
1620 a2m1f = t * (1 - eps) - eps
1625 subroutine c2f(eps, c)
1628 parameter(ord = 6, nc2 = ord)
1631 double precision eps
1633 double precision c(nc2)
1635 double precision eps2, d
1637 double precision polval, coeff((nc2**2 + 7*nc2 - 2*(nc2/2))/4)
1640 + 35, 64, 384, 2048,
1651 c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
1659 subroutine a3cof(n, A3x)
1661 integer ord, nA3, nA3x
1662 parameter(ord = 6, na3 = ord, na3x = na3)
1667 double precision A3x(0:na3x-1)
1670 double precision polval, coeff((na3**2 + 7*na3 - 2*(na3/2))/4)
1681 do 10 j = na3 - 1, 0, -1
1682 m = min(na3 - j - 1, j)
1683 a3x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
1691 subroutine c3cof(n, C3x)
1693 integer ord, nC3, nC3x
1694 parameter(ord = 6, nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2)
1699 double precision C3x(0:nc3x-1)
1701 integer o, m, l, j, k
1702 double precision polval,
1703 + coeff(((nc3-1)*(nc3**2 + 7*nc3 - 2*(nc3/2)))/8)
1723 do 20 l = 1, nc3 - 1
1724 do 10 j = nc3 - 1, l, -1
1725 m = min(nc3 - j - 1, j)
1726 c3x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
1735 subroutine c4cof(n, C4x)
1737 integer ord, nC4, nC4x
1738 parameter(ord = 6, nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
1743 double precision C4x(0:nc4x-1)
1745 integer o, m, l, j, k
1746 double precision polval, coeff((nc4 * (nc4 + 1) * (nc4 + 5)) / 6)
1748 + 97, 15015, 1088, 156, 45045, -224, -4784, 1573, 45045,
1749 + -10656, 14144, -4576, -858, 45045,
1750 + 64, 624, -4576, 6864, -3003, 15015,
1751 + 100, 208, 572, 3432, -12012, 30030, 45045,
1752 + 1, 9009, -2944, 468, 135135, 5792, 1040, -1287, 135135,
1753 + 5952, -11648, 9152, -2574, 135135,
1754 + -64, -624, 4576, -6864, 3003, 135135,
1755 + 8, 10725, 1856, -936, 225225, -8448, 4992, -1144, 225225,
1756 + -1440, 4160, -4576, 1716, 225225,
1757 + -136, 63063, 1024, -208, 105105,
1758 + 3584, -3328, 1144, 315315,
1759 + -128, 135135, -2560, 832, 405405, 128, 99099/
1763 do 20 l = 0, nc4 - 1
1764 do 10 j = nc4 - 1, l, -1
1766 c4x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
1775 double precision function sumx(u, v, t)
1777 double precision u, v
1781 double precision up, vpp
1792 double precision function angnm(x)
1796 if (x .ge. 180)
then 1798 else if (x .lt. -180)
then 1807 double precision function angnm2(x)
1811 double precision AngNm
1812 angnm2 = mod(x, 360d0)
1813 angnm2 = angnm(angnm2)
1818 double precision function angdif(x, y)
1824 double precision x, y
1826 double precision d, t, sumx
1828 if ((d - 180d0) + t .gt. 0d0)
then 1830 else if ((d + 180d0) + t .le. 0d0)
then 1838 double precision function angrnd(x)
1848 double precision y, z
1852 if (y .lt. z) y = z - (z - y)
1853 angrnd = 0 + sign(y, x)
1858 subroutine swap(x, y)
1860 double precision x, y
1870 double precision function hypotx(x, y)
1872 double precision x, y
1874 hypotx = sqrt(x**2 + y**2)
1879 subroutine norm2(sinx, cosx)
1881 double precision sinx, cosx
1883 double precision hypotx, r
1884 r = hypotx(sinx, cosx)
1891 double precision function log1px(x)
1895 double precision csmgt, y, z
1898 log1px = csmgt(x, x * log(y) / z, z .eq. 0)
1903 double precision function atanhx(x)
1907 double precision log1px, y
1909 y = log1px(2 * y/(1 - y))/2
1915 double precision function cbrt(x)
1919 cbrt = sign(abs(x)**(1/3d0), x)
1924 double precision function csmgt(x, y, p)
1926 double precision x, y
1938 double precision function trgsum(sinp, sinx, cosx, c, n)
1947 double precision sinx, cosx, c(n)
1949 double precision ar, y0, y1
1953 ar = 2 * (cosx - sinx) * (cosx + sinx)
1955 if (mod(n, 2) .eq. 1)
then 1966 y1 = ar * y0 - y1 + c(k)
1967 y0 = ar * y1 - y0 + c(k-1)
1971 trgsum = 2 * sinx * cosx * y0
1974 trgsum = cosx * (y0 - y1)
1980 integer function trnsit(lon1, lon2)
1982 double precision lon1, lon2
1984 double precision lon1x, lon2x, lon12, AngNm, AngDif
1987 lon12 = angdif(lon1x, lon2x)
1989 if (lon1x .lt. 0 .and. lon2x .ge. 0 .and. lon12 .gt. 0)
then 1991 else if (lon2x .lt. 0 .and. lon1x .ge. 0 .and. lon12 .lt. 0)
then 1998 subroutine accini(s)
2001 double precision s(2)
2009 subroutine accadd(s, y)
2014 double precision s(2)
2016 double precision z, u, sumx
2017 z = sumx(y, s(2), u)
2018 s(1) = sumx(z, s(1), s(2))
2019 if (s(1) .eq. 0)
then 2028 double precision function polval(N, p, x)
2031 double precision p(0:n), x
2040 polval = polval * x + p(i)