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, angnm2, angrnd, trgsum, a1m1f, a2m1f, a3f
203 logical arcmod, nowrap, arcp, redlp, scalp, areap
204 double precision e2, f1, ep2, n, b, c2,
205 + lon1x, azi1x, phi, alp1, 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, s12a, t, s, c, serr,
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 nowrap = 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 azi1x = angrnd(angnm(azi1))
256 alp1 = azi1x * degree
259 salp1 = csmgt(0d0, sin(alp1), azi1x .eq. -180)
260 calp1 = csmgt(0d0, cos(alp1), abs(azi1x) .eq. 90)
264 sbet1 = f1 * sin(phi)
265 cbet1 = csmgt(tiny, cos(phi), abs(lat1) .eq. 90)
266 call norm(sbet1, cbet1)
267 dn1 = sqrt(1 + ep2 * sbet1**2)
271 salp0 = salp1 * cbet1
274 calp0 = hypotx(calp1, salp1 * sbet1)
285 somg1 = salp0 * sbet1
286 csig1 = csmgt(cbet1 * calp1, 1d0, sbet1 .ne. 0 .or. calp1 .ne. 0)
289 call norm(ssig1, csig1)
293 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
297 b11 = trgsum(.true., ssig1, csig1, c1a, nc1)
301 stau1 = ssig1 * c + csig1 * s
302 ctau1 = csig1 * c - ssig1 * s
306 if (.not. arcmod)
call c1pf(eps, c1pa)
308 if (redlp .or. scalp)
then
311 b21 = trgsum(.true., ssig1, csig1, c2a, nc2)
318 call c3f(eps, c3x, c3a)
319 a3c = -f * salp0 * a3f(eps, a3x)
320 b31 = trgsum(.true., ssig1, csig1, c3a, nc3-1)
323 call c4f(eps, c4x, c4a)
325 a4 = a**2 * calp0 * salp0 * e2
326 b41 = trgsum(.false., ssig1, csig1, c4a, nc4)
335 sig12 = s12a12 * degree
337 s12a = s12a - 180 * aint(s12a / 180)
338 ssig12 = csmgt(0d0, sin(sig12), s12a .eq. 0)
339 csig12 = csmgt(0d0, cos(sig12), s12a .eq. 90)
344 tau12 = s12a12 / (b * (1 + a1m1))
348 b12 = - trgsum(.true.,
349 + stau1 * c + ctau1 * s, ctau1 * c - stau1 * s, c1pa, nc1p)
350 sig12 = tau12 - (b12 - b11)
353 if (abs(f) .gt. 0.01d0)
then
375 ssig2 = ssig1 * csig12 + csig1 * ssig12
376 csig2 = csig1 * csig12 - ssig1 * ssig12
377 b12 = trgsum(.true., ssig2, csig2, c1a, nc1)
378 serr = (1 + a1m1) * (sig12 + (b12 - b11)) - s12a12 / b
379 sig12 = sig12 - serr / sqrt(1 + k2 * ssig2**2)
387 ssig2 = ssig1 * csig12 + csig1 * ssig12
388 csig2 = csig1 * csig12 - ssig1 * ssig12
389 dn2 = sqrt(1 + k2 * ssig2**2)
390 if (arcmod .or. abs(f) .gt. 0.01d0)
391 + b12 = trgsum(.true., ssig2, csig2, c1a, nc1)
392 ab1 = (1 + a1m1) * (b12 - b11)
395 sbet2 = calp0 * ssig2
397 cbet2 = hypotx(salp0, calp0 * csig2)
398 if (cbet2 .eq. 0)
then
405 somg2 = salp0 * ssig2
410 calp2 = calp0 * csig2
413 + - (atan2(ssig2, csig2) - atan2(ssig1, csig1))
414 + + (atan2(somg2, comg2) - atan2(somg1, comg1)),
415 + atan2(somg2 * comg1 - comg2 * somg1,
416 + comg2 * comg1 + somg2 * somg1),
419 lam12 = omg12 + a3c *
420 + ( sig12 + (trgsum(.true., ssig2, csig2, c3a, nc3-1)
422 lon12 = lam12 / degree
425 lon2 = csmgt(lon1 + lon12, angnm(lon1x + angnm2(lon12)), nowrap)
426 lat2 = atan2(sbet2, f1 * cbet2) / degree
428 azi2 = 0 - atan2(-salp2, calp2) / degree
430 if (redlp .or. scalp)
then
431 b22 = trgsum(.true., ssig2, csig2, c2a, nc2)
432 ab2 = (1 + a2m1) * (b22 - b21)
433 j12 = (a1m1 - a2m1) * sig12 + (ab1 - ab2)
437 if (redlp) m12 = b * ((dn2 * (csig1 * ssig2) -
438 + dn1 * (ssig1 * csig2)) - csig1 * csig2 * j12)
440 t = k2 * (ssig2 - ssig1) * (ssig2 + ssig1) / (dn1 + dn2)
441 mm12 = csig12 + (t * ssig2 - csig2 * j12) * ssig1 / dn1
442 mm21 = csig12 - (t * ssig1 - csig1 * j12) * ssig2 / dn2
446 b42 = trgsum(.false., ssig2, csig2, c4a, nc4)
447 if (calp0 .eq. 0 .or. salp0 .eq. 0)
then
449 salp12 = salp2 * calp1 - calp2 * salp1
450 calp12 = calp2 * calp1 + salp2 * salp1
454 if (salp12 .eq. 0 .and. calp12 .lt. 0)
then
455 salp12 = tiny * calp1
467 salp12 = calp0 * salp0 *
468 + csmgt(csig1 * (1 - csig12) + ssig12 * ssig1,
469 + ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1),
471 calp12 = salp0**2 + calp0**2 * csig1 * csig2
473 ss12 = c2 * atan2(salp12, calp12) + a4 * (b42 - b41)
476 if (arcp) a12s12 = csmgt(b * ((1 + a1m1) * sig12 + ab1),
477 + sig12 / degree, arcmod)
525 subroutine invers(a, f, lat1, lon1, lat2, lon2,
526 + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
528 double precision a, f, lat1, lon1, lat2, lon2
531 double precision s12, azi1, azi2
533 double precision a12, m12, MM12, MM21, SS12
535 integer ord, nC1, nC2, nA3, nA3x, nC3, nC3x, nC4, nC4x
536 parameter(ord = 6, nc1 = ord, nc2 = ord, na3 = ord, na3x = na3,
537 + nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2,
538 + nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
539 double precision A3x(0:na3x-1), C3x(0:nc3x-1), C4x(0:nc4x-1),
540 + c1a(nc1), c2a(nc2), c3a(nc3-1), c4a(0:nc4-1)
542 double precision csmgt, atanhx, hypotx,
543 + angnm, angdif, angrnd, trgsum, lam12f, invsta
544 integer latsgn, lonsgn, swapp, numit
545 logical arcp, redlp, scalp, areap, merid, tripn, tripb
547 double precision e2, f1, ep2, n, b, c2,
548 + lat1x, lat2x, phi, salp0, calp0, k2, eps,
549 + salp1, calp1, ssig1, csig1, cbet1, sbet1, dbet1, dn1,
550 + salp2, calp2, ssig2, csig2, sbet2, cbet2, dbet2, dn2,
551 + slam12, clam12, salp12, calp12, omg12, lam12, lon12,
552 + salp1a, calp1a, salp1b, calp1b,
553 + dalp1, sdalp1, cdalp1, nsalp1, alp12, somg12, domg12,
554 + sig12, v, dv, dnm, dummy,
555 + a4, b41, b42, s12x, m12x, a12x
557 double precision dblmin, dbleps, pi, degree, tiny,
558 + tol0, tol1, tol2, tolb, xthrsh
559 integer digits, maxit1, maxit2
561 common /geocom/ dblmin, dbleps, pi, degree, tiny,
562 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
564 if (.not.init)
call geoini
573 arcp = mod(omask/1, 2) == 1
574 redlp = mod(omask/2, 2) == 1
575 scalp = mod(omask/4, 2) == 1
576 areap = mod(omask/8, 2) == 1
581 else if (e2 .gt. 0)
then
582 c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
584 c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
590 if (areap)
call c4cof(n, c4x)
595 lon12 = angdif(angnm(lon1), angnm(lon2))
597 lon12 = angrnd(lon12)
599 if (lon12 .ge. 0)
then
604 lon12 = lon12 * lonsgn
609 if (abs(lat1x) .ge. abs(lat2x))
then
614 if (swapp .lt. 0)
then
616 call swap(lat1x, lat2x)
619 if (lat1x .lt. 0)
then
624 lat1x = lat1x * latsgn
625 lat2x = lat2x * latsgn
640 sbet1 = f1 * sin(phi)
641 cbet1 = csmgt(tiny, cos(phi), lat1x .eq. -90)
642 call norm(sbet1, cbet1)
646 sbet2 = f1 * sin(phi)
647 cbet2 = csmgt(tiny, cos(phi), abs(lat2x) .eq. 90)
648 call norm(sbet2, cbet2)
659 if (cbet1 .lt. -sbet1)
then
660 if (cbet2 .eq. cbet1) sbet2 = sign(sbet1, sbet2)
662 if (abs(sbet2) .eq. -sbet1) cbet2 = cbet1
665 dn1 = sqrt(1 + ep2 * sbet1**2)
666 dn2 = sqrt(1 + ep2 * sbet2**2)
668 lam12 = lon12 * degree
670 if (lon12 .eq. 180) slam12 = 0
676 merid = lat1x .eq. -90 .or. slam12 == 0
692 csig1 = calp1 * cbet1
694 csig2 = calp2 * cbet2
697 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, 0d0),
698 + csig1 * csig2 + ssig1 * ssig2)
699 call lengs(n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
700 + cbet1, cbet2, s12x, m12x, dummy,
701 + scalp, mm12, mm21, ep2, c1a, c2a)
710 if (sig12 .lt. 1 .or. m12x .ge. 0)
then
713 a12x = sig12 / degree
721 if (.not. merid .and. sbet1 .eq. 0 .and.
722 + (f .le. 0 .or. lam12 .le. pi - f * pi))
then
732 m12x = b * sin(sig12)
738 else if (.not. merid)
then
743 sig12 = invsta(sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12,
744 + f, a3x, salp1, calp1, salp2, calp2, dnm, c1a, c2a)
746 if (sig12 .ge. 0)
then
748 s12x = sig12 * b * dnm
749 m12x = dnm**2 * b * sin(sig12 / dnm)
751 mm12 = cos(sig12 / dnm)
754 a12x = sig12 / degree
755 omg12 = lam12 / (f1 * dnm)
777 do 10 numit = 0, maxit2-1
780 v = lam12f(sbet1, cbet1, dn1, sbet2, cbet2, dn2,
781 + salp1, calp1, f, a3x, c3x, salp2, calp2, sig12,
782 + ssig1, csig1, ssig2, csig2,
783 + eps, omg12, numit .lt. maxit1, dv,
784 + c1a, c2a, c3a) - lam12
788 + .not. (abs(v) .ge. csmgt(8d0, 2d0, tripn) * tol0))
791 if (v .gt. 0 .and. (numit .gt. maxit1 .or.
792 + calp1/salp1 .gt. calp1b/salp1b))
then
795 else if (v .lt. 0 .and. (numit .gt. maxit1 .or.
796 + calp1/salp1 .lt. calp1a/salp1a))
then
800 if (numit .lt. maxit1 .and. dv .gt. 0)
then
804 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1
805 if (nsalp1 .gt. 0 .and. abs(dalp1) .lt. pi)
then
806 calp1 = calp1 * cdalp1 - salp1 * sdalp1
808 call norm(salp1, calp1)
812 tripn = abs(v) .le. 16 * tol0
824 salp1 = (salp1a + salp1b)/2
825 calp1 = (calp1a + calp1b)/2
826 call norm(salp1, calp1)
828 tripb = abs(salp1a - salp1) + (calp1a - calp1) .lt. tolb
829 + .or. abs(salp1 - salp1b) + (calp1 - calp1b) .lt. tolb
832 call lengs(eps, sig12, ssig1, csig1, dn1,
833 + ssig2, csig2, dn2, cbet1, cbet2, s12x, m12x, dummy,
834 + scalp, mm12, mm21, ep2, c1a, c2a)
837 a12x = sig12 / degree
838 omg12 = lam12 - omg12
844 if (redlp) m12 = 0 + m12x
848 salp0 = salp1 * cbet1
849 calp0 = hypotx(calp1, salp1 * sbet1)
850 if (calp0 .ne. 0 .and. salp0 .ne. 0)
then
853 csig1 = calp1 * cbet1
855 csig2 = calp2 * cbet2
857 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
859 a4 = a**2 * calp0 * salp0 * e2
860 call norm(ssig1, csig1)
861 call norm(ssig2, csig2)
862 call c4f(eps, c4x, c4a)
863 b41 = trgsum(.false., ssig1, csig1, c4a, nc4)
864 b42 = trgsum(.false., ssig2, csig2, c4a, nc4)
865 ss12 = a4 * (b42 - b41)
871 if (.not. merid .and. omg12 .lt. 0.75d0 * pi
872 + .and. sbet2 - sbet1 .lt. 1.75d0)
then
877 domg12 = 1 + cos(omg12)
880 alp12 = 2 * atan2(somg12 * (sbet1 * dbet2 + sbet2 * dbet1),
881 + domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) )
884 salp12 = salp2 * calp1 - calp2 * salp1
885 calp12 = calp2 * calp1 + salp2 * salp1
890 if (salp12 .eq. 0 .and. calp12 .lt. 0)
then
891 salp12 = tiny * calp1
894 alp12 = atan2(salp12, calp12)
896 ss12 = ss12 + c2 * alp12
897 ss12 = ss12 * swapp * lonsgn * latsgn
903 if (swapp .lt. 0)
then
904 call swap(salp1, salp2)
905 call swap(calp1, calp2)
906 if (scalp)
call swap(mm12, mm21)
909 salp1 = salp1 * swapp * lonsgn
910 calp1 = calp1 * swapp * latsgn
911 salp2 = salp2 * swapp * lonsgn
912 calp2 = calp2 * swapp * latsgn
915 azi1 = 0 - atan2(-salp1, calp1) / degree
916 azi2 = 0 - atan2(-salp2, calp2) / degree
942 subroutine area(a, f, lats, lons, n, AA, PP)
945 double precision a, f, lats(n), lons(n)
947 double precision AA, PP
949 integer i, omask, cross, trnsit
950 double precision s12, azi1, azi2, dummy, SS12, b, e2, c2, area0,
951 + atanhx, aacc(2), pacc(2)
953 double precision dblmin, dbleps, pi, degree, tiny,
954 + tol0, tol1, tol2, tolb, xthrsh
955 integer digits, maxit1, maxit2
957 common /geocom/ dblmin, dbleps, pi, degree, tiny,
958 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
965 call invers(a, f, lats(i+1), lons(i+1),
966 + lats(mod(i+1,n)+1), lons(mod(i+1,n)+1),
967 + s12, azi1, azi2, omask, dummy, dummy, dummy, dummy, ss12)
968 call accadd(pacc, s12)
969 call accadd(aacc, -ss12)
970 cross = cross + trnsit(lons(i+1), lons(mod(i+1,n)+1))
977 else if (e2 .gt. 0)
then
978 c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
980 c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
983 if (mod(abs(cross), 2) .eq. 1)
then
984 if (aacc(1) .lt. 0)
then
985 call accadd(aacc, +area0/2)
987 call accadd(aacc, -area0/2)
990 if (aacc(1) .gt. area0/2)
then
991 call accadd(aacc, -area0)
992 else if (aacc(1) .le. -area0/2)
then
993 call accadd(aacc, +area0)
1003 double precision dblmin, dbleps, pi, degree, tiny,
1004 + tol0, tol1, tol2, tolb, xthrsh
1005 integer digits, maxit1, maxit2
1008 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1009 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1013 double precision dblmin, dbleps, pi, degree, tiny,
1014 + tol0, tol1, tol2, tolb, xthrsh
1015 integer digits, maxit1, maxit2
1017 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1018 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1021 dblmin = 0.5d0**1022
1022 dbleps = 0.5d0**(digits-1)
1024 pi = atan2(0.0d0, -1.0d0)
1035 xthrsh = 1000 * tol2
1037 maxit2 = maxit1 + digits + 10
1044 subroutine lengs(eps, sig12,
1045 + ssig1, csig1, dn1, ssig2, csig2, dn2,
1046 + cbet1, cbet2, s12b, m12b, m0,
1047 + scalp, mm12, mm21, ep2, c1a, c2a)
1049 double precision eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1053 double precision s12b, m12b, m0
1055 double precision MM12, MM21
1057 double precision C1a(*), C2a(*)
1059 integer ord, nC1, nC2
1060 parameter(ord = 6, nc1 = ord, nc2 = ord)
1062 double precision A1m1f, A2m1f, TrgSum
1063 double precision A1m1, AB1, A2m1, AB2, J12, csig12, t
1071 ab1 = (1 + a1m1) * (trgsum(.true., ssig2, csig2, c1a, nc1) -
1072 + trgsum(.true., ssig1, csig1, c1a, nc1))
1074 ab2 = (1 + a2m1) * (trgsum(.true., ssig2, csig2, c2a, nc2) -
1075 + trgsum(.true., ssig1, csig1, c2a, nc2))
1077 j12 = m0 * sig12 + (ab1 - ab2)
1081 m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
1082 + csig1 * csig2 * j12
1084 s12b = (1 + a1m1) * sig12 + ab1
1086 csig12 = csig1 * csig2 + ssig1 * ssig2
1087 t = ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2)
1088 mm12 = csig12 + (t * ssig2 - csig2 * j12) * ssig1 / dn1
1089 mm21 = csig12 - (t * ssig1 - csig1 * j12) * ssig2 / dn2
1095 double precision function astrd(x, y)
1099 double precision x, y
1101 double precision cbrt, csmgt
1102 double precision k, p, q, r, S, r2, r3, disc, u,
1103 + t3, t, ang, v, uv, w
1108 if ( .not. (q .eq. 0 .and. r .lt. 0) )
then
1117 disc = s * (s + 2 * r3)
1119 if (disc .ge. 0)
then
1125 t3 = t3 + csmgt(-sqrt(disc), sqrt(disc), t3 .lt. 0)
1130 if (t .ne. 0) u = u + t + r2 / t
1133 ang = atan2(sqrt(-disc), -(s + r3))
1136 u = u + 2 * r * cos(ang / 3)
1142 uv = csmgt(q / (v - u), u + v, u .lt. 0)
1144 w = (uv - q) / (2 * v)
1148 k = uv / (sqrt(uv + w**2) + w)
1160 double precision function invsta(sbet1, cbet1, dn1,
1161 + sbet2, cbet2, dn2, lam12, f, a3x,
1162 + salp1, calp1, salp2, calp2, dnm,
1168 double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12,
1171 double precision salp1, calp1, salp2, calp2, dnm
1173 double precision C1a(*), C2a(*)
1175 double precision csmgt, hypotx, A3f, Astrd
1177 double precision f1, e2, ep2, n, etol2, k2, eps, sig12,
1178 + sbet12, cbet12, sbt12a, omg12, somg12, comg12, ssig12, csig12,
1179 + x, y, lamscl, betscl, cbt12a, bt12a, m12b, m0, dummy,
1182 double precision dblmin, dbleps, pi, degree, tiny,
1183 + tol0, tol1, tol2, tolb, xthrsh
1184 integer digits, maxit1, maxit2
1186 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1187 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1203 etol2 = 0.1d0 * tol2 /
1204 + sqrt( max(0.001d0, abs(f)) * min(1d0, 1 - f/2) / 2 )
1209 sbet12 = sbet2 * cbet1 - cbet2 * sbet1
1210 cbet12 = cbet2 * cbet1 + sbet2 * sbet1
1211 sbt12a = sbet2 * cbet1 + cbet2 * sbet1
1213 shortp = cbet12 .ge. 0 .and. sbet12 .lt. 0.5d0 .and.
1214 + cbet2 * lam12 .lt. 0.5d0
1218 sbetm2 = (sbet1 + sbet2)**2
1221 sbetm2 = sbetm2 / (sbetm2 + (cbet1 + cbet2)**2)
1222 dnm = sqrt(1 + ep2 * sbetm2)
1223 omg12 = omg12 / (f1 * dnm)
1228 salp1 = cbet2 * somg12
1229 calp1 = csmgt(sbet12 + cbet2 * sbet1 * somg12**2 / (1 + comg12),
1230 + sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12),
1233 ssig12 = hypotx(salp1, calp1)
1234 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12
1236 if (shortp .and. ssig12 .lt. etol2)
then
1238 salp2 = cbet1 * somg12
1239 calp2 = sbet12 - cbet1 * sbet2 *
1240 + csmgt(somg12**2 / (1 + comg12), 1 - comg12, comg12 .ge. 0)
1241 call norm(salp2, calp2)
1243 sig12 = atan2(ssig12, csig12)
1244 else if (abs(n) .gt. 0.1d0 .or. csig12 .ge. 0 .or.
1245 + ssig12 .ge. 6 * abs(n) * pi * cbet1**2)
then
1254 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
1255 lamscl = f * cbet1 * a3f(eps, a3x) * pi
1256 betscl = lamscl * cbet1
1257 x = (lam12 - pi) / lamscl
1261 cbt12a = cbet2 * cbet1 - sbet2 * sbet1
1262 bt12a = atan2(sbt12a, cbt12a)
1265 call lengs(n, pi + bt12a,
1266 + sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
1267 + cbet1, cbet2, dummy, m12b, m0, .false.,
1268 + dummy, dummy, ep2, c1a, c2a)
1269 x = -1 + m12b / (cbet1 * cbet2 * m0 * pi)
1270 betscl = csmgt(sbt12a / x, -f * cbet1**2 * pi,
1272 lamscl = betscl / cbet1
1273 y = (lam12 - pi) / lamscl
1276 if (y .gt. -tol1 .and. x .gt. -1 - xthrsh)
then
1279 salp1 = min(1d0, -x)
1280 calp1 = - sqrt(1 - salp1**2)
1282 calp1 = max(csmgt(0d0, 1d0, x .gt. -tol1), x)
1283 salp1 = sqrt(1 - calp1**2)
1322 + csmgt(-x * k/(1 + k), -y * (1 + k)/k, f .ge. 0)
1323 somg12 = sin(omg12a)
1324 comg12 = -cos(omg12a)
1326 salp1 = cbet2 * somg12
1327 calp1 = sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12)
1331 if (.not. (salp1 .le. 0))
then
1332 call norm(salp1, calp1)
1342 double precision function lam12f(sbet1, cbet1, dn1,
1343 + sbet2, cbet2, dn2, salp1, calp1, f, a3x, c3x, salp2, calp2,
1344 + sig12, ssig1, csig1, ssig2, csig2, eps, domg12, diffp, dlam12,
1347 double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2,
1348 + salp1, calp1, f, a3x(*), c3x(*)
1351 double precision salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
1354 double precision dlam12
1356 double precision C1a(*), C2a(*), C3a(*)
1359 parameter(ord = 6, nc3 = ord)
1361 double precision csmgt, hypotx, A3f, TrgSum
1363 double precision f1, e2, ep2, salp0, calp0,
1364 + somg1, comg1, somg2, comg2, omg12, lam12, b312, h0, k2, dummy
1366 double precision dblmin, dbleps, pi, degree, tiny,
1367 + tol0, tol1, tol2, tolb, xthrsh
1368 integer digits, maxit1, maxit2
1370 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1371 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1378 if (sbet1 .eq. 0 .and. calp1 .eq. 0) calp1 = -tiny
1381 salp0 = salp1 * cbet1
1383 calp0 = hypotx(calp1, salp1 * sbet1)
1388 somg1 = salp0 * sbet1
1389 csig1 = calp1 * cbet1
1391 call norm(ssig1, csig1)
1398 salp2 = csmgt(salp0 / cbet2, salp1, cbet2 .ne. cbet1)
1403 calp2 = csmgt(sqrt((calp1 * cbet1)**2 +
1404 + csmgt((cbet2 - cbet1) * (cbet1 + cbet2),
1405 + (sbet1 - sbet2) * (sbet1 + sbet2),
1406 + cbet1 .lt. -sbet1)) / cbet2,
1407 + abs(calp1), cbet2 .ne. cbet1 .or. abs(sbet2) .ne. -sbet1)
1411 somg2 = salp0 * sbet2
1412 csig2 = calp2 * cbet2
1414 call norm(ssig2, csig2)
1418 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, 0d0),
1419 + csig1 * csig2 + ssig1 * ssig2)
1422 omg12 = atan2(max(comg1 * somg2 - somg1 * comg2, 0d0),
1423 + comg1 * comg2 + somg1 * somg2)
1425 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
1426 call c3f(eps, c3x, c3a)
1427 b312 = (trgsum(.true., ssig2, csig2, c3a, nc3-1) -
1428 + trgsum(.true., ssig1, csig1, c3a, nc3-1))
1429 h0 = -f * a3f(eps, a3x)
1430 domg12 = salp0 * h0 * (sig12 + b312)
1431 lam12 = omg12 + domg12
1434 if (calp2 .eq. 0)
then
1435 dlam12 = - 2 * f1 * dn1 / sbet1
1437 call lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1438 + cbet1, cbet2, dummy, dlam12, dummy,
1439 + .false., dummy, dummy, ep2, c1a, c2a)
1440 dlam12 = dlam12 * f1 / (calp2 * cbet2)
1448 double precision function a3f(eps, A3x)
1450 integer ord, nA3, nA3x
1451 parameter(ord = 6, na3 = ord, na3x = na3)
1454 double precision eps
1456 double precision A3x(0: na3x-1)
1460 do 10 i = na3x-1, 0, -1
1461 a3f = eps * a3f + a3x(i)
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 t, mult
1482 do 20 k = nc3-1, 1 , -1
1484 do 10 i = nc3 - k, 1, -1
1486 t = eps * t + c3x(j)
1500 subroutine c4f(eps, C4x, c)
1503 integer ord, nC4, nC4x
1504 parameter(ord = 6, nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
1507 double precision eps, C4x(0:nc4x-1)
1509 double precision c(0:nc4-1)
1512 double precision t, mult
1515 do 20 k = nc4-1, 0, -1
1517 do 10 i = nc4 - k, 1, -1
1519 t = eps * t + c4x(j)
1535 double precision function a1m1f(eps)
1538 double precision eps
1540 double precision eps2, t
1543 t = eps2*(eps2*(eps2+4)+64)/256
1544 a1m1f = (t + eps) / (1 - eps)
1549 subroutine c1f(eps, c)
1552 parameter(ord = 6, nc1 = ord)
1555 double precision eps
1557 double precision c(nc1)
1559 double precision eps2, d
1563 c(1) = d*((6-eps2)*eps2-16)/32
1565 c(2) = d*((64-9*eps2)*eps2-128)/2048
1567 c(3) = d*(9*eps2-16)/768
1569 c(4) = d*(3*eps2-5)/512
1578 subroutine c1pf(eps, c)
1581 parameter(ord = 6, nc1p = ord)
1584 double precision eps
1586 double precision c(nc1p)
1588 double precision eps2, d
1592 c(1) = d*(eps2*(205*eps2-432)+768)/1536
1594 c(2) = d*(eps2*(4005*eps2-4736)+3840)/12288
1596 c(3) = d*(116-225*eps2)/384
1598 c(4) = d*(2695-7173*eps2)/7680
1602 c(6) = 38081*d/61440
1608 double precision function a2m1f(eps)
1610 double precision eps
1612 double precision eps2, t
1615 t = eps2*(eps2*(25*eps2+36)+64)/256
1616 a2m1f = t * (1 - eps) - eps
1621 subroutine c2f(eps, c)
1624 parameter(ord = 6, nc2 = ord)
1627 double precision eps
1629 double precision c(nc2)
1631 double precision eps2, d
1635 c(1) = d*(eps2*(eps2+2)+16)/32
1637 c(2) = d*(eps2*(35*eps2+64)+384)/2048
1639 c(3) = d*(15*eps2+80)/768
1641 c(4) = d*(7*eps2+35)/512
1650 subroutine a3cof(n, A3x)
1652 integer ord, nA3, nA3x
1653 parameter(ord = 6, na3 = ord, na3x = na3)
1658 double precision A3x(0:na3x-1)
1662 a3x(2) = (n*(3*n-1)-2)/8
1663 a3x(3) = ((-n-3)*n-1)/16
1664 a3x(4) = (-2*n-3)/64
1670 subroutine c3cof(n, C3x)
1672 integer ord, nC3, nC3x
1673 parameter(ord = 6, nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2)
1678 double precision C3x(0:nc3x-1)
1682 c3x(2) = ((3-n)*n+3)/64
1683 c3x(3) = (2*n+5)/128
1685 c3x(5) = ((n-3)*n+2)/32
1686 c3x(6) = ((-3*n-2)*n+3)/64
1689 c3x(9) = (n*(5*n-9)+5)/192
1690 c3x(10) = (9-10*n)/384
1692 c3x(12) = (7-14*n)/512
1701 subroutine c4cof(n, C4x)
1703 integer ord, nC4, nC4x
1704 parameter(ord = 6, nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
1709 double precision C4x(0:nc4x-1)
1711 c4x(0) = (n*(n*(n*(n*(100*n+208)+572)+3432)-12012)+30030)/45045
1712 c4x(1) = (n*(n*(n*(64*n+624)-4576)+6864)-3003)/15015
1713 c4x(2) = (n*((14144-10656*n)*n-4576)-858)/45045
1714 c4x(3) = ((-224*n-4784)*n+1573)/45045
1715 c4x(4) = (1088*n+156)/45045
1717 c4x(6) = (n*(n*((-64*n-624)*n+4576)-6864)+3003)/135135
1718 c4x(7) = (n*(n*(5952*n-11648)+9152)-2574)/135135
1719 c4x(8) = (n*(5792*n+1040)-1287)/135135
1720 c4x(9) = (468-2944*n)/135135
1722 c4x(11) = (n*((4160-1440*n)*n-4576)+1716)/225225
1723 c4x(12) = ((4992-8448*n)*n-1144)/225225
1724 c4x(13) = (1856*n-936)/225225
1726 c4x(15) = (n*(3584*n-3328)+1144)/315315
1727 c4x(16) = (1024*n-208)/105105
1728 c4x(17) = -136/63063d0
1729 c4x(18) = (832-2560*n)/405405
1730 c4x(19) = -128/135135d0
1731 c4x(20) = 128/99099d0
1736 double precision function sumx(u, v, t)
1738 double precision u, v
1742 double precision up, vpp
1753 double precision function angnm(x)
1757 if (x .ge. 180)
then
1759 else if (x .lt. -180)
then
1768 double precision function angnm2(x)
1772 double precision AngNm
1773 angnm2 = mod(x, 360d0)
1774 angnm2 = angnm(angnm2)
1779 double precision function angdif(x, y)
1785 double precision x, y
1787 double precision d, t, sumx
1789 if ((d - 180d0) + t .gt. 0d0)
then
1791 else if ((d + 180d0) + t .le. 0d0)
then
1799 double precision function angrnd(x)
1808 double precision y, z
1812 if (y .lt. z) y = z - (z - y)
1818 subroutine swap(x, y)
1820 double precision x, y
1830 double precision function hypotx(x, y)
1832 double precision x, y
1834 hypotx = sqrt(x**2 + y**2)
1839 subroutine norm(sinx, cosx)
1841 double precision sinx, cosx
1843 double precision hypotx, r
1844 r = hypotx(sinx, cosx)
1851 double precision function log1px(x)
1855 double precision csmgt, y, z
1858 log1px = csmgt(x, x * log(y) / z, z .eq. 0)
1863 double precision function atanhx(x)
1867 double precision log1px, y
1869 y = log1px(2 * y/(1 - y))/2
1875 double precision function cbrt(x)
1879 cbrt = sign(abs(x)**(1/3d0), x)
1884 double precision function csmgt(x, y, p)
1886 double precision x, y
1898 double precision function trgsum(sinp, sinx, cosx, c, n)
1907 double precision sinx, cosx, c(n)
1909 double precision ar, y0, y1
1913 ar = 2 * (cosx - sinx) * (cosx + sinx)
1915 if (mod(n, 2) .eq. 1)
then
1926 y1 = ar * y0 - y1 + c(k)
1927 y0 = ar * y1 - y0 + c(k-1)
1931 trgsum = 2 * sinx * cosx * y0
1934 trgsum = cosx * (y0 - y1)
1940 integer function trnsit(lon1, lon2)
1942 double precision lon1, lon2
1944 double precision lon1x, lon2x, lon12, AngNm, AngDif
1947 lon12 = angdif(lon1x, lon2x)
1949 if (lon1x .lt. 0 .and. lon2x .ge. 0 .and. lon12 .gt. 0)
then
1951 else if (lon2x .lt. 0 .and. lon1x .ge. 0 .and. lon12 .lt. 0)
then
1958 subroutine accini(s)
1961 double precision s(2)
1969 subroutine accadd(s, y)
1974 double precision s(2)
1976 double precision z, u, sumx
1977 z = sumx(y, s(2), u)
1978 s(1) = sumx(z, s(1), s(2))
1979 if (s(1) .eq. 0)
then