170 subroutine direct(a, f, lat1, lon1, azi1, s12a12, arcmod,
171 + lat2, lon2, azi2, omask, a12s12, m12, mm12, mm21, ss12)
173 double precision a, f, lat1, lon1, azi1, s12a12
177 double precision lat2, lon2, azi2
179 double precision a12s12, m12, MM12, MM21, SS12
181 integer ord, nC1, nC1p, nC2, nA3, nA3x, nC3, nC3x, nC4, nC4x
182 parameter(ord = 6, nc1 = ord, nc1p = ord,
183 + nc2 = ord, na3 = ord, na3x = na3,
184 + nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2,
185 + nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
186 double precision A3x(0:na3x-1), C3x(0:nc3x-1), C4x(0:nc4x-1),
187 + c1a(nc1), c1pa(nc1p), c2a(nc2), c3a(nc3-1), c4a(0:nc4-1)
189 double precision csmgt, atanhx, hypotx,
190 + angnm, angnm2, angrnd, trgsum, a1m1f, a2m1f, a3f
191 logical arcp, redlp, scalp, areap
192 double precision e2, f1, ep2, n, b, c2,
193 + lon1x, azi1x, phi, alp1, salp0, calp0, k2, eps,
194 + salp1, calp1, ssig1, csig1, cbet1, sbet1, dn1, somg1, comg1,
195 + salp2, calp2, ssig2, csig2, sbet2, cbet2, dn2, somg2, comg2,
196 + ssig12, csig12, salp12, calp12, omg12, lam12, lon12,
197 + sig12, stau1, ctau1, tau12, s12a, t, s, c, serr,
198 + a1m1, a2m1, a3c, a4, ab1, ab2,
199 + b11, b12, b21, b22, b31, b41, b42, j12
201 double precision dblmin, dbleps, pi, degree, tiny,
202 + tol0, tol1, tol2, tolb, xthrsh
203 integer digits, maxit1, maxit2
205 common /geocom/ dblmin, dbleps, pi, degree, tiny,
206 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
208 if (.not.init)
call geoini
217 arcp = mod(omask/1, 2) == 1
218 redlp = mod(omask/2, 2) == 1
219 scalp = mod(omask/4, 2) == 1
220 areap = mod(omask/8, 2) == 1
225 else if (e2 .gt. 0)
then
226 c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
228 c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
234 if (areap)
call c4cof(n, c4x)
237 azi1x = angrnd(angnm(azi1))
241 alp1 = azi1x * degree
244 salp1 = csmgt(0d0, sin(alp1), azi1x .eq. -180)
245 calp1 = csmgt(0d0, cos(alp1), abs(azi1x) .eq. 90)
249 sbet1 = f1 * sin(phi)
250 cbet1 = csmgt(tiny, cos(phi), abs(lat1) .eq. 90)
251 call norm(sbet1, cbet1)
252 dn1 = sqrt(1 + ep2 * sbet1**2)
256 salp0 = salp1 * cbet1
259 calp0 = hypotx(calp1, salp1 * sbet1)
270 somg1 = salp0 * sbet1
271 csig1 = csmgt(cbet1 * calp1, 1d0, sbet1 .ne. 0 .or. calp1 .ne. 0)
274 call norm(ssig1, csig1)
278 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
282 b11 = trgsum(.true., ssig1, csig1, c1a, nc1)
286 stau1 = ssig1 * c + csig1 * s
287 ctau1 = csig1 * c - ssig1 * s
291 if (.not. arcmod)
call c1pf(eps, c1pa)
293 if (redlp .or. scalp)
then
296 b21 = trgsum(.true., ssig1, csig1, c2a, nc2)
303 call c3f(eps, c3x, c3a)
304 a3c = -f * salp0 * a3f(eps, a3x)
305 b31 = trgsum(.true., ssig1, csig1, c3a, nc3-1)
308 call c4f(eps, c4x, c4a)
310 a4 = a**2 * calp0 * salp0 * e2
311 b41 = trgsum(.false., ssig1, csig1, c4a, nc4)
320 sig12 = s12a12 * degree
322 s12a = s12a - 180 * aint(s12a / 180)
323 ssig12 = csmgt(0d0, sin(sig12), s12a .eq. 0)
324 csig12 = csmgt(0d0, cos(sig12), s12a .eq. 90)
329 tau12 = s12a12 / (b * (1 + a1m1))
333 b12 = - trgsum(.true.,
334 + stau1 * c + ctau1 * s, ctau1 * c - stau1 * s, c1pa, nc1p)
335 sig12 = tau12 - (b12 - b11)
338 if (abs(f) .gt. 0.01d0)
then
360 ssig2 = ssig1 * csig12 + csig1 * ssig12
361 csig2 = csig1 * csig12 - ssig1 * ssig12
362 b12 = trgsum(.true., ssig2, csig2, c1a, nc1)
363 serr = (1 + a1m1) * (sig12 + (b12 - b11)) - s12a12 / b
364 sig12 = sig12 - serr / sqrt(1 + k2 * ssig2**2)
372 ssig2 = ssig1 * csig12 + csig1 * ssig12
373 csig2 = csig1 * csig12 - ssig1 * ssig12
374 dn2 = sqrt(1 + k2 * ssig2**2)
375 if (arcmod .or. abs(f) .gt. 0.01d0)
376 + b12 = trgsum(.true., ssig2, csig2, c1a, nc1)
377 ab1 = (1 + a1m1) * (b12 - b11)
380 sbet2 = calp0 * ssig2
382 cbet2 = hypotx(salp0, calp0 * csig2)
383 if (cbet2 .eq. 0)
then
390 somg2 = salp0 * ssig2
395 calp2 = calp0 * csig2
397 omg12 = atan2(somg2 * comg1 - comg2 * somg1,
398 + comg2 * comg1 + somg2 * somg1)
400 lam12 = omg12 + a3c *
401 + ( sig12 + (trgsum(.true., ssig2, csig2, c3a, nc3-1)
403 lon12 = lam12 / degree
406 lon12 = angnm2(lon12)
407 lon2 = angnm(lon1x + lon12)
408 lat2 = atan2(sbet2, f1 * cbet2) / degree
410 azi2 = 0 - atan2(-salp2, calp2) / degree
412 if (redlp .or. scalp)
then
413 b22 = trgsum(.true., ssig2, csig2, c2a, nc2)
414 ab2 = (1 + a2m1) * (b22 - b21)
415 j12 = (a1m1 - a2m1) * sig12 + (ab1 - ab2)
419 if (redlp) m12 = b * ((dn2 * (csig1 * ssig2) -
420 + dn1 * (ssig1 * csig2)) - csig1 * csig2 * j12)
422 t = k2 * (ssig2 - ssig1) * (ssig2 + ssig1) / (dn1 + dn2)
423 mm12 = csig12 + (t * ssig2 - csig2 * j12) * ssig1 / dn1
424 mm21 = csig12 - (t * ssig1 - csig1 * j12) * ssig2 / dn2
428 b42 = trgsum(.false., ssig2, csig2, c4a, nc4)
429 if (calp0 .eq. 0 .or. salp0 .eq. 0)
then
431 salp12 = salp2 * calp1 - calp2 * salp1
432 calp12 = calp2 * calp1 + salp2 * salp1
436 if (salp12 .eq. 0 .and. calp12 .lt. 0)
then
437 salp12 = tiny * calp1
449 salp12 = calp0 * salp0 *
450 + csmgt(csig1 * (1 - csig12) + ssig12 * ssig1,
451 + ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1),
453 calp12 = salp0**2 + calp0**2 * csig1 * csig2
455 ss12 = c2 * atan2(salp12, calp12) + a4 * (b42 - b41)
458 if (arcp) a12s12 = csmgt(b * ((1 + a1m1) * sig12 + ab1),
459 + sig12 / degree, arcmod)
507 subroutine invers(a, f, lat1, lon1, lat2, lon2,
508 + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
510 double precision a, f, lat1, lon1, lat2, lon2
513 double precision s12, azi1, azi2
515 double precision a12, m12, MM12, MM21, SS12
517 integer ord, nC1, nC2, nA3, nA3x, nC3, nC3x, nC4, nC4x
518 parameter(ord = 6, nc1 = ord, nc2 = ord, na3 = ord, na3x = na3,
519 + nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2,
520 + nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
521 double precision A3x(0:na3x-1), C3x(0:nc3x-1), C4x(0:nc4x-1),
522 + c1a(nc1), c2a(nc2), c3a(nc3-1), c4a(0:nc4-1)
524 double precision csmgt, atanhx, hypotx,
525 + angnm, angdif, angrnd, trgsum, lam12f, invsta
526 integer latsgn, lonsgn, swapp, numit
527 logical arcp, redlp, scalp, areap, merid, tripn, tripb
529 double precision e2, f1, ep2, n, b, c2,
530 + lat1x, lat2x, phi, salp0, calp0, k2, eps,
531 + salp1, calp1, ssig1, csig1, cbet1, sbet1, dbet1, dn1,
532 + salp2, calp2, ssig2, csig2, sbet2, cbet2, dbet2, dn2,
533 + slam12, clam12, salp12, calp12, omg12, lam12, lon12,
534 + salp1a, calp1a, salp1b, calp1b,
535 + dalp1, sdalp1, cdalp1, nsalp1, alp12, somg12, domg12,
536 + sig12, v, dv, dnm, dummy,
537 + a4, b41, b42, s12x, m12x, a12x
539 double precision dblmin, dbleps, pi, degree, tiny,
540 + tol0, tol1, tol2, tolb, xthrsh
541 integer digits, maxit1, maxit2
543 common /geocom/ dblmin, dbleps, pi, degree, tiny,
544 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
546 if (.not.init)
call geoini
555 arcp = mod(omask/1, 2) == 1
556 redlp = mod(omask/2, 2) == 1
557 scalp = mod(omask/4, 2) == 1
558 areap = mod(omask/8, 2) == 1
563 else if (e2 .gt. 0)
then
564 c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
566 c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
572 if (areap)
call c4cof(n, c4x)
577 lon12 = angdif(angnm(lon1), angnm(lon2))
579 lon12 = angrnd(lon12)
581 if (lon12 .ge. 0)
then
586 lon12 = lon12 * lonsgn
591 if (abs(lat1x) .ge. abs(lat2x))
then
596 if (swapp .lt. 0)
then
598 call swap(lat1x, lat2x)
601 if (lat1x .lt. 0)
then
606 lat1x = lat1x * latsgn
607 lat2x = lat2x * latsgn
622 sbet1 = f1 * sin(phi)
623 cbet1 = csmgt(tiny, cos(phi), lat1x .eq. -90)
624 call norm(sbet1, cbet1)
628 sbet2 = f1 * sin(phi)
629 cbet2 = csmgt(tiny, cos(phi), abs(lat2x) .eq. 90)
630 call norm(sbet2, cbet2)
641 if (cbet1 .lt. -sbet1)
then
642 if (cbet2 .eq. cbet1) sbet2 = sign(sbet1, sbet2)
644 if (abs(sbet2) .eq. -sbet1) cbet2 = cbet1
647 dn1 = sqrt(1 + ep2 * sbet1**2)
648 dn2 = sqrt(1 + ep2 * sbet2**2)
650 lam12 = lon12 * degree
652 if (lon12 .eq. 180) slam12 = 0
658 merid = lat1x .eq. -90 .or. slam12 == 0
674 csig1 = calp1 * cbet1
676 csig2 = calp2 * cbet2
679 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, 0d0),
680 + csig1 * csig2 + ssig1 * ssig2)
681 call lengs(n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
682 + cbet1, cbet2, s12x, m12x, dummy,
683 + scalp, mm12, mm21, ep2, c1a, c2a)
692 if (sig12 .lt. 1 .or. m12x .ge. 0)
then
695 a12x = sig12 / degree
703 if (.not. merid .and. sbet1 .eq. 0 .and.
704 + (f .le. 0 .or. lam12 .le. pi - f * pi))
then
714 m12x = b * sin(sig12)
720 else if (.not. merid)
then
725 sig12 = invsta(sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12,
726 + f, a3x, salp1, calp1, salp2, calp2, dnm, c1a, c2a)
728 if (sig12 .ge. 0)
then
730 s12x = sig12 * b * dnm
731 m12x = dnm**2 * b * sin(sig12 / dnm)
733 mm12 = cos(sig12 / dnm)
736 a12x = sig12 / degree
737 omg12 = lam12 / (f1 * dnm)
759 do 10 numit = 0, maxit2-1
762 v = lam12f(sbet1, cbet1, dn1, sbet2, cbet2, dn2,
763 + salp1, calp1, f, a3x, c3x, salp2, calp2, sig12,
764 + ssig1, csig1, ssig2, csig2,
765 + eps, omg12, numit .lt. maxit1, dv,
766 + c1a, c2a, c3a) - lam12
770 + .not. (abs(v) .ge. csmgt(8d0, 2d0, tripn) * tol0))
773 if (v .gt. 0 .and. (numit .gt. maxit1 .or.
774 + calp1/salp1 .gt. calp1b/salp1b))
then
777 else if (v .lt. 0 .and. (numit .gt. maxit1 .or.
778 + calp1/salp1 .lt. calp1a/salp1a))
then
782 if (numit .lt. maxit1 .and. dv .gt. 0)
then
786 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1
787 if (nsalp1 .gt. 0 .and. abs(dalp1) .lt. pi)
then
788 calp1 = calp1 * cdalp1 - salp1 * sdalp1
790 call norm(salp1, calp1)
794 tripn = abs(v) .le. 16 * tol0
806 salp1 = (salp1a + salp1b)/2
807 calp1 = (calp1a + calp1b)/2
808 call norm(salp1, calp1)
810 tripb = abs(salp1a - salp1) + (calp1a - calp1) .lt. tolb
811 + .or. abs(salp1 - salp1b) + (calp1 - calp1b) .lt. tolb
814 call lengs(eps, sig12, ssig1, csig1, dn1,
815 + ssig2, csig2, dn2, cbet1, cbet2, s12x, m12x, dummy,
816 + scalp, mm12, mm21, ep2, c1a, c2a)
819 a12x = sig12 / degree
820 omg12 = lam12 - omg12
826 if (redlp) m12 = 0 + m12x
830 salp0 = salp1 * cbet1
831 calp0 = hypotx(calp1, salp1 * sbet1)
832 if (calp0 .ne. 0 .and. salp0 .ne. 0)
then
835 csig1 = calp1 * cbet1
837 csig2 = calp2 * cbet2
839 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
841 a4 = a**2 * calp0 * salp0 * e2
842 call norm(ssig1, csig1)
843 call norm(ssig2, csig2)
844 call c4f(eps, c4x, c4a)
845 b41 = trgsum(.false., ssig1, csig1, c4a, nc4)
846 b42 = trgsum(.false., ssig2, csig2, c4a, nc4)
847 ss12 = a4 * (b42 - b41)
853 if (.not. merid .and. omg12 .lt. 0.75d0 * pi
854 + .and. sbet2 - sbet1 .lt. 1.75d0)
then
859 domg12 = 1 + cos(omg12)
862 alp12 = 2 * atan2(somg12 * (sbet1 * dbet2 + sbet2 * dbet1),
863 + domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) )
866 salp12 = salp2 * calp1 - calp2 * salp1
867 calp12 = calp2 * calp1 + salp2 * salp1
872 if (salp12 .eq. 0 .and. calp12 .lt. 0)
then
873 salp12 = tiny * calp1
876 alp12 = atan2(salp12, calp12)
878 ss12 = ss12 + c2 * alp12
879 ss12 = ss12 * swapp * lonsgn * latsgn
885 if (swapp .lt. 0)
then
886 call swap(salp1, salp2)
887 call swap(calp1, calp2)
888 if (scalp)
call swap(mm12, mm21)
891 salp1 = salp1 * swapp * lonsgn
892 calp1 = calp1 * swapp * latsgn
893 salp2 = salp2 * swapp * lonsgn
894 calp2 = calp2 * swapp * latsgn
897 azi1 = 0 - atan2(-salp1, calp1) / degree
898 azi2 = 0 - atan2(-salp2, calp2) / degree
924 subroutine area(a, f, lats, lons, n, AA, PP)
927 double precision a, f, lats(n), lons(n)
929 double precision AA, PP
931 integer i, omask, cross, trnsit
932 double precision s12, azi1, azi2, dummy, SS12, b, e2, c2, area0,
933 + atanhx, aacc(2), pacc(2)
935 double precision dblmin, dbleps, pi, degree, tiny,
936 + tol0, tol1, tol2, tolb, xthrsh
937 integer digits, maxit1, maxit2
939 common /geocom/ dblmin, dbleps, pi, degree, tiny,
940 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
947 call invers(a, f, lats(i+1), lons(i+1),
948 + lats(mod(i+1,n)+1), lons(mod(i+1,n)+1),
949 + s12, azi1, azi2, omask, dummy, dummy, dummy, dummy, ss12)
950 call accadd(pacc, s12)
951 call accadd(aacc, -ss12)
952 cross = cross + trnsit(lons(i+1), lons(mod(i+1,n)+1))
959 else if (e2 .gt. 0)
then
960 c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
962 c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
965 if (mod(abs(cross), 2) .eq. 1)
then
966 if (aacc(1) .lt. 0)
then
967 call accadd(aacc, +area0/2)
969 call accadd(aacc, -area0/2)
972 if (aacc(1) .gt. area0/2)
then
973 call accadd(aacc, -area0)
974 else if (aacc(1) .le. -area0/2)
then
975 call accadd(aacc, +area0)
985 double precision dblmin, dbleps, pi, degree, tiny,
986 + tol0, tol1, tol2, tolb, xthrsh
987 integer digits, maxit1, maxit2
990 common /geocom/ dblmin, dbleps, pi, degree, tiny,
991 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
995 double precision dblmin, dbleps, pi, degree, tiny,
996 + tol0, tol1, tol2, tolb, xthrsh
997 integer digits, maxit1, maxit2
999 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1000 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1003 dblmin = 0.5d0**1022
1004 dbleps = 0.5d0**(digits-1)
1006 pi = atan2(0.0d0, -1.0d0)
1017 xthrsh = 1000 * tol2
1019 maxit2 = maxit1 + digits + 10
1026 subroutine lengs(eps, sig12,
1027 + ssig1, csig1, dn1, ssig2, csig2, dn2,
1028 + cbet1, cbet2, s12b, m12b, m0,
1029 + scalp, mm12, mm21, ep2, c1a, c2a)
1031 double precision eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1035 double precision s12b, m12b, m0
1037 double precision MM12, MM21
1039 double precision C1a(*), C2a(*)
1041 integer ord, nC1, nC2
1042 parameter(ord = 6, nc1 = ord, nc2 = ord)
1044 double precision A1m1f, A2m1f, TrgSum
1045 double precision A1m1, AB1, A2m1, AB2, J12, csig12, t
1053 ab1 = (1 + a1m1) * (trgsum(.true., ssig2, csig2, c1a, nc1) -
1054 + trgsum(.true., ssig1, csig1, c1a, nc1))
1056 ab2 = (1 + a2m1) * (trgsum(.true., ssig2, csig2, c2a, nc2) -
1057 + trgsum(.true., ssig1, csig1, c2a, nc2))
1059 j12 = m0 * sig12 + (ab1 - ab2)
1063 m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
1064 + csig1 * csig2 * j12
1066 s12b = (1 + a1m1) * sig12 + ab1
1068 csig12 = csig1 * csig2 + ssig1 * ssig2
1069 t = ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2)
1070 mm12 = csig12 + (t * ssig2 - csig2 * j12) * ssig1 / dn1
1071 mm21 = csig12 - (t * ssig1 - csig1 * j12) * ssig2 / dn2
1077 double precision function astrd(x, y)
1081 double precision x, y
1083 double precision cbrt, csmgt
1084 double precision k, p, q, r, S, r2, r3, disc, u,
1085 + t3, t, ang, v, uv, w
1090 if ( .not. (q .eq. 0 .and. r .lt. 0) )
then
1099 disc = s * (s + 2 * r3)
1101 if (disc .ge. 0)
then
1107 t3 = t3 + csmgt(-sqrt(disc), sqrt(disc), t3 .lt. 0)
1112 if (t .ne. 0) u = u + t + r2 / t
1115 ang = atan2(sqrt(-disc), -(s + r3))
1118 u = u + 2 * r * cos(ang / 3)
1124 uv = csmgt(q / (v - u), u + v, u .lt. 0)
1126 w = (uv - q) / (2 * v)
1130 k = uv / (sqrt(uv + w**2) + w)
1142 double precision function invsta(sbet1, cbet1, dn1,
1143 + sbet2, cbet2, dn2, lam12, f, a3x,
1144 + salp1, calp1, salp2, calp2, dnm,
1150 double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12,
1153 double precision salp1, calp1, salp2, calp2, dnm
1155 double precision C1a(*), C2a(*)
1157 double precision csmgt, hypotx, A3f, Astrd
1159 double precision f1, e2, ep2, n, etol2, k2, eps, sig12,
1160 + sbet12, cbet12, sbt12a, omg12, somg12, comg12, ssig12, csig12,
1161 + x, y, lamscl, betscl, cbt12a, bt12a, m12b, m0, dummy,
1164 double precision dblmin, dbleps, pi, degree, tiny,
1165 + tol0, tol1, tol2, tolb, xthrsh
1166 integer digits, maxit1, maxit2
1168 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1169 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1185 etol2 = 0.1d0 * tol2 /
1186 + sqrt( max(0.001d0, abs(f)) * min(1d0, 1 - f/2) / 2 )
1191 sbet12 = sbet2 * cbet1 - cbet2 * sbet1
1192 cbet12 = cbet2 * cbet1 + sbet2 * sbet1
1193 sbt12a = sbet2 * cbet1 + cbet2 * sbet1
1195 shortp = cbet12 .ge. 0 .and. sbet12 .lt. 0.5d0 .and.
1196 + cbet2 * lam12 .lt. 0.5d0
1200 sbetm2 = (sbet1 + sbet2)**2
1203 sbetm2 = sbetm2 / (sbetm2 + (cbet1 + cbet2)**2)
1204 dnm = sqrt(1 + ep2 * sbetm2)
1205 omg12 = omg12 / (f1 * dnm)
1210 salp1 = cbet2 * somg12
1211 calp1 = csmgt(sbet12 + cbet2 * sbet1 * somg12**2 / (1 + comg12),
1212 + sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12),
1215 ssig12 = hypotx(salp1, calp1)
1216 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12
1218 if (shortp .and. ssig12 .lt. etol2)
then
1220 salp2 = cbet1 * somg12
1221 calp2 = sbet12 - cbet1 * sbet2 *
1222 + csmgt(somg12**2 / (1 + comg12), 1 - comg12, comg12 .ge. 0)
1223 call norm(salp2, calp2)
1225 sig12 = atan2(ssig12, csig12)
1226 else if (abs(n) .gt. 0.1d0 .or. csig12 .ge. 0 .or.
1227 + ssig12 .ge. 6 * abs(n) * pi * cbet1**2)
then
1236 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
1237 lamscl = f * cbet1 * a3f(eps, a3x) * pi
1238 betscl = lamscl * cbet1
1239 x = (lam12 - pi) / lamscl
1243 cbt12a = cbet2 * cbet1 - sbet2 * sbet1
1244 bt12a = atan2(sbt12a, cbt12a)
1247 call lengs(n, pi + bt12a,
1248 + sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
1249 + cbet1, cbet2, dummy, m12b, m0, .false.,
1250 + dummy, dummy, ep2, c1a, c2a)
1251 x = -1 + m12b / (cbet1 * cbet2 * m0 * pi)
1252 betscl = csmgt(sbt12a / x, -f * cbet1**2 * pi,
1254 lamscl = betscl / cbet1
1255 y = (lam12 - pi) / lamscl
1258 if (y .gt. -tol1 .and. x .gt. -1 - xthrsh)
then
1261 salp1 = min(1d0, -x)
1262 calp1 = - sqrt(1 - salp1**2)
1264 calp1 = max(csmgt(0d0, 1d0, x .gt. -tol1), x)
1265 salp1 = sqrt(1 - calp1**2)
1304 + csmgt(-x * k/(1 + k), -y * (1 + k)/k, f .ge. 0)
1305 somg12 = sin(omg12a)
1306 comg12 = -cos(omg12a)
1308 salp1 = cbet2 * somg12
1309 calp1 = sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12)
1313 if (salp1 .gt. 0)
then
1314 call norm(salp1, calp1)
1324 double precision function lam12f(sbet1, cbet1, dn1,
1325 + sbet2, cbet2, dn2, salp1, calp1, f, a3x, c3x, salp2, calp2,
1326 + sig12, ssig1, csig1, ssig2, csig2, eps, domg12, diffp, dlam12,
1329 double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2,
1330 + salp1, calp1, f, a3x(*), c3x(*)
1333 double precision salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
1336 double precision dlam12
1338 double precision C1a(*), C2a(*), C3a(*)
1341 parameter(ord = 6, nc3 = ord)
1343 double precision csmgt, hypotx, A3f, TrgSum
1345 double precision f1, e2, ep2, salp0, calp0,
1346 + somg1, comg1, somg2, comg2, omg12, lam12, b312, h0, k2, dummy
1348 double precision dblmin, dbleps, pi, degree, tiny,
1349 + tol0, tol1, tol2, tolb, xthrsh
1350 integer digits, maxit1, maxit2
1352 common /geocom/ dblmin, dbleps, pi, degree, tiny,
1353 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1360 if (sbet1 .eq. 0 .and. calp1 .eq. 0) calp1 = -tiny
1363 salp0 = salp1 * cbet1
1365 calp0 = hypotx(calp1, salp1 * sbet1)
1370 somg1 = salp0 * sbet1
1371 csig1 = calp1 * cbet1
1373 call norm(ssig1, csig1)
1380 salp2 = csmgt(salp0 / cbet2, salp1, cbet2 .ne. cbet1)
1385 calp2 = csmgt(sqrt((calp1 * cbet1)**2 +
1386 + csmgt((cbet2 - cbet1) * (cbet1 + cbet2),
1387 + (sbet1 - sbet2) * (sbet1 + sbet2),
1388 + cbet1 .lt. -sbet1)) / cbet2,
1389 + abs(calp1), cbet2 .ne. cbet1 .or. abs(sbet2) .ne. -sbet1)
1393 somg2 = salp0 * sbet2
1394 csig2 = calp2 * cbet2
1396 call norm(ssig2, csig2)
1400 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, 0d0),
1401 + csig1 * csig2 + ssig1 * ssig2)
1404 omg12 = atan2(max(comg1 * somg2 - somg1 * comg2, 0d0),
1405 + comg1 * comg2 + somg1 * somg2)
1407 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
1408 call c3f(eps, c3x, c3a)
1409 b312 = (trgsum(.true., ssig2, csig2, c3a, nc3-1) -
1410 + trgsum(.true., ssig1, csig1, c3a, nc3-1))
1411 h0 = -f * a3f(eps, a3x)
1412 domg12 = salp0 * h0 * (sig12 + b312)
1413 lam12 = omg12 + domg12
1416 if (calp2 .eq. 0)
then
1417 dlam12 = - 2 * f1 * dn1 / sbet1
1419 call lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1420 + cbet1, cbet2, dummy, dlam12, dummy,
1421 + .false., dummy, dummy, ep2, c1a, c2a)
1422 dlam12 = dlam12 * f1 / (calp2 * cbet2)
1430 double precision function a3f(eps, A3x)
1432 integer ord, nA3, nA3x
1433 parameter(ord = 6, na3 = ord, na3x = na3)
1436 double precision eps
1438 double precision A3x(0: na3x-1)
1442 do 10 i = na3x-1, 0, -1
1443 a3f = eps * a3f + a3x(i)
1449 subroutine c3f(eps, C3x, c)
1452 integer ord, nC3, nC3x
1453 parameter(ord = 6, nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2)
1456 double precision eps, C3x(0:nc3x-1)
1458 double precision c(nc3-1)
1461 double precision t, mult
1464 do 20 k = nc3-1, 1 , -1
1466 do 10 i = nc3 - k, 1, -1
1468 t = eps * t + c3x(j)
1482 subroutine c4f(eps, C4x, c)
1485 integer ord, nC4, nC4x
1486 parameter(ord = 6, nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
1489 double precision eps, C4x(0:nc4x-1)
1491 double precision c(0:nc4-1)
1494 double precision t, mult
1497 do 20 k = nc4-1, 0, -1
1499 do 10 i = nc4 - k, 1, -1
1501 t = eps * t + c4x(j)
1517 double precision function a1m1f(eps)
1520 double precision eps
1522 double precision eps2, t
1525 t = eps2*(eps2*(eps2+4)+64)/256
1526 a1m1f = (t + eps) / (1 - eps)
1531 subroutine c1f(eps, c)
1534 parameter(ord = 6, nc1 = ord)
1537 double precision eps
1539 double precision c(nc1)
1541 double precision eps2, d
1545 c(1) = d*((6-eps2)*eps2-16)/32
1547 c(2) = d*((64-9*eps2)*eps2-128)/2048
1549 c(3) = d*(9*eps2-16)/768
1551 c(4) = d*(3*eps2-5)/512
1560 subroutine c1pf(eps, c)
1563 parameter(ord = 6, nc1p = ord)
1566 double precision eps
1568 double precision c(nc1p)
1570 double precision eps2, d
1574 c(1) = d*(eps2*(205*eps2-432)+768)/1536
1576 c(2) = d*(eps2*(4005*eps2-4736)+3840)/12288
1578 c(3) = d*(116-225*eps2)/384
1580 c(4) = d*(2695-7173*eps2)/7680
1584 c(6) = 38081*d/61440
1590 double precision function a2m1f(eps)
1592 double precision eps
1594 double precision eps2, t
1597 t = eps2*(eps2*(25*eps2+36)+64)/256
1598 a2m1f = t * (1 - eps) - eps
1603 subroutine c2f(eps, c)
1606 parameter(ord = 6, nc2 = ord)
1609 double precision eps
1611 double precision c(nc2)
1613 double precision eps2, d
1617 c(1) = d*(eps2*(eps2+2)+16)/32
1619 c(2) = d*(eps2*(35*eps2+64)+384)/2048
1621 c(3) = d*(15*eps2+80)/768
1623 c(4) = d*(7*eps2+35)/512
1632 subroutine a3cof(n, A3x)
1634 integer ord, nA3, nA3x
1635 parameter(ord = 6, na3 = ord, na3x = na3)
1640 double precision A3x(0:na3x-1)
1644 a3x(2) = (n*(3*n-1)-2)/8
1645 a3x(3) = ((-n-3)*n-1)/16
1646 a3x(4) = (-2*n-3)/64
1652 subroutine c3cof(n, C3x)
1654 integer ord, nC3, nC3x
1655 parameter(ord = 6, nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2)
1660 double precision C3x(0:nc3x-1)
1664 c3x(2) = ((3-n)*n+3)/64
1665 c3x(3) = (2*n+5)/128
1667 c3x(5) = ((n-3)*n+2)/32
1668 c3x(6) = ((-3*n-2)*n+3)/64
1671 c3x(9) = (n*(5*n-9)+5)/192
1672 c3x(10) = (9-10*n)/384
1674 c3x(12) = (7-14*n)/512
1683 subroutine c4cof(n, C4x)
1685 integer ord, nC4, nC4x
1686 parameter(ord = 6, nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
1691 double precision C4x(0:nc4x-1)
1693 c4x(0) = (n*(n*(n*(n*(100*n+208)+572)+3432)-12012)+30030)/45045
1694 c4x(1) = (n*(n*(n*(64*n+624)-4576)+6864)-3003)/15015
1695 c4x(2) = (n*((14144-10656*n)*n-4576)-858)/45045
1696 c4x(3) = ((-224*n-4784)*n+1573)/45045
1697 c4x(4) = (1088*n+156)/45045
1699 c4x(6) = (n*(n*((-64*n-624)*n+4576)-6864)+3003)/135135
1700 c4x(7) = (n*(n*(5952*n-11648)+9152)-2574)/135135
1701 c4x(8) = (n*(5792*n+1040)-1287)/135135
1702 c4x(9) = (468-2944*n)/135135
1704 c4x(11) = (n*((4160-1440*n)*n-4576)+1716)/225225
1705 c4x(12) = ((4992-8448*n)*n-1144)/225225
1706 c4x(13) = (1856*n-936)/225225
1708 c4x(15) = (n*(3584*n-3328)+1144)/315315
1709 c4x(16) = (1024*n-208)/105105
1710 c4x(17) = -136/63063d0
1711 c4x(18) = (832-2560*n)/405405
1712 c4x(19) = -128/135135d0
1713 c4x(20) = 128/99099d0
1718 double precision function sumx(u, v, t)
1720 double precision u, v
1724 double precision up, vpp
1735 double precision function angnm(x)
1739 if (x .ge. 180)
then
1741 else if (x .lt. -180)
then
1750 double precision function angnm2(x)
1754 double precision AngNm
1755 angnm2 = mod(x, 360d0)
1756 angnm2 = angnm(angnm2)
1761 double precision function angdif(x, y)
1767 double precision x, y
1769 double precision d, t, sumx
1771 if ((d - 180d0) + t .gt. 0d0)
then
1773 else if ((d + 180d0) + t .le. 0d0)
then
1781 double precision function angrnd(x)
1790 double precision y, z
1794 if (y .lt. z) y = z - (z - y)
1800 subroutine swap(x, y)
1802 double precision x, y
1812 double precision function hypotx(x, y)
1814 double precision x, y
1816 hypotx = sqrt(x**2 + y**2)
1821 subroutine norm(sinx, cosx)
1823 double precision sinx, cosx
1825 double precision hypotx, r
1826 r = hypotx(sinx, cosx)
1833 double precision function log1px(x)
1837 double precision csmgt, y, z
1840 log1px = csmgt(x, x * log(y) / z, z .eq. 0)
1845 double precision function atanhx(x)
1849 double precision log1px, y
1851 y = log1px(2 * y/(1 - y))/2
1857 double precision function cbrt(x)
1861 cbrt = sign(abs(x)**(1/3d0), x)
1866 double precision function csmgt(x, y, p)
1868 double precision x, y
1880 double precision function trgsum(sinp, sinx, cosx, c, n)
1889 double precision sinx, cosx, c(n)
1891 double precision ar, y0, y1
1895 ar = 2 * (cosx - sinx) * (cosx + sinx)
1897 if (mod(n, 2) .eq. 1)
then
1908 y1 = ar * y0 - y1 + c(k)
1909 y0 = ar * y1 - y0 + c(k-1)
1913 trgsum = 2 * sinx * cosx * y0
1916 trgsum = cosx * (y0 - y1)
1922 integer function trnsit(lon1, lon2)
1924 double precision lon1, lon2
1926 double precision lon1x, lon2x, lon12, AngNm, AngDif
1929 lon12 = angdif(lon1x, lon2x)
1931 if (lon1x .lt. 0 .and. lon2x .ge. 0 .and. lon12 .gt. 0)
then
1933 else if (lon2x .lt. 0 .and. lon1x .ge. 0 .and. lon12 .lt. 0)
then
1940 subroutine accini(s)
1943 double precision s(2)
1951 subroutine accadd(s, y)
1956 double precision s(2)
1958 double precision z, u, sumx
1959 z = sumx(y, s(2), u)
1960 s(1) = sumx(z, s(1), s(2))
1961 if (s(1) .eq. 0)
then