1 /* -- translated by f2c (version 20100827).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
7 cc *.o -lf2c -lm
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10 http://www.netlib.org/f2c/libf2c.zip
11 */
12
13 #include "libtinyf2c.h"
14
15 /* ----------------------------------------------------------------------- */
16 /* INV MAP PROJECTION (NEWTON METHOD ON SPHERE) 2007-11-10 E. TOYODA */
17 /* ----------------------------------------------------------------------- */
18 /* Copyright (C) 2000-2007 GFD Dennou Club. All rights reserved. */
19 /* ----------------------------------------------------------------------- */
mpnwtn_(real * x,real * y,real * xlon,real * ylat,S_fp func,logical * xcyc)20 /* Subroutine */ int mpnwtn_(real *x, real *y, real *xlon, real *ylat, S_fp
21 func, logical *xcyc)
22 {
23 /* System generated locals */
24 real r__1, r__2;
25
26 /* Builtin functions */
27 double atan(doublereal);
28
29 /* Local variables */
30 static integer i__;
31 static real x0, y0, x1, y1, x2, y2, xb, yb, dx, dy, cx, cy, xl, yl, det,
32 rna, rdet, dist, ylat2, xlon2, accel;
33 extern /* Subroutine */ int glrget_(char *, real *, ftnlen);
34
35 accel = 1.f;
36 i__ = 0;
37 L100:
38 if (*xlon > 0.f) {
39 xlon2 = *xlon - 1e-4f;
40 } else {
41 xlon2 = *xlon + 1e-4f;
42 }
43 if (*ylat > 0.f) {
44 ylat2 = *ylat - 1e-4f;
45 } else {
46 ylat2 = *ylat + 1e-4f;
47 }
48 (*func)(xlon, ylat, &x0, &y0);
49 (*func)(&xlon2, ylat, &x1, &y1);
50 (*func)(xlon, &ylat2, &x2, &y2);
51 /* print "('foot:',4f13.6)", degree(xlon), degree(xlon2), */
52 /* + degree(ylat), degree(ylat2) */
53 xl = (x1 - x0) / (xlon2 - *xlon);
54 xb = (x2 - x0) / (ylat2 - *ylat);
55 yl = (y1 - y0) / (xlon2 - *xlon);
56 yb = (y2 - y0) / (ylat2 - *ylat);
57 dx = *x - x0;
58 dy = *y - y0;
59 det = xl * yb - xb * yl;
60 if (abs(det) < 1e-7f) {
61 /* print "('singularity')" */
62 goto L110;
63 }
64 rdet = accel / det;
65 /* print "('det=',f15.8)", det */
66 cx = rdet * (yb * dx - xb * dy);
67 cy = rdet * (-yl * dx + xl * dy);
68 /* print "('corr=',2f13.6)", degree(cx), degree(cy) */
69 *xlon += cx;
70 if (*xlon > atan(1.f) * 4.f + 1e-4f) {
71 /* print "(' xlon > 180')" */
72 if (*xcyc) {
73 *xlon -= atan(1.f) * 4.f * 2.f;
74 } else {
75 *xlon = atan(1.f) * 4.f - .01f;
76 }
77 accel = .5f;
78 } else if (*xlon < -(atan(1.f) * 4.f) - 1e-4f) {
79 /* print "(' xlon < -180')" */
80 if (*xcyc) {
81 *xlon += atan(1.f) * 4.f * 2.f;
82 } else {
83 *xlon = -(atan(1.f) * 4.f) + .01f;
84 }
85 accel = .5f;
86 }
87 *ylat += cy;
88 if (*ylat > atan(1.f) * 4.f * .5f) {
89 /* print "(' ylat > 90')" */
90 *ylat = atan(1.f) * 4.f * .5f - .01f;
91 accel = .5f;
92 } else if (*ylat < atan(1.f) * 4.f * -.5f) {
93 /* print "(' ylat < -90')" */
94 *ylat = atan(1.f) * 4.f * -.5f + .01f;
95 accel = .5f;
96 }
97 /* Computing MAX */
98 r__1 = abs(cx), r__2 = abs(cy);
99 dist = max(r__1,r__2);
100 if (dist < 3e-5f) {
101 goto L110;
102 }
103 ++i__;
104 if (i__ > 40) {
105 goto L911;
106 }
107 goto L100;
108 L110:
109 return 0;
110 L911:
111 glrget_("RUNDEF", &rna, (ftnlen)6);
112 *xlon = rna;
113 *ylat = rna;
114 return 0;
115 } /* mpnwtn_ */
116
117