183 subroutine direct(a, f, lat1, lon1, azi1, s12a12, flags,
184 + lat2, lon2, azi2, omask, a12s12, m12, mm12, mm21, ss12)
186 double precision a, f, lat1, lon1, azi1, s12a12
189 double precision lat2, lon2, azi2
191 double precision a12s12, m12, MM12, MM21, SS12
193 integer ord, nC1, nC1p, nC2, nA3, nA3x, nC3, nC3x, nC4, nC4x
194 parameter(ord = 6, nc1 = ord, nc1p = ord,
195 + nc2 = ord, na3 = ord, na3x = na3,
196 + nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2,
197 + nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
198 double precision A3x(0:na3x-1), C3x(0:nc3x-1), C4x(0:nc4x-1),
199 + c1a(nc1), c1pa(nc1p), c2a(nc2), c3a(nc3-1), c4a(0:nc4-1)
201 double precision csmgt, atanhx, hypotx,
202 + angnm, angrnd, trgsum, a1m1f, a2m1f, a3f, atn2dx
203 logical arcmod, unroll, arcp, redlp, scalp, areap
204 double precision e2, f1, ep2, n, b, c2,
205 + salp0, calp0, k2, eps,
206 + salp1, calp1, ssig1, csig1, cbet1, sbet1, dn1, somg1, comg1,
207 + salp2, calp2, ssig2, csig2, sbet2, cbet2, dn2, somg2, comg2,
208 + ssig12, csig12, salp12, calp12, omg12, lam12, lon12,
209 + sig12, stau1, ctau1, tau12, t, s, c, serr, e,
210 + a1m1, a2m1, a3c, a4, ab1, ab2,
211 + b11, b12, b21, b22, b31, b41, b42, j12
213 double precision dblmin, dbleps, pi, degree, tiny,
214 + tol0, tol1, tol2, tolb, xthrsh
215 integer digits, maxit1, maxit2
217 common /geocom/ dblmin, dbleps, pi, degree, tiny,
218 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
220 if (.not.init)
call geoini
229 arcmod = mod(flags/1, 2) == 1
230 unroll = mod(flags/2, 2) == 1
232 arcp = mod(omask/1, 2) == 1
233 redlp = mod(omask/2, 2) == 1
234 scalp = mod(omask/4, 2) == 1
235 areap = mod(omask/8, 2) == 1
240 else if (e2 .gt. 0)
then
241 c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
243 c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
249 if (areap)
call c4cof(n, c4x)
252 call sncsdx(angrnd(azi1), salp1, calp1)
254 call sncsdx(angrnd(lat1), sbet1, cbet1)
256 call norm2(sbet1, cbet1)
258 cbet1 = max(tiny, cbet1)
259 dn1 = sqrt(1 + ep2 * sbet1**2)
263 salp0 = salp1 * cbet1
266 calp0 = hypotx(calp1, salp1 * sbet1)
277 somg1 = salp0 * sbet1
278 csig1 = csmgt(cbet1 * calp1, 1d0, sbet1 .ne. 0 .or. calp1 .ne. 0)
281 call norm2(ssig1, csig1)
285 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
289 b11 = trgsum(.true., ssig1, csig1, c1a, nc1)
293 stau1 = ssig1 * c + csig1 * s
294 ctau1 = csig1 * c - ssig1 * s
298 if (.not. arcmod)
call c1pf(eps, c1pa)
300 if (redlp .or. scalp)
then
303 b21 = trgsum(.true., ssig1, csig1, c2a, nc2)
310 call c3f(eps, c3x, c3a)
311 a3c = -f * salp0 * a3f(eps, a3x)
312 b31 = trgsum(.true., ssig1, csig1, c3a, nc3-1)
315 call c4f(eps, c4x, c4a)
317 a4 = a**2 * calp0 * salp0 * e2
318 b41 = trgsum(.false., ssig1, csig1, c4a, nc4)
327 sig12 = s12a12 * degree
328 call sncsdx(s12a12, ssig12, csig12)
333 tau12 = s12a12 / (b * (1 + a1m1))
337 b12 = - trgsum(.true.,
338 + stau1 * c + ctau1 * s, ctau1 * c - stau1 * s, c1pa, nc1p)
339 sig12 = tau12 - (b12 - b11)
342 if (abs(f) .gt. 0.01d0)
then
364 ssig2 = ssig1 * csig12 + csig1 * ssig12
365 csig2 = csig1 * csig12 - ssig1 * ssig12
366 b12 = trgsum(.true., ssig2, csig2, c1a, nc1)
367 serr = (1 + a1m1) * (sig12 + (b12 - b11)) - s12a12 / b
368 sig12 = sig12 - serr / sqrt(1 + k2 * ssig2**2)
376 ssig2 = ssig1 * csig12 + csig1 * ssig12
377 csig2 = csig1 * csig12 - ssig1 * ssig12
378 dn2 = sqrt(1 + k2 * ssig2**2)
379 if (arcmod .or. abs(f) .gt. 0.01d0)
380 + b12 = trgsum(.true., ssig2, csig2, c1a, nc1)
381 ab1 = (1 + a1m1) * (b12 - b11)
384 sbet2 = calp0 * ssig2
386 cbet2 = hypotx(salp0, calp0 * csig2)
387 if (cbet2 .eq. 0)
then
394 somg2 = salp0 * ssig2
399 calp2 = calp0 * csig2
405 + - (atan2( ssig2, csig2) - atan2( ssig1, csig1))
406 + + (atan2(e * somg2, comg2) - atan2(e * somg1, comg1)))
408 omg12 = atan2(somg2 * comg1 - comg2 * somg1,
409 + comg2 * comg1 + somg2 * somg1)
412 lam12 = omg12 + a3c *
413 + ( sig12 + (trgsum(.true., ssig2, csig2, c3a, nc3-1)
415 lon12 = lam12 / degree
419 lon2 = angnm(angnm(lon1) + angnm(lon12))
421 lat2 = atn2dx(sbet2, f1 * cbet2)
423 azi2 = atn2dx(salp2, calp2)
425 if (redlp .or. scalp)
then
426 b22 = trgsum(.true., ssig2, csig2, c2a, nc2)
427 ab2 = (1 + a2m1) * (b22 - b21)
428 j12 = (a1m1 - a2m1) * sig12 + (ab1 - ab2)
432 if (redlp) m12 = b * ((dn2 * (csig1 * ssig2) -
433 + dn1 * (ssig1 * csig2)) - csig1 * csig2 * j12)
435 t = k2 * (ssig2 - ssig1) * (ssig2 + ssig1) / (dn1 + dn2)
436 mm12 = csig12 + (t * ssig2 - csig2 * j12) * ssig1 / dn1
437 mm21 = csig12 - (t * ssig1 - csig1 * j12) * ssig2 / dn2
441 b42 = trgsum(.false., ssig2, csig2, c4a, nc4)
442 if (calp0 .eq. 0 .or. salp0 .eq. 0)
then
444 salp12 = salp2 * calp1 - calp2 * salp1
445 calp12 = calp2 * calp1 + salp2 * salp1
449 if (salp12 .eq. 0 .and. calp12 .lt. 0)
then
450 salp12 = tiny * calp1
462 salp12 = calp0 * salp0 *
463 + csmgt(csig1 * (1 - csig12) + ssig12 * ssig1,
464 + ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1),
466 calp12 = salp0**2 + calp0**2 * csig1 * csig2
468 ss12 = c2 * atan2(salp12, calp12) + a4 * (b42 - b41)
471 if (arcp) a12s12 = csmgt(b * ((1 + a1m1) * sig12 + ab1),
472 + sig12 / degree, arcmod)
522 subroutine invers(a, f, lat1, lon1, lat2, lon2,
523 + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
525 double precision a, f, lat1, lon1, lat2, lon2
528 double precision s12, azi1, azi2
530 double precision a12, m12, MM12, MM21, SS12
532 integer ord, nA3, nA3x, nC3, nC3x, nC4, nC4x, nC
533 parameter(ord = 6, na3 = ord, na3x = na3,
534 + nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2,
535 + nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2,
537 double precision A3x(0:na3x-1), C3x(0:nc3x-1), C4x(0:nc4x-1),
540 double precision csmgt, atanhx, hypotx,
541 + angdif, angrnd, trgsum, lam12f, invsta, atn2dx
542 integer latsgn, lonsgn, swapp, numit
543 logical arcp, redlp, scalp, areap, merid, tripn, tripb
545 double precision e2, f1, ep2, n, b, c2,
546 + lat1x, lat2x, salp0, calp0, k2, eps,
547 + salp1, calp1, ssig1, csig1, cbet1, sbet1, dbet1, dn1,
548 + salp2, calp2, ssig2, csig2, sbet2, cbet2, dbet2, dn2,
549 + slam12, clam12, salp12, calp12, omg12, lam12, lon12,
550 + salp1a, calp1a, salp1b, calp1b,
551 + dalp1, sdalp1, cdalp1, nsalp1, alp12, somg12, domg12,
552 + sig12, v, dv, dnm, dummy,
553 + a4, b41, b42, s12x, m12x, a12x
555 double precision dblmin, dbleps, pi, degree, tiny,
556 + tol0, tol1, tol2, tolb, xthrsh
557 integer digits, maxit1, maxit2, lmask
559 common /geocom/ dblmin, dbleps, pi, degree, tiny,
560 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
562 if (.not.init)
call geoini
571 arcp = mod(omask/1, 2) == 1
572 redlp = mod(omask/2, 2) == 1
573 scalp = mod(omask/4, 2) == 1
574 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)
599 lon12 = angrnd(angdif(lon1, lon2))
601 if (lon12 .ge. 0)
then
606 lon12 = lon12 * lonsgn
611 if (abs(lat1x) .ge. abs(lat2x))
then
616 if (swapp .lt. 0)
then
618 call swap(lat1x, lat2x)
621 if (lat1x .lt. 0)
then
626 lat1x = lat1x * latsgn
627 lat2x = lat2x * latsgn
640 call sncsdx(lat1x, sbet1, cbet1)
642 call norm2(sbet1, cbet1)
644 cbet1 = max(tiny, cbet1)
646 call sncsdx(lat2x, sbet2, cbet2)
648 call norm2(sbet2, cbet2)
650 cbet2 = max(tiny, cbet2)
661 if (cbet1 .lt. -sbet1)
then
662 if (cbet2 .eq. cbet1) sbet2 = sign(sbet1, sbet2)
664 if (abs(sbet2) .eq. -sbet1) cbet2 = cbet1
667 dn1 = sqrt(1 + ep2 * sbet1**2)
668 dn2 = sqrt(1 + ep2 * sbet2**2)
670 lam12 = lon12 * degree
671 call sncsdx(lon12, slam12, clam12)
675 merid = lat1x .eq. -90 .or. slam12 == 0
691 csig1 = calp1 * cbet1
693 csig2 = calp2 * cbet2
696 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, 0d0),
697 + csig1 * csig2 + ssig1 * ssig2)
698 call lengs(n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
699 + cbet1, cbet2, lmask,
700 + s12x, m12x, dummy, mm12, mm21, ep2, ca)
709 if (sig12 .lt. 1 .or. m12x .ge. 0)
then
710 if (sig12 .lt. 3 * tiny)
then
717 a12x = sig12 / degree
725 if (.not. merid .and. sbet1 .eq. 0 .and.
726 + (f .le. 0 .or. lam12 .le. pi - f * pi))
then
736 m12x = b * sin(sig12)
742 else if (.not. merid)
then
747 sig12 = invsta(sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12,
748 + f, a3x, salp1, calp1, salp2, calp2, dnm, ca)
750 if (sig12 .ge. 0)
then
752 s12x = sig12 * b * dnm
753 m12x = dnm**2 * b * sin(sig12 / dnm)
755 mm12 = cos(sig12 / dnm)
758 a12x = sig12 / degree
759 omg12 = lam12 / (f1 * dnm)
781 do 10 numit = 0, maxit2-1
784 v = lam12f(sbet1, cbet1, dn1, sbet2, cbet2, dn2,
785 + salp1, calp1, f, a3x, c3x, salp2, calp2, sig12,
786 + ssig1, csig1, ssig2, csig2,
787 + eps, omg12, numit .lt. maxit1, dv,
792 + .not. (abs(v) .ge. csmgt(8d0, 2d0, tripn) * tol0))
795 if (v .gt. 0 .and. (numit .gt. maxit1 .or.
796 + calp1/salp1 .gt. calp1b/salp1b))
then
799 else if (v .lt. 0 .and. (numit .gt. maxit1 .or.
800 + calp1/salp1 .lt. calp1a/salp1a))
then
804 if (numit .lt. maxit1 .and. dv .gt. 0)
then
808 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1
809 if (nsalp1 .gt. 0 .and. abs(dalp1) .lt. pi)
then
810 calp1 = calp1 * cdalp1 - salp1 * sdalp1
812 call norm2(salp1, calp1)
816 tripn = abs(v) .le. 16 * tol0
828 salp1 = (salp1a + salp1b)/2
829 calp1 = (calp1a + calp1b)/2
830 call norm2(salp1, calp1)
832 tripb = abs(salp1a - salp1) + (calp1a - calp1) .lt. tolb
833 + .or. abs(salp1 - salp1b) + (calp1 - calp1b) .lt. tolb
836 call lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
837 + cbet1, cbet2, lmask,
838 + s12x, m12x, dummy, mm12, mm21, ep2, ca)
841 a12x = sig12 / degree
842 omg12 = lam12 - omg12
848 if (redlp) m12 = 0 + m12x
852 salp0 = salp1 * cbet1
853 calp0 = hypotx(calp1, salp1 * sbet1)
854 if (calp0 .ne. 0 .and. salp0 .ne. 0)
then
857 csig1 = calp1 * cbet1
859 csig2 = calp2 * cbet2
861 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
863 a4 = a**2 * calp0 * salp0 * e2
864 call norm2(ssig1, csig1)
865 call norm2(ssig2, csig2)
866 call c4f(eps, c4x, ca)
867 b41 = trgsum(.false., ssig1, csig1, ca, nc4)
868 b42 = trgsum(.false., ssig2, csig2, ca, nc4)
869 ss12 = a4 * (b42 - b41)
875 if (.not. merid .and. omg12 .lt. 0.75d0 * pi
876 + .and. sbet2 - sbet1 .lt. 1.75d0)
then
881 domg12 = 1 + cos(omg12)
884 alp12 = 2 * atan2(somg12 * (sbet1 * dbet2 + sbet2 * dbet1),
885 + domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) )
888 salp12 = salp2 * calp1 - calp2 * salp1
889 calp12 = calp2 * calp1 + salp2 * salp1
894 if (salp12 .eq. 0 .and. calp12 .lt. 0)
then
895 salp12 = tiny * calp1
898 alp12 = atan2(salp12, calp12)
900 ss12 = ss12 + c2 * alp12
901 ss12 = ss12 * swapp * lonsgn * latsgn
907 if (swapp .lt. 0)
then
908 call swap(salp1, salp2)
909 call swap(calp1, calp2)
910 if (scalp)
call swap(mm12, mm21)
913 salp1 = salp1 * swapp * lonsgn
914 calp1 = calp1 * swapp * latsgn
915 salp2 = salp2 * swapp * lonsgn
916 calp2 = calp2 * swapp * latsgn
919 azi1 = atn2dx(salp1, calp1)
920 azi2 = atn2dx(salp2, calp2)
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)
1011 subroutine geover(major, minor, patch)
1013 integer major, minor, patch
1025 double precision dblmin, dbleps, pi, degree, tiny,
1026 + tol0, tol1, tol2, tolb, xthrsh
1027 integer digits, maxit1, maxit2
1030 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1031 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1035 double precision dblmin, dbleps, pi, degree, tiny,
1036 + tol0, tol1, tol2, tolb, xthrsh
1037 integer digits, maxit1, maxit2
1039 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1040 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1043 dblmin = 0.5d0**1022
1044 dbleps = 0.5d0**(digits-1)
1046 pi = atan2(0d0, -1d0)
1057 xthrsh = 1000 * tol2
1059 maxit2 = maxit1 + digits + 10
1066 subroutine lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1067 + cbet1, cbet2, omask,
1068 + s12b, m12b, m0, mm12, mm21, ep2, ca)
1070 double precision eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1074 double precision s12b, m12b, m0, MM12, MM21
1076 double precision Ca(*)
1078 integer ord, nC1, nC2
1079 parameter(ord = 6, nc1 = ord, nc2 = ord)
1081 double precision A1m1f, A2m1f, TrgSum
1082 double precision m0x, J12, A1, A2, B1, B2, csig12, t, Cb(nc2)
1083 logical distp, redlp, scalp
1089 distp = mod(omask/16, 2) == 1
1090 redlp = mod(omask/2, 2) == 1
1091 scalp = mod(omask/4, 2) == 1
1098 if (distp .or. redlp .or. scalp)
then
1101 if (redlp .or. scalp)
then
1110 b1 = trgsum(.true., ssig2, csig2, ca, nc1) -
1111 + trgsum(.true., ssig1, csig1, ca, nc1)
1113 s12b = a1 * (sig12 + b1)
1114 if (redlp .or. scalp)
then
1115 b2 = trgsum(.true., ssig2, csig2, cb, nc2) -
1116 + trgsum(.true., ssig1, csig1, cb, nc2)
1117 j12 = m0x * sig12 + (a1 * b1 - a2 * b2)
1119 else if (redlp .or. scalp)
then
1122 cb(l) = a1 * ca(l) - a2 * cb(l)
1124 j12 = m0x * sig12 + (trgsum(.true., ssig2, csig2, cb, nc2) -
1125 + trgsum(.true., ssig1, csig1, cb, nc2))
1132 m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
1133 + csig1 * csig2 * j12
1136 csig12 = csig1 * csig2 + ssig1 * ssig2
1137 t = ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2)
1138 mm12 = csig12 + (t * ssig2 - csig2 * j12) * ssig1 / dn1
1139 mm21 = csig12 - (t * ssig1 - csig1 * j12) * ssig2 / dn2
1145 double precision function astrd(x, y)
1149 double precision x, y
1151 double precision cbrt, csmgt
1152 double precision k, p, q, r, S, r2, r3, disc, u,
1153 + t3, t, ang, v, uv, w
1158 if ( .not. (q .eq. 0 .and. r .lt. 0) )
then
1167 disc = s * (s + 2 * r3)
1169 if (disc .ge. 0)
then
1175 t3 = t3 + csmgt(-sqrt(disc), sqrt(disc), t3 .lt. 0)
1180 if (t .ne. 0) u = u + t + r2 / t
1183 ang = atan2(sqrt(-disc), -(s + r3))
1186 u = u + 2 * r * cos(ang / 3)
1192 uv = csmgt(q / (v - u), u + v, u .lt. 0)
1194 w = (uv - q) / (2 * v)
1198 k = uv / (sqrt(uv + w**2) + w)
1210 double precision function invsta(sbet1, cbet1, dn1,
1211 + sbet2, cbet2, dn2, lam12, f, a3x,
1212 + salp1, calp1, salp2, calp2, dnm,
1218 double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12,
1221 double precision salp1, calp1, salp2, calp2, dnm
1223 double precision Ca(*)
1225 double precision csmgt, hypotx, A3f, Astrd
1227 double precision f1, e2, ep2, n, etol2, k2, eps, sig12,
1228 + sbet12, cbet12, sbt12a, omg12, somg12, comg12, ssig12, csig12,
1229 + x, y, lamscl, betscl, cbt12a, bt12a, m12b, m0, dummy,
1232 double precision dblmin, dbleps, pi, degree, tiny,
1233 + tol0, tol1, tol2, tolb, xthrsh
1234 integer digits, maxit1, maxit2
1236 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1237 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1253 etol2 = 0.1d0 * tol2 /
1254 + sqrt( max(0.001d0, abs(f)) * min(1d0, 1 - f/2) / 2 )
1259 sbet12 = sbet2 * cbet1 - cbet2 * sbet1
1260 cbet12 = cbet2 * cbet1 + sbet2 * sbet1
1261 sbt12a = sbet2 * cbet1 + cbet2 * sbet1
1263 shortp = cbet12 .ge. 0 .and. sbet12 .lt. 0.5d0 .and.
1264 + cbet2 * lam12 .lt. 0.5d0
1268 sbetm2 = (sbet1 + sbet2)**2
1271 sbetm2 = sbetm2 / (sbetm2 + (cbet1 + cbet2)**2)
1272 dnm = sqrt(1 + ep2 * sbetm2)
1273 omg12 = omg12 / (f1 * dnm)
1278 salp1 = cbet2 * somg12
1279 calp1 = csmgt(sbet12 + cbet2 * sbet1 * somg12**2 / (1 + comg12),
1280 + sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12),
1283 ssig12 = hypotx(salp1, calp1)
1284 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12
1286 if (shortp .and. ssig12 .lt. etol2)
then
1288 salp2 = cbet1 * somg12
1289 calp2 = sbet12 - cbet1 * sbet2 *
1290 + csmgt(somg12**2 / (1 + comg12), 1 - comg12, comg12 .ge. 0)
1291 call norm2(salp2, calp2)
1293 sig12 = atan2(ssig12, csig12)
1294 else if (abs(n) .gt. 0.1d0 .or. csig12 .ge. 0 .or.
1295 + ssig12 .ge. 6 * abs(n) * pi * cbet1**2)
then
1304 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
1305 lamscl = f * cbet1 * a3f(eps, a3x) * pi
1306 betscl = lamscl * cbet1
1307 x = (lam12 - pi) / lamscl
1311 cbt12a = cbet2 * cbet1 - sbet2 * sbet1
1312 bt12a = atan2(sbt12a, cbt12a)
1315 call lengs(n, pi + bt12a,
1316 + sbet1, -cbet1, dn1, sbet2, cbet2, dn2, cbet1, cbet2, 2,
1317 + dummy, m12b, m0, dummy, dummy, ep2, ca)
1318 x = -1 + m12b / (cbet1 * cbet2 * m0 * pi)
1319 betscl = csmgt(sbt12a / x, -f * cbet1**2 * pi,
1321 lamscl = betscl / cbet1
1322 y = (lam12 - pi) / lamscl
1325 if (y .gt. -tol1 .and. x .gt. -1 - xthrsh)
then
1328 salp1 = min(1d0, -x)
1329 calp1 = - sqrt(1 - salp1**2)
1331 calp1 = max(csmgt(0d0, 1d0, x .gt. -tol1), x)
1332 salp1 = sqrt(1 - calp1**2)
1371 + csmgt(-x * k/(1 + k), -y * (1 + k)/k, f .ge. 0)
1372 somg12 = sin(omg12a)
1373 comg12 = -cos(omg12a)
1375 salp1 = cbet2 * somg12
1376 calp1 = sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12)
1380 if (.not. (salp1 .le. 0))
then
1381 call norm2(salp1, calp1)
1391 double precision function lam12f(sbet1, cbet1, dn1,
1392 + sbet2, cbet2, dn2, salp1, calp1, f, a3x, c3x, salp2, calp2,
1393 + sig12, ssig1, csig1, ssig2, csig2, eps, domg12, diffp, dlam12,
1396 double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2,
1397 + salp1, calp1, f, a3x(*), c3x(*)
1400 double precision salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
1403 double precision dlam12
1405 double precision Ca(*)
1408 parameter(ord = 6, nc3 = ord)
1410 double precision csmgt, hypotx, A3f, TrgSum
1412 double precision f1, e2, ep2, salp0, calp0,
1413 + somg1, comg1, somg2, comg2, omg12, lam12, b312, h0, k2, dummy
1415 double precision dblmin, dbleps, pi, degree, tiny,
1416 + tol0, tol1, tol2, tolb, xthrsh
1417 integer digits, maxit1, maxit2
1419 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1420 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1427 if (sbet1 .eq. 0 .and. calp1 .eq. 0) calp1 = -tiny
1430 salp0 = salp1 * cbet1
1432 calp0 = hypotx(calp1, salp1 * sbet1)
1437 somg1 = salp0 * sbet1
1438 csig1 = calp1 * cbet1
1440 call norm2(ssig1, csig1)
1447 salp2 = csmgt(salp0 / cbet2, salp1, cbet2 .ne. cbet1)
1452 calp2 = csmgt(sqrt((calp1 * cbet1)**2 +
1453 + csmgt((cbet2 - cbet1) * (cbet1 + cbet2),
1454 + (sbet1 - sbet2) * (sbet1 + sbet2),
1455 + cbet1 .lt. -sbet1)) / cbet2,
1456 + abs(calp1), cbet2 .ne. cbet1 .or. abs(sbet2) .ne. -sbet1)
1460 somg2 = salp0 * sbet2
1461 csig2 = calp2 * cbet2
1463 call norm2(ssig2, csig2)
1467 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, 0d0),
1468 + csig1 * csig2 + ssig1 * ssig2)
1471 omg12 = atan2(max(comg1 * somg2 - somg1 * comg2, 0d0),
1472 + comg1 * comg2 + somg1 * somg2)
1474 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
1475 call c3f(eps, c3x, ca)
1476 b312 = (trgsum(.true., ssig2, csig2, ca, nc3-1) -
1477 + trgsum(.true., ssig1, csig1, ca, nc3-1))
1478 h0 = -f * a3f(eps, a3x)
1479 domg12 = salp0 * h0 * (sig12 + b312)
1480 lam12 = omg12 + domg12
1483 if (calp2 .eq. 0)
then
1484 dlam12 = - 2 * f1 * dn1 / sbet1
1486 call lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1488 + dummy, dlam12, dummy, dummy, dummy, ep2, ca)
1489 dlam12 = dlam12 * f1 / (calp2 * cbet2)
1497 double precision function a3f(eps, A3x)
1499 integer ord, nA3, nA3x
1500 parameter(ord = 6, na3 = ord, na3x = na3)
1503 double precision eps
1505 double precision A3x(0: na3x-1)
1507 double precision polval
1508 a3f = polval(na3 - 1, a3x, eps)
1513 subroutine c3f(eps, C3x, c)
1516 integer ord, nC3, nC3x
1517 parameter(ord = 6, nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2)
1520 double precision eps, C3x(0:nc3x-1)
1522 double precision c(nc3-1)
1525 double precision mult, polval
1529 do 10 l = 1, nc3 - 1
1532 c(l) = mult * polval(m, c3x(o), eps)
1539 subroutine c4f(eps, C4x, c)
1542 integer ord, nC4, nC4x
1543 parameter(ord = 6, nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
1546 double precision eps, C4x(0:nc4x-1)
1548 double precision c(0:nc4-1)
1551 double precision mult, polval
1555 do 10 l = 0, nc4 - 1
1557 c(l) = mult * polval(m, c4x(o), eps)
1565 double precision function a1m1f(eps)
1568 double precision eps
1571 integer ord, nA1, o, m
1572 parameter(ord = 6, na1 = ord)
1573 double precision polval, coeff(na1/2 + 2)
1574 data coeff /1, 4, 64, 0, 256/
1578 t = polval(m, coeff(o), eps**2) / coeff(o + m + 1)
1579 a1m1f = (t + eps) / (1 - eps)
1584 subroutine c1f(eps, c)
1587 parameter(ord = 6, nc1 = ord)
1590 double precision eps
1592 double precision c(nc1)
1594 double precision eps2, d
1596 double precision polval, coeff((nc1**2 + 7*nc1 - 2*(nc1/2))/4)
1599 + -9, 64, -128, 2048,
1610 c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
1618 subroutine c1pf(eps, c)
1621 parameter(ord = 6, nc1p = ord)
1624 double precision eps
1626 double precision c(nc1p)
1628 double precision eps2, d
1630 double precision polval, coeff((nc1p**2 + 7*nc1p - 2*(nc1p/2))/4)
1632 + 205, -432, 768, 1536,
1633 + 4005, -4736, 3840, 12288,
1635 + -7173, 2695, 7680,
1644 c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
1653 double precision function a2m1f(eps)
1655 double precision eps
1658 integer ord, nA2, o, m
1659 parameter(ord = 6, na2 = ord)
1660 double precision polval, coeff(na2/2 + 2)
1661 data coeff /-11, -28, -192, 0, 256/
1665 t = polval(m, coeff(o), eps**2) / coeff(o + m + 1)
1666 a2m1f = (t - eps) / (1 + eps)
1671 subroutine c2f(eps, c)
1674 parameter(ord = 6, nc2 = ord)
1677 double precision eps
1679 double precision c(nc2)
1681 double precision eps2, d
1683 double precision polval, coeff((nc2**2 + 7*nc2 - 2*(nc2/2))/4)
1686 + 35, 64, 384, 2048,
1697 c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
1705 subroutine a3cof(n, A3x)
1707 integer ord, nA3, nA3x
1708 parameter(ord = 6, na3 = ord, na3x = na3)
1713 double precision A3x(0:na3x-1)
1716 double precision polval, coeff((na3**2 + 7*na3 - 2*(na3/2))/4)
1727 do 10 j = na3 - 1, 0, -1
1728 m = min(na3 - j - 1, j)
1729 a3x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
1737 subroutine c3cof(n, C3x)
1739 integer ord, nC3, nC3x
1740 parameter(ord = 6, nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2)
1745 double precision C3x(0:nc3x-1)
1747 integer o, m, l, j, k
1748 double precision polval,
1749 + coeff(((nc3-1)*(nc3**2 + 7*nc3 - 2*(nc3/2)))/8)
1769 do 20 l = 1, nc3 - 1
1770 do 10 j = nc3 - 1, l, -1
1771 m = min(nc3 - j - 1, j)
1772 c3x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
1781 subroutine c4cof(n, C4x)
1783 integer ord, nC4, nC4x
1784 parameter(ord = 6, nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
1789 double precision C4x(0:nc4x-1)
1791 integer o, m, l, j, k
1792 double precision polval, coeff((nc4 * (nc4 + 1) * (nc4 + 5)) / 6)
1794 + 97, 15015, 1088, 156, 45045, -224, -4784, 1573, 45045,
1795 + -10656, 14144, -4576, -858, 45045,
1796 + 64, 624, -4576, 6864, -3003, 15015,
1797 + 100, 208, 572, 3432, -12012, 30030, 45045,
1798 + 1, 9009, -2944, 468, 135135, 5792, 1040, -1287, 135135,
1799 + 5952, -11648, 9152, -2574, 135135,
1800 + -64, -624, 4576, -6864, 3003, 135135,
1801 + 8, 10725, 1856, -936, 225225, -8448, 4992, -1144, 225225,
1802 + -1440, 4160, -4576, 1716, 225225,
1803 + -136, 63063, 1024, -208, 105105,
1804 + 3584, -3328, 1144, 315315,
1805 + -128, 135135, -2560, 832, 405405, 128, 99099/
1809 do 20 l = 0, nc4 - 1
1810 do 10 j = nc4 - 1, l, -1
1812 c4x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
1821 double precision function sumx(u, v, t)
1823 double precision u, v
1827 double precision up, vpp
1838 double precision function angnm(x)
1842 angnm = mod(x, 360d0)
1843 if (angnm .lt. -180)
then
1845 else if (angnm .ge. 180)
then
1852 double precision function angdif(x, y)
1858 double precision x, y
1860 double precision d, t, sumx, AngNm
1861 d = - angnm(sumx(angnm(x), angnm(-y), t))
1862 if (d .eq. 180 .and. t .lt. 0)
then
1870 double precision function angrnd(x)
1880 double precision y, z
1884 if (y .lt. z) y = z - (z - y)
1885 angrnd = 0 + sign(y, x)
1890 subroutine swap(x, y)
1892 double precision x, y
1902 double precision function hypotx(x, y)
1904 double precision x, y
1906 hypotx = sqrt(x**2 + y**2)
1911 subroutine norm2(x, y)
1913 double precision x, y
1915 double precision hypotx, r
1923 double precision function log1px(x)
1927 double precision csmgt, y, z
1930 log1px = csmgt(x, x * log(y) / z, z .eq. 0)
1935 double precision function atanhx(x)
1939 double precision log1px, y
1941 y = log1px(2 * y/(1 - y))/2
1947 double precision function cbrt(x)
1951 cbrt = sign(abs(x)**(1/3d0), x)
1956 double precision function csmgt(x, y, p)
1958 double precision x, y
1970 double precision function trgsum(sinp, sinx, cosx, c, n)
1979 double precision sinx, cosx, c(n)
1981 double precision ar, y0, y1
1985 ar = 2 * (cosx - sinx) * (cosx + sinx)
1987 if (mod(n, 2) .eq. 1)
then
1998 y1 = ar * y0 - y1 + c(k)
1999 y0 = ar * y1 - y0 + c(k-1)
2003 trgsum = 2 * sinx * cosx * y0
2006 trgsum = cosx * (y0 - y1)
2012 integer function trnsit(lon1, lon2)
2014 double precision lon1, lon2
2016 double precision lon1x, lon2x, lon12, AngNm, AngDif
2019 lon12 = angdif(lon1x, lon2x)
2021 if (lon1x .lt. 0 .and. lon2x .ge. 0 .and. lon12 .gt. 0)
then
2023 else if (lon2x .lt. 0 .and. lon1x .ge. 0 .and. lon12 .lt. 0)
then
2030 subroutine accini(s)
2033 double precision s(2)
2041 subroutine accadd(s, y)
2046 double precision s(2)
2048 double precision z, u, sumx
2049 z = sumx(y, s(2), u)
2050 s(1) = sumx(z, s(1), s(2))
2051 if (s(1) .eq. 0)
then
2060 subroutine sncsdx(x, sinx, cosx)
2065 double precision sinx, cosx
2067 double precision dblmin, dbleps, pi, degree, tiny,
2068 + tol0, tol1, tol2, tolb, xthrsh
2069 integer digits, maxit1, maxit2
2071 common /geocom/ dblmin, dbleps, pi, degree, tiny,
2072 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
2074 double precision r, s, c
2078 r = (r - 90 * q) * degree
2085 else if (q .eq. 1)
then
2088 else if (q .eq. 2)
then
2100 double precision function atn2dx(y, x)
2102 double precision x, y
2104 double precision dblmin, dbleps, pi, degree, tiny,
2105 + tol0, tol1, tol2, tolb, xthrsh
2106 integer digits, maxit1, maxit2
2108 common /geocom/ dblmin, dbleps, pi, degree, tiny,
2109 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
2111 double precision xx, yy
2113 if (abs(y) .gt. abs(x))
then
2126 atn2dx = atan2(yy, xx) / degree
2129 atn2dx = 180 - atn2dx
2131 atn2dx = -180 - atn2dx
2133 else if (q .eq. 2)
then
2134 atn2dx = 90 - atn2dx
2135 else if (q .eq. 3)
then
2136 atn2dx = -90 + atn2dx
2142 double precision function polval(N, p, x)
2145 double precision p(0:n), x
2154 polval = polval * x + p(i)