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 atanhx, hypotx,
203 + angnm, angrnd, trgsum, a1m1f, a2m1f, a3f, atn2dx, latfix
204 logical arcmod, unroll, arcp, redlp, scalp, areap
205 double precision e2, f1, ep2, n, b, c2,
206 + 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, 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) .eq. 1
231 unroll = mod(flags/2, 2) .eq. 1
233 arcp = mod(omask/1, 2) .eq. 1
234 redlp = mod(omask/2, 2) .eq. 1
235 scalp = mod(omask/4, 2) .eq. 1
236 areap = mod(omask/8, 2) .eq. 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 call sncsdx(angrnd(azi1), salp1, calp1)
255 call sncsdx(angrnd(latfix(lat1)), sbet1, cbet1)
257 call norm2x(sbet1, cbet1)
259 cbet1 = max(tiny, cbet1)
260 dn1 = sqrt(1 + ep2 * sbet1**2)
264 salp0 = salp1 * cbet1
267 calp0 = hypotx(calp1, salp1 * sbet1)
278 somg1 = salp0 * sbet1
279 if (sbet1 .ne. 0 .or. calp1 .ne. 0)
then
280 csig1 = cbet1 * calp1
286 call norm2x(ssig1, csig1)
290 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
294 b11 = trgsum(.true., ssig1, csig1, c1a, nc1)
298 stau1 = ssig1 * c + csig1 * s
299 ctau1 = csig1 * c - ssig1 * s
303 if (.not. arcmod)
call c1pf(eps, c1pa)
305 if (redlp .or. scalp)
then
308 b21 = trgsum(.true., ssig1, csig1, c2a, nc2)
315 call c3f(eps, c3x, c3a)
316 a3c = -f * salp0 * a3f(eps, a3x)
317 b31 = trgsum(.true., ssig1, csig1, c3a, nc3-1)
320 call c4f(eps, c4x, c4a)
322 a4 = a**2 * calp0 * salp0 * e2
323 b41 = trgsum(.false., ssig1, csig1, c4a, nc4)
332 sig12 = s12a12 * degree
333 call sncsdx(s12a12, ssig12, csig12)
338 tau12 = s12a12 / (b * (1 + a1m1))
342 b12 = - trgsum(.true.,
343 + stau1 * c + ctau1 * s, ctau1 * c - stau1 * s, c1pa, nc1p)
344 sig12 = tau12 - (b12 - b11)
347 if (abs(f) .gt. 0.01d0)
then
369 ssig2 = ssig1 * csig12 + csig1 * ssig12
370 csig2 = csig1 * csig12 - ssig1 * ssig12
371 b12 = trgsum(.true., ssig2, csig2, c1a, nc1)
372 serr = (1 + a1m1) * (sig12 + (b12 - b11)) - s12a12 / b
373 sig12 = sig12 - serr / sqrt(1 + k2 * ssig2**2)
381 ssig2 = ssig1 * csig12 + csig1 * ssig12
382 csig2 = csig1 * csig12 - ssig1 * ssig12
383 dn2 = sqrt(1 + k2 * ssig2**2)
384 if (arcmod .or. abs(f) .gt. 0.01d0)
385 + b12 = trgsum(.true., ssig2, csig2, c1a, nc1)
386 ab1 = (1 + a1m1) * (b12 - b11)
389 sbet2 = calp0 * ssig2
391 cbet2 = hypotx(salp0, calp0 * csig2)
392 if (cbet2 .eq. 0)
then
399 somg2 = salp0 * ssig2
404 calp2 = calp0 * csig2
410 + - (atan2( ssig2, csig2) - atan2( ssig1, csig1))
411 + + (atan2(e * somg2, comg2) - atan2(e * somg1, comg1)))
413 omg12 = atan2(somg2 * comg1 - comg2 * somg1,
414 + comg2 * comg1 + somg2 * somg1)
417 lam12 = omg12 + a3c *
418 + ( sig12 + (trgsum(.true., ssig2, csig2, c3a, nc3-1)
420 lon12 = lam12 / degree
424 lon2 = angnm(angnm(lon1) + angnm(lon12))
426 lat2 = atn2dx(sbet2, f1 * cbet2)
427 azi2 = atn2dx(salp2, calp2)
429 if (redlp .or. scalp)
then
430 b22 = trgsum(.true., ssig2, csig2, c2a, nc2)
431 ab2 = (1 + a2m1) * (b22 - b21)
432 j12 = (a1m1 - a2m1) * sig12 + (ab1 - ab2)
436 if (redlp) m12 = b * ((dn2 * (csig1 * ssig2) -
437 + dn1 * (ssig1 * csig2)) - csig1 * csig2 * j12)
439 t = k2 * (ssig2 - ssig1) * (ssig2 + ssig1) / (dn1 + dn2)
440 mm12 = csig12 + (t * ssig2 - csig2 * j12) * ssig1 / dn1
441 mm21 = csig12 - (t * ssig1 - csig1 * j12) * ssig2 / dn2
445 b42 = trgsum(.false., ssig2, csig2, c4a, nc4)
446 if (calp0 .eq. 0 .or. salp0 .eq. 0)
then
448 salp12 = salp2 * calp1 - calp2 * salp1
449 calp12 = calp2 * calp1 + salp2 * salp1
459 if (csig12 .le. 0)
then
460 salp12 = csig1 * (1 - csig12) + ssig12 * ssig1
462 salp12 = ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1)
464 salp12 = calp0 * salp0 * salp12
465 calp12 = salp0**2 + calp0**2 * csig1 * csig2
467 ss12 = c2 * atan2(salp12, calp12) + a4 * (b42 - b41)
472 a12s12 = b * ((1 + a1m1) * sig12 + ab1)
474 a12s12 = sig12 / degree
526 subroutine invers(a, f, lat1, lon1, lat2, lon2,
527 + s12, azi1, azi2, omask, a12, m12, MM12, MM21, SS12)
529 double precision a, f, lat1, lon1, lat2, lon2
532 double precision s12, azi1, azi2
534 double precision a12, m12, MM12, MM21, SS12
536 integer ord, nA3, nA3x, nC3, nC3x, nC4, nC4x, nC
537 parameter (ord = 6, na3 = ord, na3x = na3,
538 + nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2,
539 + nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2,
541 double precision A3x(0:nA3x-1), C3x(0:nC3x-1), C4x(0:nC4x-1),
544 double precision atanhx, hypotx,
545 + AngDif, AngRnd, TrgSum, Lam12f, InvSta, atn2dx, LatFix
546 integer latsgn, lonsgn, swapp, numit
547 logical arcp, redlp, scalp, areap, merid, tripn, tripb
549 double precision e2, f1, ep2, n, b, c2,
550 + lat1x, lat2x, salp0, calp0, k2, eps,
551 + salp1, calp1, ssig1, csig1, cbet1, sbet1, dbet1, dn1,
552 + salp2, calp2, ssig2, csig2, sbet2, cbet2, dbet2, dn2,
553 + slam12, clam12, salp12, calp12, omg12, lam12, lon12, lon12s,
554 + salp1a, calp1a, salp1b, calp1b,
555 + dalp1, sdalp1, cdalp1, nsalp1, alp12, somg12, comg12, domg12,
556 + sig12, v, dv, dnm, dummy,
557 + a4, b41, b42, s12x, m12x, a12x, sdomg12, cdomg12
559 double precision dblmin, dbleps, pi, degree, tiny,
560 + tol0, tol1, tol2, tolb, xthrsh
561 integer digits, maxit1, maxit2, lmask
563 common /geocom/ dblmin, dbleps, pi, degree, tiny,
564 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
566 if (.not.init)
call geoini
575 arcp = mod(omask/1, 2) .eq. 1
576 redlp = mod(omask/2, 2) .eq. 1
577 scalp = mod(omask/4, 2) .eq. 1
578 areap = mod(omask/8, 2) .eq. 1
588 else if (e2 .gt. 0)
then
589 c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
591 c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
597 if (areap)
call c4cof(n, c4x)
603 lon12 = angdif(lon1, lon2, lon12s)
605 if (lon12 .ge. 0)
then
610 lon12 = lonsgn * angrnd(lon12)
611 lon12s = angrnd((180 - lon12) - lonsgn * lon12s)
612 lam12 = lon12 * degree
613 if (lon12 .gt. 90)
then
614 call sncsdx(lon12s, slam12, clam12)
617 call sncsdx(lon12, slam12, clam12)
621 lat1x = angrnd(latfix(lat1))
622 lat2x = angrnd(latfix(lat2))
625 if (abs(lat1x) .lt. abs(lat2x))
then
630 if (swapp .lt. 0)
then
632 call swap(lat1x, lat2x)
635 if (lat1x .lt. 0)
then
640 lat1x = lat1x * latsgn
641 lat2x = lat2x * latsgn
654 call sncsdx(lat1x, sbet1, cbet1)
656 call norm2x(sbet1, cbet1)
658 cbet1 = max(tiny, cbet1)
660 call sncsdx(lat2x, sbet2, cbet2)
662 call norm2x(sbet2, cbet2)
664 cbet2 = max(tiny, cbet2)
675 if (cbet1 .lt. -sbet1)
then
676 if (cbet2 .eq. cbet1) sbet2 = sign(sbet1, sbet2)
678 if (abs(sbet2) .eq. -sbet1) cbet2 = cbet1
681 dn1 = sqrt(1 + ep2 * sbet1**2)
682 dn2 = sqrt(1 + ep2 * sbet2**2)
686 merid = lat1x .eq. -90 .or. slam12 .eq. 0
702 csig1 = calp1 * cbet1
704 csig2 = calp2 * cbet2
707 sig12 = atan2(0d0 + max(0d0, csig1 * ssig2 - ssig1 * csig2),
708 + csig1 * csig2 + ssig1 * ssig2)
709 call lengs(n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
710 + cbet1, cbet2, lmask,
711 + s12x, m12x, dummy, mm12, mm21, ep2, ca)
720 if (sig12 .lt. 1 .or. m12x .ge. 0)
then
721 if (sig12 .lt. 3 * tiny)
then
728 a12x = sig12 / degree
739 if (.not. merid .and. sbet1 .eq. 0 .and.
740 + (f .le. 0 .or. lon12s .ge. f * 180))
then
750 m12x = b * sin(sig12)
756 else if (.not. merid)
then
761 sig12 = invsta(sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12,
762 + slam12, clam12, f, a3x, salp1, calp1, salp2, calp2, dnm, ca)
764 if (sig12 .ge. 0)
then
766 s12x = sig12 * b * dnm
767 m12x = dnm**2 * b * sin(sig12 / dnm)
769 mm12 = cos(sig12 / dnm)
772 a12x = sig12 / degree
773 omg12 = lam12 / (f1 * dnm)
795 do 10 numit = 0, maxit2-1
798 v = lam12f(sbet1, cbet1, dn1, sbet2, cbet2, dn2,
799 + salp1, calp1, slam12, clam12, f, a3x, c3x, salp2, calp2,
800 + sig12, ssig1, csig1, ssig2, csig2,
801 + eps, domg12, numit .lt. maxit1, dv, ca)
808 if (tripb .or. .not. (abs(v) .ge. dummy * tol0))
811 if (v .gt. 0 .and. (numit .gt. maxit1 .or.
812 + calp1/salp1 .gt. calp1b/salp1b))
then
815 else if (v .lt. 0 .and. (numit .gt. maxit1 .or.
816 + calp1/salp1 .lt. calp1a/salp1a))
then
820 if (numit .lt. maxit1 .and. dv .gt. 0)
then
824 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1
825 if (nsalp1 .gt. 0 .and. abs(dalp1) .lt. pi)
then
826 calp1 = calp1 * cdalp1 - salp1 * sdalp1
828 call norm2x(salp1, calp1)
832 tripn = abs(v) .le. 16 * tol0
844 salp1 = (salp1a + salp1b)/2
845 calp1 = (calp1a + calp1b)/2
846 call norm2x(salp1, calp1)
848 tripb = abs(salp1a - salp1) + (calp1a - calp1) .lt. tolb
849 + .or. abs(salp1 - salp1b) + (calp1 - calp1b) .lt. tolb
852 call lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
853 + cbet1, cbet2, lmask,
854 + s12x, m12x, dummy, mm12, mm21, ep2, ca)
857 a12x = sig12 / degree
859 sdomg12 = sin(domg12)
860 cdomg12 = cos(domg12)
861 somg12 = slam12 * cdomg12 - clam12 * sdomg12
862 comg12 = clam12 * cdomg12 + slam12 * sdomg12
869 if (redlp) m12 = 0 + m12x
873 salp0 = salp1 * cbet1
874 calp0 = hypotx(calp1, salp1 * sbet1)
875 if (calp0 .ne. 0 .and. salp0 .ne. 0)
then
878 csig1 = calp1 * cbet1
880 csig2 = calp2 * cbet2
882 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
884 a4 = a**2 * calp0 * salp0 * e2
885 call norm2x(ssig1, csig1)
886 call norm2x(ssig2, csig2)
887 call c4f(eps, c4x, ca)
888 b41 = trgsum(.false., ssig1, csig1, ca, nc4)
889 b42 = trgsum(.false., ssig2, csig2, ca, nc4)
890 ss12 = a4 * (b42 - b41)
896 if (.not. merid .and. somg12 .gt. 1)
then
901 if (.not. merid .and. comg12 .ge. 0.7071d0
902 + .and. sbet2 - sbet1 .lt. 1.75d0)
then
909 alp12 = 2 * atan2(somg12 * (sbet1 * dbet2 + sbet2 * dbet1),
910 + domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) )
913 salp12 = salp2 * calp1 - calp2 * salp1
914 calp12 = calp2 * calp1 + salp2 * salp1
919 if (salp12 .eq. 0 .and. calp12 .lt. 0)
then
920 salp12 = tiny * calp1
923 alp12 = atan2(salp12, calp12)
925 ss12 = ss12 + c2 * alp12
926 ss12 = ss12 * swapp * lonsgn * latsgn
932 if (swapp .lt. 0)
then
933 call swap(salp1, salp2)
934 call swap(calp1, calp2)
935 if (scalp)
call swap(mm12, mm21)
938 salp1 = salp1 * swapp * lonsgn
939 calp1 = calp1 * swapp * latsgn
940 salp2 = salp2 * swapp * lonsgn
941 calp2 = calp2 * swapp * latsgn
943 azi1 = atn2dx(salp1, calp1)
944 azi2 = atn2dx(salp2, calp2)
971 subroutine area(a, f, lats, lons, n, AA, PP)
974 double precision a, f, lats(n), lons(n)
976 double precision AA, PP
978 integer i, omask, cross, trnsit
979 double precision s12, azi1, azi2, dummy, SS12, b, e2, c2, area0,
980 + atanhx, aacc(2), pacc(2)
982 double precision dblmin, dbleps, pi, degree, tiny,
983 + tol0, tol1, tol2, tolb, xthrsh
984 integer digits, maxit1, maxit2
986 common /geocom/ dblmin, dbleps, pi, degree, tiny,
987 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
994 call invers(a, f, lats(i+1), lons(i+1),
995 + lats(mod(i+1,n)+1), lons(mod(i+1,n)+1),
996 + s12, azi1, azi2, omask, dummy, dummy, dummy, dummy, ss12)
997 call accadd(pacc, s12)
998 call accadd(aacc, -ss12)
999 cross = cross + trnsit(lons(i+1), lons(mod(i+1,n)+1))
1006 else if (e2 .gt. 0)
then
1007 c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
1009 c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
1012 call accrem(aacc, area0)
1013 if (mod(abs(cross), 2) .eq. 1)
then
1014 if (aacc(1) .lt. 0)
then
1015 call accadd(aacc, +area0/2)
1017 call accadd(aacc, -area0/2)
1020 if (aacc(1) .gt. area0/2)
then
1021 call accadd(aacc, -area0)
1022 else if (aacc(1) .le. -area0/2)
then
1023 call accadd(aacc, +area0)
1038 subroutine geover(major, minor, patch)
1040 integer major, minor, patch
1052 double precision dblmin, dbleps, pi, degree, tiny,
1053 + tol0, tol1, tol2, tolb, xthrsh
1054 integer digits, maxit1, maxit2
1057 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1058 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1062 double precision dblmin, dbleps, pi, degree, tiny,
1063 + tol0, tol1, tol2, tolb, xthrsh
1064 integer digits, maxit1, maxit2
1066 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1067 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1070 dblmin = 0.5d0**1022
1071 dbleps = 0.5d0**(digits-1)
1073 pi = atan2(0d0, -1d0)
1079 tiny = 0.5d0**((1022+1)/3)
1088 xthrsh = 1000 * tol2
1090 maxit2 = maxit1 + digits + 10
1097 subroutine lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1098 + cbet1, cbet2, omask,
1099 + s12b, m12b, m0, MM12, MM21, ep2, Ca)
1101 double precision eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1105 double precision s12b, m12b, m0, MM12, MM21
1107 double precision Ca(*)
1109 integer ord, nC1, nC2
1110 parameter (ord = 6, nc1 = ord, nc2 = ord)
1112 double precision A1m1f, A2m1f, TrgSum
1113 double precision m0x, J12, A1, A2, B1, B2, csig12, t, Cb(nC2)
1114 logical distp, redlp, scalp
1120 distp = (mod(omask/16, 2) .eq. 1)
1121 redlp = (mod(omask/2, 2) .eq. 1)
1122 scalp = (mod(omask/4, 2) .eq. 1)
1129 if (distp .or. redlp .or. scalp)
then
1132 if (redlp .or. scalp)
then
1141 b1 = trgsum(.true., ssig2, csig2, ca, nc1) -
1142 + trgsum(.true., ssig1, csig1, ca, nc1)
1144 s12b = a1 * (sig12 + b1)
1145 if (redlp .or. scalp)
then
1146 b2 = trgsum(.true., ssig2, csig2, cb, nc2) -
1147 + trgsum(.true., ssig1, csig1, cb, nc2)
1148 j12 = m0x * sig12 + (a1 * b1 - a2 * b2)
1150 else if (redlp .or. scalp)
then
1153 cb(l) = a1 * ca(l) - a2 * cb(l)
1155 j12 = m0x * sig12 + (trgsum(.true., ssig2, csig2, cb, nc2) -
1156 + trgsum(.true., ssig1, csig1, cb, nc2))
1163 m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
1164 + csig1 * csig2 * j12
1167 csig12 = csig1 * csig2 + ssig1 * ssig2
1168 t = ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2)
1169 mm12 = csig12 + (t * ssig2 - csig2 * j12) * ssig1 / dn1
1170 mm21 = csig12 - (t * ssig1 - csig1 * j12) * ssig2 / dn2
1176 double precision function astrd(x, y)
1180 double precision x, y
1182 double precision cbrt
1183 double precision k, p, q, r, S, r2, r3, disc, u,
1184 + t3, t, ang, v, uv, w
1189 if ( .not. (q .eq. 0 .and. r .lt. 0) )
then
1198 disc = s * (s + 2 * r3)
1200 if (disc .ge. 0)
then
1216 if (t .ne. 0) u = u + t + r2 / t
1219 ang = atan2(sqrt(-disc), -(s + r3))
1222 u = u + 2 * r * cos(ang / 3)
1234 w = (uv - q) / (2 * v)
1238 k = uv / (sqrt(uv + w**2) + w)
1250 double precision function invsta(sbet1, cbet1, dn1,
1251 + sbet2, cbet2, dn2, lam12, slam12, clam12, f, A3x,
1252 + salp1, calp1, salp2, calp2, dnm,
1258 double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2,
1259 + lam12, slam12, clam12, f, A3x(*)
1261 double precision salp1, calp1, salp2, calp2, dnm
1263 double precision Ca(*)
1265 double precision hypotx, A3f, Astrd
1267 double precision f1, e2, ep2, n, etol2, k2, eps, sig12,
1268 + sbet12, cbet12, sbt12a, omg12, somg12, comg12, ssig12, csig12,
1269 + x, y, lamscl, betscl, cbt12a, bt12a, m12b, m0, dummy,
1270 + k, omg12a, sbetm2, lam12x
1272 double precision dblmin, dbleps, pi, degree, tiny,
1273 + tol0, tol1, tol2, tolb, xthrsh
1274 integer digits, maxit1, maxit2
1276 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1277 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1293 etol2 = 0.1d0 * tol2 /
1294 + sqrt( max(0.001d0, abs(f)) * min(1d0, 1 - f/2) / 2 )
1299 sbet12 = sbet2 * cbet1 - cbet2 * sbet1
1300 cbet12 = cbet2 * cbet1 + sbet2 * sbet1
1301 sbt12a = sbet2 * cbet1 + cbet2 * sbet1
1303 shortp = cbet12 .ge. 0 .and. sbet12 .lt. 0.5d0 .and.
1304 + cbet2 * lam12 .lt. 0.5d0
1307 sbetm2 = (sbet1 + sbet2)**2
1310 sbetm2 = sbetm2 / (sbetm2 + (cbet1 + cbet2)**2)
1311 dnm = sqrt(1 + ep2 * sbetm2)
1312 omg12 = lam12 / (f1 * dnm)
1320 salp1 = cbet2 * somg12
1321 if (comg12 .ge. 0)
then
1322 calp1 = sbet12 + cbet2 * sbet1 * somg12**2 / (1 + comg12)
1324 calp1 = sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12)
1327 ssig12 = hypotx(salp1, calp1)
1328 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12
1330 if (shortp .and. ssig12 .lt. etol2)
then
1332 salp2 = cbet1 * somg12
1333 if (comg12 .ge. 0)
then
1334 calp2 = somg12**2 / (1 + comg12)
1338 calp2 = sbet12 - cbet1 * sbet2 * calp2
1339 call norm2x(salp2, calp2)
1341 sig12 = atan2(ssig12, csig12)
1342 else if (abs(n) .gt. 0.1d0 .or. csig12 .ge. 0 .or.
1343 + ssig12 .ge. 6 * abs(n) * pi * cbet1**2)
then
1348 lam12x = atan2(-slam12, -clam12)
1354 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
1355 lamscl = f * cbet1 * a3f(eps, a3x) * pi
1356 betscl = lamscl * cbet1
1361 cbt12a = cbet2 * cbet1 - sbet2 * sbet1
1362 bt12a = atan2(sbt12a, cbt12a)
1365 call lengs(n, pi + bt12a,
1366 + sbet1, -cbet1, dn1, sbet2, cbet2, dn2, cbet1, cbet2, 2,
1367 + dummy, m12b, m0, dummy, dummy, ep2, ca)
1368 x = -1 + m12b / (cbet1 * cbet2 * m0 * pi)
1369 if (x .lt. -0.01d0)
then
1372 betscl = -f * cbet1**2 * pi
1374 lamscl = betscl / cbet1
1378 if (y .gt. -tol1 .and. x .gt. -1 - xthrsh)
then
1381 salp1 = min(1d0, -x)
1382 calp1 = - sqrt(1 - salp1**2)
1384 if (x .gt. -tol1)
then
1389 calp1 = max(calp1, x)
1390 salp1 = sqrt(1 - calp1**2)
1429 omg12a = -x * k/(1 + k)
1431 omg12a = -y * (1 + k)/k
1433 omg12a = lamscl * omg12a
1434 somg12 = sin(omg12a)
1435 comg12 = -cos(omg12a)
1437 salp1 = cbet2 * somg12
1438 calp1 = sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12)
1442 if (.not. (salp1 .le. 0))
then
1443 call norm2x(salp1, calp1)
1453 double precision function lam12f(sbet1, cbet1, dn1,
1454 + sbet2, cbet2, dn2, salp1, calp1, slm120, clm120, f, A3x, C3x,
1455 + salp2, calp2, sig12, ssig1, csig1, ssig2, csig2, eps,
1456 + domg12, diffp, dlam12, Ca)
1458 double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2,
1459 + salp1, calp1, slm120, clm120, f, A3x(*), C3x(*)
1462 double precision salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
1465 double precision dlam12
1467 double precision Ca(*)
1470 parameter(ord = 6, nc3 = ord)
1472 double precision hypotx, A3f, TrgSum
1474 double precision f1, e2, ep2, salp0, calp0,
1475 + somg1, comg1, somg2, comg2, somg12, comg12,
1476 + lam12, eta, b312, k2, dummy
1478 double precision dblmin, dbleps, pi, degree, tiny,
1479 + tol0, tol1, tol2, tolb, xthrsh
1480 integer digits, maxit1, maxit2
1482 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1483 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1490 if (sbet1 .eq. 0 .and. calp1 .eq. 0) calp1 = -tiny
1493 salp0 = salp1 * cbet1
1495 calp0 = hypotx(calp1, salp1 * sbet1)
1500 somg1 = salp0 * sbet1
1501 csig1 = calp1 * cbet1
1503 call norm2x(ssig1, csig1)
1510 if (cbet2 .ne. cbet1)
then
1511 salp2 = salp0 / cbet2
1519 if (cbet2 .ne. cbet1 .or. abs(sbet2) .ne. -sbet1)
then
1520 if (cbet1 .lt. -sbet1)
then
1521 calp2 = (cbet2 - cbet1) * (cbet1 + cbet2)
1523 calp2 = (sbet1 - sbet2) * (sbet1 + sbet2)
1525 calp2 = sqrt((calp1 * cbet1)**2 + calp2) / cbet2
1532 somg2 = salp0 * sbet2
1533 csig2 = calp2 * cbet2
1535 call norm2x(ssig2, csig2)
1539 sig12 = atan2(0d0 + max(0d0, csig1 * ssig2 - ssig1 * csig2),
1540 + csig1 * csig2 + ssig1 * ssig2)
1543 somg12 = 0d0 + max(0d0, comg1 * somg2 - somg1 * comg2)
1544 comg12 = comg1 * comg2 + somg1 * somg2
1546 eta = atan2(somg12 * clm120 - comg12 * slm120,
1547 + comg12 * clm120 + somg12 * slm120)
1549 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
1550 call c3f(eps, c3x, ca)
1551 b312 = (trgsum(.true., ssig2, csig2, ca, nc3-1) -
1552 + trgsum(.true., ssig1, csig1, ca, nc3-1))
1553 domg12 = -f * a3f(eps, a3x) * salp0 * (sig12 + b312)
1554 lam12 = eta + domg12
1557 if (calp2 .eq. 0)
then
1558 dlam12 = - 2 * f1 * dn1 / sbet1
1560 call lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1562 + dummy, dlam12, dummy, dummy, dummy, ep2, ca)
1563 dlam12 = dlam12 * f1 / (calp2 * cbet2)
1571 double precision function a3f(eps, A3x)
1573 integer ord, nA3, nA3x
1574 parameter(ord = 6, na3 = ord, na3x = na3)
1577 double precision eps
1579 double precision A3x(0: nA3x-1)
1581 double precision polval
1582 A3f = polval(na3 - 1, a3x, eps)
1587 subroutine c3f(eps, C3x, c)
1590 integer ord, nC3, nC3x
1591 parameter(ord = 6, nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2)
1594 double precision eps, C3x(0:nC3x-1)
1596 double precision c(nC3-1)
1599 double precision mult, polval
1603 do 10 l = 1, nc3 - 1
1606 c(l) = mult * polval(m, c3x(o), eps)
1613 subroutine c4f(eps, C4x, c)
1616 integer ord, nC4, nC4x
1617 parameter(ord = 6, nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
1620 double precision eps, C4x(0:nC4x-1)
1622 double precision c(0:nC4-1)
1625 double precision mult, polval
1629 do 10 l = 0, nc4 - 1
1631 c(l) = mult * polval(m, c4x(o), eps)
1639 double precision function a1m1f(eps)
1642 double precision eps
1645 integer ord, nA1, o, m
1646 parameter(ord = 6, na1 = ord)
1647 double precision polval, coeff(nA1/2 + 2)
1648 data coeff /1, 4, 64, 0, 256/
1652 t = polval(m, coeff(o), eps**2) / coeff(o + m + 1)
1653 a1m1f = (t + eps) / (1 - eps)
1658 subroutine c1f(eps, c)
1661 parameter(ord = 6, nc1 = ord)
1664 double precision eps
1666 double precision c(nC1)
1668 double precision eps2, d
1670 double precision polval, coeff((nC1**2 + 7*nC1 - 2*(nC1/2))/4)
1673 + -9, 64, -128, 2048,
1684 c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
1692 subroutine c1pf(eps, c)
1695 parameter(ord = 6, nc1p = ord)
1698 double precision eps
1700 double precision c(nC1p)
1702 double precision eps2, d
1704 double precision polval, coeff((nC1p**2 + 7*nC1p - 2*(nC1p/2))/4)
1706 + 205, -432, 768, 1536,
1707 + 4005, -4736, 3840, 12288,
1709 + -7173, 2695, 7680,
1718 c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
1727 double precision function a2m1f(eps)
1729 double precision eps
1732 integer ord, nA2, o, m
1733 parameter(ord = 6, na2 = ord)
1734 double precision polval, coeff(nA2/2 + 2)
1735 data coeff /-11, -28, -192, 0, 256/
1739 t = polval(m, coeff(o), eps**2) / coeff(o + m + 1)
1740 a2m1f = (t - eps) / (1 + eps)
1745 subroutine c2f(eps, c)
1748 parameter(ord = 6, nc2 = ord)
1751 double precision eps
1753 double precision c(nC2)
1755 double precision eps2, d
1757 double precision polval, coeff((nC2**2 + 7*nC2 - 2*(nC2/2))/4)
1760 + 35, 64, 384, 2048,
1771 c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
1779 subroutine a3cof(n, A3x)
1781 integer ord, nA3, nA3x
1782 parameter(ord = 6, na3 = ord, na3x = na3)
1787 double precision A3x(0:nA3x-1)
1790 double precision polval, coeff((nA3**2 + 7*nA3 - 2*(nA3/2))/4)
1801 do 10 j = na3 - 1, 0, -1
1802 m = min(na3 - j - 1, j)
1803 a3x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
1811 subroutine c3cof(n, C3x)
1813 integer ord, nC3, nC3x
1814 parameter(ord = 6, nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2)
1819 double precision C3x(0:nC3x-1)
1821 integer o, m, l, j, k
1822 double precision polval,
1823 + coeff(((nc3-1)*(nc3**2 + 7*nc3 - 2*(nc3/2)))/8)
1843 do 20 l = 1, nc3 - 1
1844 do 10 j = nc3 - 1, l, -1
1845 m = min(nc3 - j - 1, j)
1846 c3x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
1855 subroutine c4cof(n, C4x)
1857 integer ord, nC4, nC4x
1858 parameter(ord = 6, nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
1863 double precision C4x(0:nC4x-1)
1865 integer o, m, l, j, k
1866 double precision polval, coeff((nC4 * (nC4 + 1) * (nC4 + 5)) / 6)
1868 + 97, 15015, 1088, 156, 45045, -224, -4784, 1573, 45045,
1869 + -10656, 14144, -4576, -858, 45045,
1870 + 64, 624, -4576, 6864, -3003, 15015,
1871 + 100, 208, 572, 3432, -12012, 30030, 45045,
1872 + 1, 9009, -2944, 468, 135135, 5792, 1040, -1287, 135135,
1873 + 5952, -11648, 9152, -2574, 135135,
1874 + -64, -624, 4576, -6864, 3003, 135135,
1875 + 8, 10725, 1856, -936, 225225, -8448, 4992, -1144, 225225,
1876 + -1440, 4160, -4576, 1716, 225225,
1877 + -136, 63063, 1024, -208, 105105,
1878 + 3584, -3328, 1144, 315315,
1879 + -128, 135135, -2560, 832, 405405, 128, 99099/
1883 do 20 l = 0, nc4 - 1
1884 do 10 j = nc4 - 1, l, -1
1886 c4x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
1895 double precision function sumx(u, v, t)
1897 double precision u, v
1901 double precision up, vpp
1912 double precision function remx(x, y)
1916 double precision x, y
1919 if (remx .lt. -y/2)
then
1921 else if (remx .gt. +y/2)
then
1928 double precision function angnm(x)
1932 double precision remx
1933 angnm = remx(x, 360d0)
1934 if (angnm .eq. -180)
then
1941 double precision function latfix(x)
1946 if (.not. (abs(x) .gt. 90))
return
1948 latfix = sqrt(90 - abs(x))
1953 double precision function angdif(x, y, e)
1960 double precision x, y
1964 double precision d, t, sumx, AngNm
1965 d = angnm(sumx(angnm(-x), angnm(y), t))
1966 if (d .eq. 180 .and. t .gt. 0)
then
1969 angdif = sumx(d, t, e)
1974 double precision function angrnd(x)
1984 double precision y, z
1988 if (y .lt. z) y = z - (z - y)
1990 if (x .eq. 0) angrnd = 0
1995 subroutine swap(x, y)
1997 double precision x, y
2007 double precision function hypotx(x, y)
2009 double precision x, y
2012 hypotx = sqrt(x**2 + y**2)
2017 subroutine norm2x(x, y)
2019 double precision x, y
2021 double precision hypotx, r
2029 double precision function log1px(x)
2033 double precision y, z
2039 log1px = x * log(y) / z
2045 double precision function atanhx(x)
2050 double precision log1px, y
2052 y = log1px(2 * y/(1 - y))/2
2058 double precision function cbrt(x)
2062 cbrt = sign(abs(x)**(1/3d0), x)
2067 double precision function trgsum(sinp, sinx, cosx, c, n)
2076 double precision sinx, cosx, c(n)
2078 double precision ar, y0, y1
2082 ar = 2 * (cosx - sinx) * (cosx + sinx)
2084 if (mod(n, 2) .eq. 1)
then
2095 y1 = ar * y0 - y1 + c(k)
2096 y0 = ar * y1 - y0 + c(k-1)
2100 trgsum = 2 * sinx * cosx * y0
2103 trgsum = cosx * (y0 - y1)
2109 integer function trnsit(lon1, lon2)
2111 double precision lon1, lon2
2113 double precision lon1x, lon2x, lon12, AngNm, AngDif, e
2116 lon12 = angdif(lon1x, lon2x, e)
2118 if (lon1x .le. 0 .and. lon2x .gt. 0 .and. lon12 .gt. 0)
then
2120 else if (lon2x .le. 0 .and. lon1x .gt. 0 .and. lon12 .lt. 0)
then
2127 subroutine accini(s)
2130 double precision s(2)
2138 subroutine accadd(s, y)
2143 double precision s(2)
2145 double precision z, u, sumx
2146 z = sumx(y, s(2), u)
2147 s(1) = sumx(z, s(1), s(2))
2148 if (s(1) .eq. 0)
then
2157 subroutine accrem(s, y)
2162 double precision s(2)
2164 double precision remx
2165 s(1) = remx(s(1), y)
2171 subroutine sncsdx(x, sinx, cosx)
2176 double precision sinx, cosx
2178 double precision dblmin, dbleps, pi, degree, tiny,
2179 + tol0, tol1, tol2, tolb, xthrsh
2180 integer digits, maxit1, maxit2
2182 common /geocom/ dblmin, dbleps, pi, degree, tiny,
2183 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
2185 double precision r, s, c
2189 r = (r - 90 * q) * degree
2198 else if (q .eq. 1)
then
2201 else if (q .eq. 2)
then
2218 double precision function atn2dx(y, x)
2220 double precision x, y
2222 double precision dblmin, dbleps, pi, degree, tiny,
2223 + tol0, tol1, tol2, tolb, xthrsh
2224 integer digits, maxit1, maxit2
2226 common /geocom/ dblmin, dbleps, pi, degree, tiny,
2227 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
2229 double precision xx, yy
2231 if (abs(y) .gt. abs(x)) then
2244 atn2dx = atan2(yy, xx) / degree
2247 atn2dx = 180 - atn2dx
2249 atn2dx = -180 - atn2dx
2251 else if (q .eq. 2)
then
2252 atn2dx = 90 - atn2dx
2253 else if (q .eq. 3)
then
2254 atn2dx = -90 + atn2dx
2260 double precision function polval(N, p, x)
2263 double precision p(0:N), x
2272 polval = polval * x + p(i)