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)
969 subroutine area(a, f, lats, lons, n, AA, PP)
972 double precision a, f, lats(n), lons(n)
974 double precision AA, PP
976 integer i, omask, cross, trnsit
977 double precision s12, azi1, azi2, dummy, SS12, b, e2, c2, area0,
978 + atanhx, aacc(2), pacc(2)
980 double precision dblmin, dbleps, pi, degree, tiny,
981 + tol0, tol1, tol2, tolb, xthrsh
982 integer digits, maxit1, maxit2
984 common /geocom/ dblmin, dbleps, pi, degree, tiny,
985 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
992 call invers(a, f, lats(i+1), lons(i+1),
993 + lats(mod(i+1,n)+1), lons(mod(i+1,n)+1),
994 + s12, azi1, azi2, omask, dummy, dummy, dummy, dummy, ss12)
995 call accadd(pacc, s12)
996 call accadd(aacc, -ss12)
997 cross = cross + trnsit(lons(i+1), lons(mod(i+1,n)+1))
1004 else if (e2 .gt. 0)
then 1005 c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
1007 c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
1010 if (mod(abs(cross), 2) .eq. 1)
then 1011 if (aacc(1) .lt. 0)
then 1012 call accadd(aacc, +area0/2)
1014 call accadd(aacc, -area0/2)
1017 if (aacc(1) .gt. area0/2)
then 1018 call accadd(aacc, -area0)
1019 else if (aacc(1) .le. -area0/2)
then 1020 call accadd(aacc, +area0)
1035 subroutine geover(major, minor, patch)
1037 integer major, minor, patch
1049 double precision dblmin, dbleps, pi, degree, tiny,
1050 + tol0, tol1, tol2, tolb, xthrsh
1051 integer digits, maxit1, maxit2
1054 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1055 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1059 double precision dblmin, dbleps, pi, degree, tiny,
1060 + tol0, tol1, tol2, tolb, xthrsh
1061 integer digits, maxit1, maxit2
1063 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1064 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1067 dblmin = 0.5d0**1022
1068 dbleps = 0.5d0**(digits-1)
1070 pi = atan2(0d0, -1d0)
1076 tiny = 0.5d0**((1022+1)/3)
1085 xthrsh = 1000 * tol2
1087 maxit2 = maxit1 + digits + 10
1094 subroutine lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1095 + cbet1, cbet2, omask,
1096 + s12b, m12b, m0, MM12, MM21, ep2, Ca)
1098 double precision eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1102 double precision s12b, m12b, m0, MM12, MM21
1104 double precision Ca(*)
1106 integer ord, nC1, nC2
1107 parameter (ord = 6, nc1 = ord, nc2 = ord)
1109 double precision A1m1f, A2m1f, TrgSum
1110 double precision m0x, J12, A1, A2, B1, B2, csig12, t, Cb(nC2)
1111 logical distp, redlp, scalp
1117 distp = (mod(omask/16, 2) .eq. 1)
1118 redlp = (mod(omask/2, 2) .eq. 1)
1119 scalp = (mod(omask/4, 2) .eq. 1)
1126 if (distp .or. redlp .or. scalp)
then 1129 if (redlp .or. scalp)
then 1138 b1 = trgsum(.true., ssig2, csig2, ca, nc1) -
1139 + trgsum(.true., ssig1, csig1, ca, nc1)
1141 s12b = a1 * (sig12 + b1)
1142 if (redlp .or. scalp)
then 1143 b2 = trgsum(.true., ssig2, csig2, cb, nc2) -
1144 + trgsum(.true., ssig1, csig1, cb, nc2)
1145 j12 = m0x * sig12 + (a1 * b1 - a2 * b2)
1147 else if (redlp .or. scalp)
then 1150 cb(l) = a1 * ca(l) - a2 * cb(l)
1152 j12 = m0x * sig12 + (trgsum(.true., ssig2, csig2, cb, nc2) -
1153 + trgsum(.true., ssig1, csig1, cb, nc2))
1160 m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
1161 + csig1 * csig2 * j12
1164 csig12 = csig1 * csig2 + ssig1 * ssig2
1165 t = ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2)
1166 mm12 = csig12 + (t * ssig2 - csig2 * j12) * ssig1 / dn1
1167 mm21 = csig12 - (t * ssig1 - csig1 * j12) * ssig2 / dn2
1173 double precision function astrd(x, y)
1177 double precision x, y
1179 double precision cbrt
1180 double precision k, p, q, r, S, r2, r3, disc, u,
1181 + t3, t, ang, v, uv, w
1186 if ( .not. (q .eq. 0 .and. r .lt. 0) )
then 1195 disc = s * (s + 2 * r3)
1197 if (disc .ge. 0)
then 1213 if (t .ne. 0) u = u + t + r2 / t
1216 ang = atan2(sqrt(-disc), -(s + r3))
1219 u = u + 2 * r * cos(ang / 3)
1231 w = (uv - q) / (2 * v)
1235 k = uv / (sqrt(uv + w**2) + w)
1247 double precision function invsta(sbet1, cbet1, dn1,
1248 + sbet2, cbet2, dn2, lam12, slam12, clam12, f, A3x,
1249 + salp1, calp1, salp2, calp2, dnm,
1255 double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2,
1256 + lam12, slam12, clam12, f, A3x(*)
1258 double precision salp1, calp1, salp2, calp2, dnm
1260 double precision Ca(*)
1262 double precision hypotx, A3f, Astrd
1264 double precision f1, e2, ep2, n, etol2, k2, eps, sig12,
1265 + sbet12, cbet12, sbt12a, omg12, somg12, comg12, ssig12, csig12,
1266 + x, y, lamscl, betscl, cbt12a, bt12a, m12b, m0, dummy,
1267 + k, omg12a, sbetm2, lam12x
1269 double precision dblmin, dbleps, pi, degree, tiny,
1270 + tol0, tol1, tol2, tolb, xthrsh
1271 integer digits, maxit1, maxit2
1273 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1274 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1290 etol2 = 0.1d0 * tol2 /
1291 + sqrt( max(0.001d0, abs(f)) * min(1d0, 1 - f/2) / 2 )
1296 sbet12 = sbet2 * cbet1 - cbet2 * sbet1
1297 cbet12 = cbet2 * cbet1 + sbet2 * sbet1
1298 sbt12a = sbet2 * cbet1 + cbet2 * sbet1
1300 shortp = cbet12 .ge. 0 .and. sbet12 .lt. 0.5d0 .and.
1301 + cbet2 * lam12 .lt. 0.5d0
1304 sbetm2 = (sbet1 + sbet2)**2
1307 sbetm2 = sbetm2 / (sbetm2 + (cbet1 + cbet2)**2)
1308 dnm = sqrt(1 + ep2 * sbetm2)
1309 omg12 = lam12 / (f1 * dnm)
1317 salp1 = cbet2 * somg12
1318 if (comg12 .ge. 0)
then 1319 calp1 = sbet12 + cbet2 * sbet1 * somg12**2 / (1 + comg12)
1321 calp1 = sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12)
1324 ssig12 = hypotx(salp1, calp1)
1325 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12
1327 if (shortp .and. ssig12 .lt. etol2)
then 1329 salp2 = cbet1 * somg12
1330 if (comg12 .ge. 0)
then 1331 calp2 = somg12**2 / (1 + comg12)
1335 calp2 = sbet12 - cbet1 * sbet2 * calp2
1336 call norm2x(salp2, calp2)
1338 sig12 = atan2(ssig12, csig12)
1339 else if (abs(n) .gt. 0.1d0 .or. csig12 .ge. 0 .or.
1340 + ssig12 .ge. 6 * abs(n) * pi * cbet1**2)
then 1345 lam12x = atan2(-slam12, -clam12)
1351 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
1352 lamscl = f * cbet1 * a3f(eps, a3x) * pi
1353 betscl = lamscl * cbet1
1358 cbt12a = cbet2 * cbet1 - sbet2 * sbet1
1359 bt12a = atan2(sbt12a, cbt12a)
1362 call lengs(n, pi + bt12a,
1363 + sbet1, -cbet1, dn1, sbet2, cbet2, dn2, cbet1, cbet2, 2,
1364 + dummy, m12b, m0, dummy, dummy, ep2, ca)
1365 x = -1 + m12b / (cbet1 * cbet2 * m0 * pi)
1366 if (x .lt. -0.01d0)
then 1369 betscl = -f * cbet1**2 * pi
1371 lamscl = betscl / cbet1
1375 if (y .gt. -tol1 .and. x .gt. -1 - xthrsh)
then 1378 salp1 = min(1d0, -x)
1379 calp1 = - sqrt(1 - salp1**2)
1381 if (x .gt. -tol1)
then 1386 calp1 = max(calp1, x)
1387 salp1 = sqrt(1 - calp1**2)
1426 omg12a = -x * k/(1 + k)
1428 omg12a = -y * (1 + k)/k
1430 omg12a = lamscl * omg12a
1431 somg12 = sin(omg12a)
1432 comg12 = -cos(omg12a)
1434 salp1 = cbet2 * somg12
1435 calp1 = sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12)
1439 if (.not. (salp1 .le. 0))
then 1440 call norm2x(salp1, calp1)
1450 double precision function lam12f(sbet1, cbet1, dn1,
1451 + sbet2, cbet2, dn2, salp1, calp1, slm120, clm120, f, A3x, C3x,
1452 + salp2, calp2, sig12, ssig1, csig1, ssig2, csig2, eps,
1453 + domg12, diffp, dlam12, Ca)
1455 double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2,
1456 + salp1, calp1, slm120, clm120, f, A3x(*), C3x(*)
1459 double precision salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
1462 double precision dlam12
1464 double precision Ca(*)
1467 parameter(ord = 6, nc3 = ord)
1469 double precision hypotx, A3f, TrgSum
1471 double precision f1, e2, ep2, salp0, calp0,
1472 + somg1, comg1, somg2, comg2, somg12, comg12,
1473 + lam12, eta, b312, k2, dummy
1475 double precision dblmin, dbleps, pi, degree, tiny,
1476 + tol0, tol1, tol2, tolb, xthrsh
1477 integer digits, maxit1, maxit2
1479 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1480 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1487 if (sbet1 .eq. 0 .and. calp1 .eq. 0) calp1 = -tiny
1490 salp0 = salp1 * cbet1
1492 calp0 = hypotx(calp1, salp1 * sbet1)
1497 somg1 = salp0 * sbet1
1498 csig1 = calp1 * cbet1
1500 call norm2x(ssig1, csig1)
1507 if (cbet2 .ne. cbet1)
then 1508 salp2 = salp0 / cbet2
1516 if (cbet2 .ne. cbet1 .or. abs(sbet2) .ne. -sbet1)
then 1517 if (cbet1 .lt. -sbet1)
then 1518 calp2 = (cbet2 - cbet1) * (cbet1 + cbet2)
1520 calp2 = (sbet1 - sbet2) * (sbet1 + sbet2)
1522 calp2 = sqrt((calp1 * cbet1)**2 + calp2) / cbet2
1529 somg2 = salp0 * sbet2
1530 csig2 = calp2 * cbet2
1532 call norm2x(ssig2, csig2)
1536 sig12 = atan2(0d0 + max(0d0, csig1 * ssig2 - ssig1 * csig2),
1537 + csig1 * csig2 + ssig1 * ssig2)
1540 somg12 = 0d0 + max(0d0, comg1 * somg2 - somg1 * comg2)
1541 comg12 = comg1 * comg2 + somg1 * somg2
1543 eta = atan2(somg12 * clm120 - comg12 * slm120,
1544 + comg12 * clm120 + somg12 * slm120);
1546 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
1547 call c3f(eps, c3x, ca)
1548 b312 = (trgsum(.true., ssig2, csig2, ca, nc3-1) -
1549 + trgsum(.true., ssig1, csig1, ca, nc3-1))
1550 domg12 = -f * a3f(eps, a3x) * salp0 * (sig12 + b312)
1551 lam12 = eta + domg12
1554 if (calp2 .eq. 0)
then 1555 dlam12 = - 2 * f1 * dn1 / sbet1
1557 call lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1559 + dummy, dlam12, dummy, dummy, dummy, ep2, ca)
1560 dlam12 = dlam12 * f1 / (calp2 * cbet2)
1568 double precision function a3f(eps, A3x)
1570 integer ord, nA3, nA3x
1571 parameter(ord = 6, na3 = ord, na3x = na3)
1574 double precision eps
1576 double precision A3x(0: nA3x-1)
1578 double precision polval
1579 A3f = polval(na3 - 1, a3x, eps)
1584 subroutine c3f(eps, C3x, c)
1587 integer ord, nC3, nC3x
1588 parameter(ord = 6, nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2)
1591 double precision eps, C3x(0:nC3x-1)
1593 double precision c(nC3-1)
1596 double precision mult, polval
1600 do 10 l = 1, nc3 - 1
1603 c(l) = mult * polval(m, c3x(o), eps)
1610 subroutine c4f(eps, C4x, c)
1613 integer ord, nC4, nC4x
1614 parameter(ord = 6, nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
1617 double precision eps, C4x(0:nC4x-1)
1619 double precision c(0:nC4-1)
1622 double precision mult, polval
1626 do 10 l = 0, nc4 - 1
1628 c(l) = mult * polval(m, c4x(o), eps)
1636 double precision function a1m1f(eps)
1639 double precision eps
1642 integer ord, nA1, o, m
1643 parameter(ord = 6, na1 = ord)
1644 double precision polval, coeff(nA1/2 + 2)
1645 data coeff /1, 4, 64, 0, 256/
1649 t = polval(m, coeff(o), eps**2) / coeff(o + m + 1)
1650 a1m1f = (t + eps) / (1 - eps)
1655 subroutine c1f(eps, c)
1658 parameter(ord = 6, nc1 = ord)
1661 double precision eps
1663 double precision c(nC1)
1665 double precision eps2, d
1667 double precision polval, coeff((nC1**2 + 7*nC1 - 2*(nC1/2))/4)
1670 + -9, 64, -128, 2048,
1681 c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
1689 subroutine c1pf(eps, c)
1692 parameter(ord = 6, nc1p = ord)
1695 double precision eps
1697 double precision c(nC1p)
1699 double precision eps2, d
1701 double precision polval, coeff((nC1p**2 + 7*nC1p - 2*(nC1p/2))/4)
1703 + 205, -432, 768, 1536,
1704 + 4005, -4736, 3840, 12288,
1706 + -7173, 2695, 7680,
1715 c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
1724 double precision function a2m1f(eps)
1726 double precision eps
1729 integer ord, nA2, o, m
1730 parameter(ord = 6, na2 = ord)
1731 double precision polval, coeff(nA2/2 + 2)
1732 data coeff /-11, -28, -192, 0, 256/
1736 t = polval(m, coeff(o), eps**2) / coeff(o + m + 1)
1737 a2m1f = (t - eps) / (1 + eps)
1742 subroutine c2f(eps, c)
1745 parameter(ord = 6, nc2 = ord)
1748 double precision eps
1750 double precision c(nC2)
1752 double precision eps2, d
1754 double precision polval, coeff((nC2**2 + 7*nC2 - 2*(nC2/2))/4)
1757 + 35, 64, 384, 2048,
1768 c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
1776 subroutine a3cof(n, A3x)
1778 integer ord, nA3, nA3x
1779 parameter(ord = 6, na3 = ord, na3x = na3)
1784 double precision A3x(0:nA3x-1)
1787 double precision polval, coeff((nA3**2 + 7*nA3 - 2*(nA3/2))/4)
1798 do 10 j = na3 - 1, 0, -1
1799 m = min(na3 - j - 1, j)
1800 a3x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
1808 subroutine c3cof(n, C3x)
1810 integer ord, nC3, nC3x
1811 parameter(ord = 6, nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2)
1816 double precision C3x(0:nC3x-1)
1818 integer o, m, l, j, k
1819 double precision polval,
1820 + coeff(((nc3-1)*(nc3**2 + 7*nc3 - 2*(nc3/2)))/8)
1840 do 20 l = 1, nc3 - 1
1841 do 10 j = nc3 - 1, l, -1
1842 m = min(nc3 - j - 1, j)
1843 c3x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
1852 subroutine c4cof(n, C4x)
1854 integer ord, nC4, nC4x
1855 parameter(ord = 6, nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
1860 double precision C4x(0:nC4x-1)
1862 integer o, m, l, j, k
1863 double precision polval, coeff((nC4 * (nC4 + 1) * (nC4 + 5)) / 6)
1865 + 97, 15015, 1088, 156, 45045, -224, -4784, 1573, 45045,
1866 + -10656, 14144, -4576, -858, 45045,
1867 + 64, 624, -4576, 6864, -3003, 15015,
1868 + 100, 208, 572, 3432, -12012, 30030, 45045,
1869 + 1, 9009, -2944, 468, 135135, 5792, 1040, -1287, 135135,
1870 + 5952, -11648, 9152, -2574, 135135,
1871 + -64, -624, 4576, -6864, 3003, 135135,
1872 + 8, 10725, 1856, -936, 225225, -8448, 4992, -1144, 225225,
1873 + -1440, 4160, -4576, 1716, 225225,
1874 + -136, 63063, 1024, -208, 105105,
1875 + 3584, -3328, 1144, 315315,
1876 + -128, 135135, -2560, 832, 405405, 128, 99099/
1880 do 20 l = 0, nc4 - 1
1881 do 10 j = nc4 - 1, l, -1
1883 c4x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
1892 double precision function sumx(u, v, t)
1894 double precision u, v
1898 double precision up, vpp
1909 double precision function angnm(x)
1913 angnm = mod(x, 360d0)
1914 if (angnm .le. -180)
then 1916 else if (angnm .gt. 180)
then 1923 double precision function latfix(x)
1928 if (.not. (abs(x) .gt. 90))
return 1930 latfix = sqrt(90 - abs(x))
1935 double precision function angdif(x, y, e)
1942 double precision x, y
1946 double precision d, t, sumx, AngNm
1947 d = angnm(sumx(angnm(-x), angnm(y), t))
1948 if (d .eq. 180 .and. t .gt. 0)
then 1951 angdif = sumx(d, t, e)
1956 double precision function angrnd(x)
1966 double precision y, z
1970 if (y .lt. z) y = z - (z - y)
1972 if (x .eq. 0) angrnd = 0
1977 subroutine swap(x, y)
1979 double precision x, y
1989 double precision function hypotx(x, y)
1991 double precision x, y
1994 hypotx = sqrt(x**2 + y**2)
1999 subroutine norm2x(x, y)
2001 double precision x, y
2003 double precision hypotx, r
2011 double precision function log1px(x)
2015 double precision y, z
2021 log1px = x * log(y) / z
2027 double precision function atanhx(x)
2032 double precision log1px, y
2034 y = log1px(2 * y/(1 - y))/2
2040 double precision function cbrt(x)
2044 cbrt = sign(abs(x)**(1/3d0), x)
2049 double precision function trgsum(sinp, sinx, cosx, c, n)
2058 double precision sinx, cosx, c(n)
2060 double precision ar, y0, y1
2064 ar = 2 * (cosx - sinx) * (cosx + sinx)
2066 if (mod(n, 2) .eq. 1)
then 2077 y1 = ar * y0 - y1 + c(k)
2078 y0 = ar * y1 - y0 + c(k-1)
2082 trgsum = 2 * sinx * cosx * y0
2085 trgsum = cosx * (y0 - y1)
2091 integer function trnsit(lon1, lon2)
2093 double precision lon1, lon2
2095 double precision lon1x, lon2x, lon12, AngNm, AngDif, e
2098 lon12 = angdif(lon1x, lon2x, e)
2100 if (lon1x .le. 0 .and. lon2x .gt. 0 .and. lon12 .gt. 0)
then 2102 else if (lon2x .le. 0 .and. lon1x .gt. 0 .and. lon12 .lt. 0)
then 2109 subroutine accini(s)
2112 double precision s(2)
2120 subroutine accadd(s, y)
2125 double precision s(2)
2127 double precision z, u, sumx
2128 z = sumx(y, s(2), u)
2129 s(1) = sumx(z, s(1), s(2))
2130 if (s(1) .eq. 0)
then 2139 subroutine sncsdx(x, sinx, cosx)
2144 double precision sinx, cosx
2146 double precision dblmin, dbleps, pi, degree, tiny,
2147 + tol0, tol1, tol2, tolb, xthrsh
2148 integer digits, maxit1, maxit2
2150 common /geocom/ dblmin, dbleps, pi, degree, tiny,
2151 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
2153 double precision r, s, c
2157 r = (r - 90 * q) * degree
2166 else if (q .eq. 1)
then 2169 else if (q .eq. 2)
then 2186 double precision function atn2dx(y, x)
2188 double precision x, y
2190 double precision dblmin, dbleps, pi, degree, tiny,
2191 + tol0, tol1, tol2, tolb, xthrsh
2192 integer digits, maxit1, maxit2
2194 common /geocom/ dblmin, dbleps, pi, degree, tiny,
2195 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
2197 double precision xx, yy
2199 if (abs(y) .gt. abs(x)) then
2212 atn2dx = atan2(yy, xx) / degree
2215 atn2dx = 180 - atn2dx
2217 atn2dx = -180 - atn2dx
2219 else if (q .eq. 2)
then 2220 atn2dx = 90 - atn2dx
2221 else if (q .eq. 3)
then 2222 atn2dx = -90 + atn2dx
2228 double precision function polval(N, p, x)
2231 double precision p(0:N), x
2240 polval = polval * x + p(i)