1 /*---------------------------------------------------------------------------
2 *
3 *
4 * File: cm4_functions.c
5 *
6 * Functions required to compute CM4 magnetic components
7 *
8 * Authors: J. Luis translated from original Fortran code
9 * P. Wessel further massaged it into this form
10 *
11 * Version: 1.0
12 * Revised: 1-MAY-2009
13 *
14 *
15 * NOTES: The original Fortran code written by Terry Sabaka from
16 * - Planetary Geodynamics Lab at Goddard Space Flight Center -
17 * can be found at
18 * denali.gsfc.nasa.gov/cm/cm4field.f
19 * This C version is a bit more limited (it doesn't allow computing the
20 * source coefficients - the GMDL array in original) and was striped of
21 * the long help/comments section. Many of those comments make no sense
22 * here since we changed the subroutine interface. With regard to this point,
23 * a substantial difference is that all is need is one single call to the
24 * function (with location arrays transmitted in input), and all selected
25 * field sources contribution add up to the final result.
26 *
27 *-------------------------------------------------------------------------*/
28
29 #include "mgd77.h"
30 #include "cm4_functions.h"
31
32 #define I_DIM(x, y) (((x) > (y)) ? (x) - (y) : 0)
33
34 static void ymdtomjd(int yearad, int month, int dayofmonth, int *mjd, int *dayofyear);
35 static void ydtomjdx(int yearad, int dayofyear, int * mjd, int *month, int *dayofmonth, int *daysinmonth);
36 static double intdst(int mjdl, int mjdh, int mjdy, int msec, double *dstx, int *cerr);
37 static double intf107(int iyrl, int imol, int iyrh, int imoh, int iyr, int imon, int idom, int *idim,
38 int msec, double *f107x, int *cerr);
39 static double getmut2(double thenmp, double phinmp, int iyear, int iday, int msec);
40 static void sun2(int iyr, int iday, double secs, double *gst, double *slong, double *srasn, double *sdec);
41 static void rmerge_(double *rmrg, double *rmlt);
42 static void tsearad(int full, int ks, int kr, int ns, int ng, double f, double *t, double *e, double *g);
43 static void tseardr(int full, int ks, int kr, int ns, int ng, double f, double *t, double *e, double *g);
44 static void mseason(int ks, int ns, int ng, double d, double *t, double *e, double *g);
45 static void iseason(int ks, int ns, int ng, double f, double *t, double *e, double *g);
46 static void mpotent(int nmax, int mmax, int nd, int nz, double cphi, double sphi, double *d, double *z);
47 static void jtbelow(int pmin, int pmax, int nmax, int mmax, double r, double rref, int nz, double *z);
48 static void jtabove(int pmin, int pmax, int nmax, int mmax, double r, double rref, int nz, double *z);
49 static void jtbcont(int pmin, int pmax, int nmax, int mmax, double rold, double rnew, int nz, double *z);
50 static void mstream(int nmax, int mmax, int nd, int nz, double cphi, double sphi, double faor, double *d, double *z);
51 static void jpoloid(int pmin, int pmax, int nmax, int mmax, double r, double rm, int nd, int nz, double *t, double *d, double *z);
52 static void blsgen(int nc, int nd, int ni, double *b, double *c, double *dldc);
53 static void getgmf(int nder, int ns, double *ep, double *tm, double *b, double *c, double *g, int *h, int *o, double *p);
54 static void dbspln_(int *l, double *t, int *n, int * d__, int *k, double *x, double *b, double *w);
55 static void getgxf(int pmin, int pmax, int nmax, int mmax, int *ng, double *e, double *g, double *t);
56 static void bfield(int rgen, int nmxi, int nmxe, int nmni, int nmne, int mmxi, int mmxe, int mmni,
57 int mmne, int grad, int ctyp, int dtyp, int ityp, int etyp, double ep, double re,
58 double rp, double rm, double tm, double clat, double elon, double h, double dst, double dstt,
59 double *rse, int *nc, int *na, double *ro, double *theta, int *atyp, int *dsti, int *bori, int *bkni,
60 double *bkpi, int *tdgi, int *dste, int *bore, int *bkne, double *bkpe, int *tdge, double *a,
61 double *b, double *c, double *p, double *r, double *t, int *u, double *w, double *dsdc,
62 double *dldc, double *dlda, int *cerr);
63 static void prebf_(int *rgen, int *ityp, int *etyp, int *dtyp, int *grad, int *nmni, int *nmxi, int *
64 nmne, int *nmxe, int *mmni, int *mmxi, int *mmne, int *mmxe, int *nmax, int *mmin, int *mmax, int *
65 ns, int *nsi, int *nse, int *nc, int *nci, int *nce, int *na, int *np, int *ii, int *ie, int *
66 atyp, int *dsti, int *bori, int *bkni, int *tdgi, int *dste, int *bore, int *bkne, int *tdge, int *u, int *cerr);
67 static void fdlds_(int *rgen, int *grad, int *ctyp, double *clat, double *phi, double *h, double *re,
68 double *rp, double *rm, double *ro, int *nsi, int *nc, int *nci, int *np, int *ii, int *ie, int *
69 nmni, int *nmxi, int *nmne, int *nmxe, int *nmax, int *mmni, int *mmxi, int *mmne, int *mmxe, int *
70 mmin, int *mmax, double *theta, double *p, double *r, double *t, int *u, double *w, double *dldc, int *cerr);
71 static void geocen(int ctyp, double re, double rp, double rm, double h, double clat, double *r, double *theta, double *sinthe, double *costhe);
72 static void schmit_(int *grad, int *rgen, int *nmax, int *mmin, int *mmax, double *sinthe, double *costhe, double *p, double *r);
73 static void srecur_(int *grad, int *nmax, int *mmin, int *mmax, int *ksm2, int *ktm2, int *npall, int *
74 nad1, int *nad2, int *nad3, int *nad4, int *nad5, int *nad6, int *nad7, int *nad8, double *r);
75 static void trigmp(int mmax, double phi, double *t);
76 static void tdc(int grad, int nc, double clat, double theta, double *dldc, double *r);
77 static void fdsdc_(int *rgen, int *ityp, int *etyp, int *nsi, int *nse, int *nc, int *nci, double *ta,
78 double *tb, double *dst, double *dstt, int *dsti, int *bori, int *bkni, double *bkpi, int *tdgi,
79 int *dste, int *bore, int *bkne, double *bkpe, int *tdge, int *u, double *w, double *dsdc, int *cerr);
80 static void taylor(int nc, int ns, double ta, double tb, int *tdeg, int *u, double *dsdt, double *dsdc);
81 static void bsplyn(int nc, int ns, double *ta, double *tb, int *bord, int *bkno, double *bkpo, int *u, double *dtdb, double *dsdc, int *cerr);
82 static void sbspln_(double *ta, double *tb, int *n, int *k, double *bkpo, double *dtdb, double *dsdc, int *cerr);
83 static void tbspln_(double *t, int *n, int *k, double *bkpo, double *dtdb, int *cerr);
84 static void dstorm(int nc, int ns, double *dst, double *dstt, int *dstm, int *u, double *dsdc);
85 static void fdldc(int grad, int nc, double *dsdc, double *dldc);
86 static void blgen(int grad, int nc, double *b, double *c, double *dldc);
87 static void bngen_(double *b);
88 static void tec(int grad, int k, int nc, double *theta, double *phi, double *b, double *dldc, double *r);
89 static void tse(int grad, int k, int nc, double *rse, double *b, double *dldc, double *r);
90 static void tms(int grad, int k, int nc, int na, int ia, double *a, double *b, double *dldc, double *dlda, double *r);
91 static void fdldeu_(int *k, int *na, int *ia, double *seulx, double *ceulx, double *seuly, double *ceuly,
92 double *seulz, double *ceulz, double *r, double *b, double *dlda);
93 static void tnm_(int *grad, int *k, int *nc, int *na, int *ia, double *a, double *b, double *dldc, double *dlda, double *r);
94 static void fdldno_(int *k, int *na, int *ia, double *schix, double *cchix, double *schiy, double *cchiy,
95 double *schiz, double *cchiz, double *r, double *b, double *dlda);
96 static void fdldsl_(int *k, int *na, int *ia, double *b, double *dlda);
97 static void tvn_(int *grad, int *k, int *nc, int *na, int *ia, double *a, double *b, double *dldc, double *dlda, double *r);
98 static void tbi_(int *k, int *na, int *ia, double *a, double *b, double *dlda);
99 static void fdldbi_(int *k, int *na, int *ia, double *dlda);
100 static void ltrans(int n, int m, double *q, double *r, double *s);
101 static void ltranv(int rfac, int n, int m, double *r, double *v);
102 static int nshx(int nmax, int nmin, int mmax, int mmin);
103 static int nlpx(int nmax, int mmax, int mmin);
104 static int i8ssum(int abeg, int alen, int *a);
105 static void i8vset(int abeg, int alen, int s, int *a);
106 static void i8vadd(int abeg, int bbeg, int cbeg, int vlen, int *a, int *b, int *c);
107 static void i8vadds(int abeg, int bbeg, int vlen, int s, int *a, int *b);
108 static void i8vcum(int abas, int abeg, int alen, int *a);
109 static void i8vdel(int abas, int abeg, int alen, int *a);
110 static void r8vset(int abeg, int alen, double s, double *a);
111 static double r8sdot(int abeg, int bbeg, int vlen, double *a, double *b);
112 static double r8ssum_(int *abeg, int *alen, double *a);
113 static void r8slt(int abeg, int alen, double s, double *a, int *j);
114 static void r8vsub(int abeg, int bbeg, int cbeg, int vlen, double *a, double *b, double *c);
115 static void r8vmul(int abeg, int bbeg, int cbeg, int vlen, double *a, double *b, double *c);
116 static void r8vscale(int abeg, int alen, double s, double *a);
117 static void r8vscats(int qbeg, int qlen, double s, int *q, double *a);
118 static void r8vlinkt(int abeg, int bbeg, int vlen, double s, double *a, double *b);
119 static void r8vlinkq(int abeg, int bbeg, int cbeg, int vlen, double s, double *a, double *b, double *c);
120 static void r8vgathp(int abeg, int ainc, int bbeg, int blen, double *a, double *b);
121 static double d_mod(double x, double y);
122 static double pow_di(double ap, int bp);
123 static int i_dnnt(double x);
124 static void clear_mem (double *mut, double *gpsq, double *gssq, double *gpmg, double *gsmg, double *hysq, double *epsq, double *essq,
125 double *ecto, double *hyto, double *hq, double *ht, double *bkpo, double *ws, double *gamf, double *epmg,
126 double *esmg, double *hymg, double *f107x, double *pleg, double *rcur, double *gcto_or, double *gcto_mg);
127
MGD77_cm4field(struct GMT_CTRL * GMT,struct MGD77_CM4 * Ctrl,double * p_lon,double * p_lat,double * p_alt,double * p_date)128 int MGD77_cm4field (struct GMT_CTRL *GMT, struct MGD77_CM4 *Ctrl, double *p_lon, double *p_lat, double *p_alt, double *p_date) {
129
130 int c__1356 = 1356, c__13680 = 13680;
131 int i, j, k, l, n, p, nu, mz, nz, mu, js, jy, nt, mt, iyr = 0, jyr = 0, jf107, cerr = 0;
132 int lum1, lum2, lum3, lum4, lum5, lum6, lum7, nsm1, nsm2, lcmf, idim[12], omdl;
133 int lsmf, lpos, lcmg, lsmg, lcsq, lssq, lcto, lsto, lrto, idoy, n_Dst_rows, i_unused = 0;
134 int *msec = NULL, *mjdy = NULL;
135 int imon = 0, idom = 0, jaft, jmon = 0, jdom, jmjd = 0, jdoy, mjdl = 0, mjdh = 0, iyrl = 0, imol = 0, iyrh = 0, imoh = 0;
136 int nout = 0, nygo = 0, nmax, nmin, nobo, nopo, nomn, nomx, noff, noga, nohq, nimf, nyto, nsto, ntay, mmdl;
137 int us[4355], bord[4355], bkno[4355], pbto, peto, csys, jdst[24];
138 double *mut = NULL, *dstx = NULL, dstt = 0., x, y, z, h, t, dumb, bmdl[21], jmdl[12], date, dst, mut_now, alt;
139 double re, xd, yd, rm, xg, ro = 0, rp, yg, zg, zd;
140 double bc[29], wb[58], trig[132], ru, rt, rse[9], doy, fyr, cego, epch;
141 double rlgm[15], rrgt[9], tsmg[6], tssq[6], tsto[6], tdmg[12], tdsq[10], tdto[10];
142 double rtay_dw, rtay_or, sinp, fsrf, rtay, frto, frho, thetas = 0, rtay_dk;
143 double cnmp, enmp, omgs, omgd, hion, cpol, epol, ctmp, stmp, cemp, semp, rion, fdoy, clat, elon;
144 double sthe, cthe, psiz, cpsi, spsi, ctgo, stgo, sego, cdip = 0, edip = 0, ctmo, stmo, cemo, semo, taus = 0, taud = 0, cosp;
145 double *hq = NULL, *ht = NULL, *pleg = NULL, *rcur = NULL; /* was hq[53040], ht[17680], pleg[4422], rcur[9104] */
146 double *bkpo = NULL, *ws = NULL, *gamf = NULL, *epmg = NULL, *esmg = NULL; /* was bkpo[12415], ws[4355], gamf[8840], epmg[1356], esmg[1356] */
147 double *f107x = NULL; /* was [100][12] */
148 double *hymg = NULL; /* was [1356][6] */
149 double *gcto_or = NULL; /* was [13680][5][2] */
150 double *gcto_mg = NULL; /* was [2736][3][2][2] */
151 double *gpsq = NULL; /* was [13680][5][2] */
152 double *gssq = NULL; /* was [13680][5] */
153 double *gpmg = NULL; /* was [1356][5][2] */
154 double *gsmg = NULL; /* was [1356][5][2] */
155 double *hysq = NULL; /* was [1356][6] */
156 double *epsq = NULL; /* was [13680] */
157 double *essq = NULL; /* was [13680] */
158 double *ecto = NULL; /* was [16416] */
159 double *hyto = NULL; /* was [49248] */
160 char line[GMT_BUFSIZ] = {""}, *c_unused = NULL;
161 void *ptr = NULL;
162
163 FILE *fp;
164 gmt_M_unused(GMT);
165
166 /* ===== FORTRAN SOUVENIRS ==============================================
167 =======================================================================
168 Main Field potential expansion parameters
169 PARAMETER (NXMF=65,NYMF=8840,NXOR=4,NXKN=19,NXPO=12415)
170 =======================================================================
171 Magnetospheric and coupled Induction potential expansion parameters
172 PARAMETER (PBMG=0,PEMG=5,PSMG=2,NXMG=11,MXMG=6,IXMG=226)
173 =======================================================================
174 Sq and coupled Induction potential expansion parameters
175 PARAMETER (PBSQ=0,PESQ=4,PSSQ=2,NXSQ=60,MXSQ=12,IXSQ=2736)
176 =======================================================================
177 Toroidal scalar or stream function expansion parameters
178 PARAMETER (PBTO_MG=0,PETO_MG=0,PBTO_OR=0,PETO_OR=4)
179 PARAMETER (PSTO=2,NXTO=60,MXTO=12,IXTO=2736)
180 PARAMETER (NTAY_MG=1,NTAY_OR=1)
181 ======================================================================= */
182
183 bkpo = calloc(12415U, sizeof(double));
184 gamf = calloc(8840U, sizeof(double));
185 f107x = calloc(1200U, sizeof(double));
186
187 if ((fp = fopen(Ctrl->CM4_M.path, "r")) == NULL) {
188 GMT_Report (GMT->parent, GMT_MSG_ERROR, "CM4: Could not open file %s\n", Ctrl->CM4_M.path);
189 gmt_M_str_free (bkpo); gmt_M_str_free (gamf); gmt_M_str_free (f107x);
190 return 1;
191 }
192
193 /* coverity[tainted_data_argument] */ /* Try to shut up Coverity about Tainted variables */
194 c_unused = fgets(line, GMT_BUFSIZ, fp);
195 while (line[0] == '#') /* Skip comments */
196 c_unused = fgets(line, GMT_BUFSIZ, fp);
197
198 if ((n = sscanf (line, "%d %d %d", &lsmf, &lpos, &lcmf)) != 3) {
199 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Failed to parse line in MGD77_cm4field\n");
200 gmt_M_str_free (bkpo); gmt_M_str_free (gamf); gmt_M_str_free (f107x);
201 fclose (fp);
202 return 1;
203 }
204
205 if (lsmf > 4355 || lpos > 12415 || lcmf > 8840) { /* Mainly to shut up CID 39232 (TAINTED variables), but it maybe right */
206 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Suspicious values in first line of umdl.CM4 file\n");
207 gmt_M_str_free (bkpo); gmt_M_str_free (gamf); gmt_M_str_free (f107x);
208 fclose (fp);
209 return 1;
210 }
211
212 c_unused = fgets(line, GMT_BUFSIZ, fp);
213 (void)(c_unused++); /* silence -Wunused-but-set-variable and PVS warning of double assignment */
214 sscanf (line, "%d", &lum1);
215 c_unused = fgets(line, GMT_BUFSIZ, fp);
216 (void)c_unused; /* silence -Wunused-but-set-variable */
217 sscanf (line, "%lf %lf %lf %lf", &epch, &re, &rp, &rm);
218 for (j = 0; j < lsmf; ++j)
219 i_unused = fscanf (fp, "%d", &bord[j]);
220 for (j = 0; j < lsmf; ++j)
221 i_unused = fscanf (fp, "%d", &bkno[j]);
222 for (j = 0; j < lpos; ++j)
223 i_unused = fscanf (fp, "%lf", &bkpo[j]);
224 for (j = 0; j < lcmf; ++j)
225 i_unused = fscanf (fp, "%lf", &gamf[j]);
226
227 i_unused = fscanf (fp, "%d %d", &lcmg, &lsmg);
228 i_unused = fscanf (fp, "%d %d %d %d %d %d", &lum1, &lum2, &lum3, &lum4, &lum5, &lum6);
229 i_unused = fscanf (fp, "%lf %lf %lf %lf %lf %lf %lf", &cnmp, &enmp, &omgs, &omgd, &re, &rp, &rm);
230 /* coverity[overflow_sink] */ /* For Coverity analysis. Do not remove this comment */
231 gpmg = calloc ((2 * (size_t)lsmg * (size_t)lcmg), sizeof(double));
232 for (k = 0; k < 2; ++k)
233 for (j = 0; j < lsmg; ++j) {
234 n = (j + k * 5) * 1356;
235 for (i = 0; i < lcmg; ++i)
236 i_unused = fscanf (fp, "%lf", &gpmg[i + n]);
237 }
238
239 /* coverity[overflow_sink] */ /* For Coverity analysis. Do not remove this comment */
240 gsmg = calloc((2 * (size_t)lsmg * (size_t)lcmg), sizeof(double));
241 for (k = 0; k < 2; ++k)
242 for (j = 0; j < lsmg; ++j) {
243 n = (j + k * 5) * 1356;
244 for (i = 0; i < lcmg; ++i)
245 i_unused = fscanf (fp, "%lf", &gsmg[i + n]);
246 }
247
248 i_unused = fscanf (fp, "%d %d", &lcsq, &lssq);
249 i_unused = fscanf (fp, "%d %d %d %d %d %d", &lum1, &lum2, &lum3, &lum4, &lum5, &lum6);
250 i_unused = fscanf (fp, "%lf %lf %lf %lf %lf %lf %lf %lf", &cnmp, &enmp, &omgs, &omgd, &re, &rp, &rm, &hion);
251 gpsq = calloc((2 * (size_t)lssq * (size_t)lcsq), sizeof(double));
252 for (k = 0; k < 2; ++k)
253 for (j = 0; j < lssq; ++j) {
254 n = (j + k * 5) * 13680;
255 for (i = 0; i < lcsq; ++i)
256 i_unused = fscanf (fp, "%lf", &gpsq[i + n]);
257 }
258
259 gssq = calloc(((size_t)lssq * (size_t)lcsq), sizeof(double));
260 for (j = 0; j < lssq; ++j) {
261 n = j * 13680;
262 for (i = 0; i < lcsq; ++i)
263 i_unused = fscanf (fp, "%lf", &gssq[i + n]);
264 }
265
266 i_unused = fscanf (fp, "%d %d %d", &lcto, &lsto, &lrto);
267 i_unused = fscanf (fp, "%d %d %d %d %d %d %d", &lum1, &lum2, &lum3, &lum4, &lum5, &lum6, &lum7);
268 i_unused = fscanf (fp, "%lf %lf %lf %lf %lf %lf %lf %lf %lf", &cnmp, &enmp, &omgs, &omgd, &re, &rp, &rm, &rtay_dw, &rtay_dk);
269 if (Ctrl->CM4_DATA.pred[3]) { /* In other cases the next coefficients are not used, so no waste time/memory with them */
270 gcto_mg = calloc((2 * (size_t)lrto * (size_t)lsto * (size_t)lcto), sizeof(double));
271 for (l = 0; l < 2; ++l)
272 for (k = 0; k < lrto; ++k)
273 for (j = 0; j < lsto; ++j) {
274 n = (j + (k + (l << 1)) * 3) * 2736;
275 for (i = 0; i < lcto; ++i)
276 i_unused = fscanf (fp, "%lf", &gcto_mg[i + n]);
277 }
278 }
279 else /* Jump the unused coeffs */
280 for (l = 0; l < 2 * lrto * lsto * lcto; ++l)
281 i_unused = fscanf (fp, "%lf", &dumb);
282
283 i_unused = fscanf (fp, "%d %d %d", &lcto, &lsto, &lrto);
284 i_unused = fscanf (fp, "%d %d %d %d %d %d %d", &lum1, &lum2, &lum3, &lum4, &lum5, &lum6, &lum7);
285 i_unused = fscanf (fp, "%lf %lf %lf %lf %lf %lf %lf %lf", &cnmp, &enmp, &omgs, &omgd, &re, &rp, &rm, &rtay_or);
286 if (Ctrl->CM4_DATA.pred[3] && !Ctrl->CM4_DATA.pred[4]) { /* In other cases the next coefficients are not used, so no waste time/memory with them */
287 gcto_or = calloc(((size_t)lrto * (size_t)lsto * (size_t)lcto), sizeof(double));
288 for (k = 0; k < lrto; ++k)
289 for (j = 0; j < lsto; ++j) {
290 n = (j + k * 5) * 13680;
291 for (i = 0; i < lcto; ++i)
292 i_unused = fscanf (fp, "%lf", &gcto_or[i + n]);
293 }
294 }
295 (void)i_unused; /* silence -Wunused-but-set-variable */
296
297 fclose(fp);
298 cpol = cnmp * D2R;
299 epol = enmp * D2R;
300 sincos(cpol, &stmp, &ctmp);
301 sincos(epol, &semp, &cemp);
302 rion = rm + hion;
303
304 mut = calloc((size_t)(Ctrl->CM4_DATA.n_times), sizeof(double));
305 msec = calloc((size_t)(Ctrl->CM4_DATA.n_times), sizeof(int));
306 mjdy = calloc((size_t)(Ctrl->CM4_DATA.n_times), sizeof(int));
307 for (n = 0; n < Ctrl->CM4_DATA.n_times; ++n) { /* If time is not constant compute the mut array */
308 iyr = (int)(p_date[n]);
309 fyr = p_date[n] - (double) iyr;
310 doy = fyr * (double) (366 - MIN(1, iyr % 4));
311 idoy = (int) doy;
312 fdoy = doy - (double) idoy;
313 ++idoy;
314 msec[n] = i_dnnt(fdoy * 8.64e7);
315 ydtomjdx(iyr, idoy, &mjdy[n], &imon, &idom, idim);
316 mut[n] = getmut2(cnmp, enmp, iyr, idoy, msec[n]);
317 }
318
319 csys = 1;
320 if (Ctrl->CM4_G.geodetic) csys = 0;
321 if (Ctrl->CM4_D.index) {
322 if (Ctrl->CM4_D.load) {
323 int k;
324 if ((fp = fopen(Ctrl->CM4_D.path, "r")) == NULL) {
325 fprintf (stderr, "CM4: Could not open file %s\n", Ctrl->CM4_D.path);
326 clear_mem (mut, gpsq, gssq, gpmg, gsmg, hysq, epsq, essq, ecto, hyto, hq, ht, bkpo, ws, gamf, epmg,
327 esmg, hymg, f107x, pleg, rcur, gcto_or, gcto_mg);
328 gmt_M_str_free (msec); gmt_M_str_free (mjdy); free (mjdy);
329 return 1;
330 }
331 jaft = 0;
332 n = 0;
333 n_Dst_rows = 21884 * 24; /* Current (13-05-2009) number of lines in Dst_all.wdc file */
334 dstx = calloc(((size_t)n_Dst_rows), sizeof(double));
335 /* One improvement would be to compute year_min/year_max and retain only the needed data in dstx */
336
337 while (fgets (line, GMT_BUFSIZ, fp)) {
338 if (line[0] == '#') continue; /* Comment line only */
339 sscanf (&line[3], "%2d %2d", &jyr, &jmon);
340 sscanf (&line[8], "%2d", &jdom);
341 for (i = 0; i < 24; ++i)
342 sscanf (&line[20+i*4],"%4d", &jdst[i]);
343 jyr += 1900;
344 if (jyr < 1957) jyr += 100;
345 ymdtomjd(jyr, jmon, jdom, &jmjd, &jdoy);
346 if (jaft == 0) {
347 jaft = 1;
348 mjdl = jmjd;
349 }
350 k = (jmjd - mjdl) * 24;
351 if (k >= n_Dst_rows) {
352 n_Dst_rows += (100 * 24);
353 ptr = (void *)dstx;
354 if ((dstx = realloc(dstx, (size_t)(n_Dst_rows * 24) * sizeof(double))) == NULL) {
355 gmt_M_str_free (ptr);
356 return 1;
357 }
358 }
359 for (j = 0; j < 24; ++j)
360 dstx[k + j] = (double)jdst[j];
361 n++;
362 }
363 fclose(fp);
364 mjdh = jmjd;
365 }
366 if (Ctrl->CM4_DATA.n_times > 1) { /* Need to re-allocate memory for all n_times in dst array */
367 ptr = (void *)Ctrl->CM4_D.dst;
368 if ((Ctrl->CM4_D.dst = realloc(Ctrl->CM4_D.dst, (size_t)(Ctrl->CM4_DATA.n_times) * sizeof(double))) == NULL) {
369 gmt_M_str_free (ptr);
370 return 1;
371 }
372 }
373
374 /* Get only one dst first so that we can test (and abort if needed) if date is out of bounds */
375 if (dstx) Ctrl->CM4_D.dst[0] = intdst(mjdl, mjdh, mjdy[0], msec[0], dstx, &cerr);
376 if (cerr > 49) {
377 gmt_M_str_free (dstx);
378 if (Ctrl->CM4_DATA.n_times > 1) gmt_M_str_free (Ctrl->CM4_D.dst);
379 clear_mem (mut, gpsq, gssq, gpmg, gsmg, hysq, epsq, essq, ecto, hyto, hq, ht, bkpo, ws, gamf, epmg,
380 esmg, hymg, f107x, pleg, rcur, gcto_or, gcto_mg);
381 gmt_M_str_free (msec); gmt_M_str_free (mjdy); free (mjdy);
382 return 1;
383 }
384
385
386 if (dstx) {
387 for (n = 1; n < Ctrl->CM4_DATA.n_times; n++)
388 Ctrl->CM4_D.dst[n] = intdst(mjdl, mjdh, mjdy[n], msec[n], dstx, &cerr);
389
390 gmt_M_str_free (dstx);
391 }
392 if (cerr > 49) {
393 clear_mem (mut, gpsq, gssq, gpmg, gsmg, hysq, epsq, essq, ecto, hyto, hq, ht, bkpo, ws, gamf, epmg,
394 esmg, hymg, f107x, pleg, rcur, gcto_or, gcto_mg);
395 gmt_M_str_free (msec); gmt_M_str_free (mjdy); free (mjdy);
396 return 1;
397 }
398 }
399 if (Ctrl->CM4_I.index) {
400 if (Ctrl->CM4_I.load) {
401 if ((fp = fopen(Ctrl->CM4_I.path, "r")) == NULL) {
402 fprintf (stderr, "CM4: Could not open file %s\n", Ctrl->CM4_I.path);
403 clear_mem (mut, gpsq, gssq, gpmg, gsmg, hysq, epsq, essq, ecto, hyto, hq, ht, bkpo, ws, gamf, epmg,
404 esmg, hymg, f107x, pleg, rcur, gcto_or, gcto_mg);
405 gmt_M_str_free (msec); gmt_M_str_free (mjdy); free (mjdy);
406 return 1;
407 }
408 jaft = 0;
409 while (fgets (line, GMT_BUFSIZ, fp)) {
410 if (line[0] == '#') continue; /* Comment line only */
411 if (line[9] != '-') {
412 sscanf (line, "%d %d %d", &jyr, &jmon, &jf107);
413 if (jaft == 0) {
414 jaft = 1;
415 iyrl = jyr;
416 imol = jmon;
417 }
418 f107x[(jyr - iyrl) * 12 + jmon-1] = (double)jf107 / 10.;
419 }
420 }
421 fclose(fp);
422 iyrh = jyr;
423 imoh = jmon;
424 }
425 /* MUST INVESTIGATE IF IT WORTH HAVING AN ARRAY OF f107 LIKE IN THE DST CASE */
426 Ctrl->CM4_I.F107 = intf107(iyrl, imol, iyrh, imoh, iyr, imon, idom, idim, msec[0], f107x, &cerr);
427 if (cerr > 49) {
428 gmt_M_str_free (msec); gmt_M_str_free (mjdy); free(mjdy);
429 clear_mem (mut, gpsq, gssq, gpmg, gsmg, hysq, epsq, essq, ecto, hyto, hq, ht, bkpo, ws, gamf, epmg,
430 esmg, hymg, f107x, pleg, rcur, gcto_or, gcto_mg);
431 return 1;
432 }
433 }
434 gmt_M_str_free (msec); gmt_M_str_free (mjdy); free (mjdy);
435
436 /* On Windows, either this or declare them as "static", otherwise ... BOOM */
437 hysq = calloc(82080U, sizeof(double));
438 epsq = calloc(13680U, sizeof(double));
439 essq = calloc(13680U, sizeof(double));
440 ecto = calloc(16416U, sizeof(double));
441 hyto = calloc(49248U, sizeof(double));
442 hq = calloc(53040U, sizeof(double));
443 ht = calloc(17680U, sizeof(double));
444 ws = calloc(4355U, sizeof(double));
445 epmg = calloc(1356U, sizeof(double));
446 esmg = calloc(1356U, sizeof(double));
447 hymg = calloc(8136U, sizeof(double));
448 pleg = calloc(4422U, sizeof(double));
449 rcur = calloc(9104U, sizeof(double));
450
451 /* LOOP over number of input points (many computations below are useless repeated - room for improvement */
452 for (n = 0; n < Ctrl->CM4_DATA.n_pts; ++n) {
453 memset(bmdl, 0, 21 * sizeof(double));
454 if (Ctrl->CM4_L.curr)
455 memset(jmdl, 0, 12 * sizeof(double));
456 clat = (90 - p_lat[n]) * D2R;
457 elon = p_lon[n] * D2R;
458
459 /* See if we are using a constant time or an array */
460 if (Ctrl->CM4_DATA.n_times > 1) {
461 date = p_date[n];
462 dst = Ctrl->CM4_D.dst[n];
463 mut_now = mut[n];
464 }
465 else {
466 date = p_date[0];
467 dst = Ctrl->CM4_D.dst[0];
468 mut_now = mut[0];
469 }
470
471 /* See if we are using a constant altitude or an array */
472 alt = (Ctrl->CM4_DATA.n_altitudes > 1) ? p_alt[n] : p_alt[0];
473
474 if (Ctrl->CM4_DATA.coef) {
475 nout = 1; nygo = 0;
476 if (Ctrl->CM4_DATA.pred[1]) nygo = MAX(nygo,113);
477 if (Ctrl->CM4_DATA.pred[2]) nygo = MAX(nygo,1368);
478 if (Ctrl->CM4_DATA.pred[3]) nygo = MAX(nygo,1368);
479 }
480 if (Ctrl->CM4_DATA.pred[0]) {
481 nmax = MAX(Ctrl->CM4_S.nhmf[0], Ctrl->CM4_S.nhmf[1]);
482 nmin = MIN(Ctrl->CM4_S.nlmf[0], Ctrl->CM4_S.nlmf[1]);
483 nobo = nshx(nmin - 1, 1, nmin - 1, 0);
484 nopo = i8ssum(1, nobo, bkno) + (nobo << 1);
485 bfield(1, nmax, 0, nmin, 1, nmax, 0, 0, 0, 0, csys, 3, 2, 0,
486 epch, re, rp, rm, date, clat, elon, alt, dst, dstt, rse, &nz,
487 &mz, &ro, &thetas, us, us, &bord[nobo], &bkno[nobo], &bkpo[nopo], us, us, us, us,
488 ws, us, gamf, bc, gamf, pleg, rcur, trig, us, ws, ht, hq, hq, &cerr);
489 if (cerr > 49) {
490 clear_mem (mut, gpsq, gssq, gpmg, gsmg, hysq, epsq, essq, ecto, hyto, hq, ht, bkpo, ws, gamf, epmg,
491 esmg, hymg, f107x, pleg, rcur, gcto_or, gcto_mg);
492 return 1;
493 }
494 nomn = nshx(Ctrl->CM4_S.nlmf[0] - 1, 1, Ctrl->CM4_S.nlmf[0] - 1, 0);
495 nomx = nshx(Ctrl->CM4_S.nhmf[0], 1, Ctrl->CM4_S.nhmf[0], 0);
496 noff = nomn - nobo;
497 nsm1 = I_DIM(nomx, nomn);
498 noga = i8ssum(1, nomn, bord) + i8ssum(1, nomn, bkno) + nomn;
499 nohq = i8ssum(nobo + 1, noff, bord) + i8ssum(nobo + 1, noff, bkno) + noff;
500 nimf = i8ssum(nobo + 1, nsm1, bord) + i8ssum(nobo + 1, nsm1, bkno) + nsm1;
501 blsgen(nimf, nz, 3, &bmdl[0], &gamf[noga], &hq[nohq]);
502 if (Ctrl->CM4_DATA.coef) {
503 nopo = i8ssum(1, nomn, bkno) + (nomn << 1);
504 getgmf(4, nsm1, &epch, &date, wb, &gamf[noga], &Ctrl->CM4_DATA.gmdl[nout-1],
505 &bkno[nomn], &bord[nomn], &bkpo[nopo]);
506 if (cerr > 49) {
507 clear_mem (mut, gpsq, gssq, gpmg, gsmg, hysq, epsq, essq, ecto, hyto, hq, ht, bkpo, ws, gamf, epmg,
508 esmg, hymg, f107x, pleg, rcur, gcto_or, gcto_mg);
509 return 1;
510 }
511 }
512 nomn = nshx(Ctrl->CM4_S.nlmf[1] - 1, 1, Ctrl->CM4_S.nlmf[1] - 1, 0);
513 nomx = nshx(Ctrl->CM4_S.nhmf[1], 1, Ctrl->CM4_S.nhmf[1], 0);
514 noff = nomn - nobo;
515 nsm2 = I_DIM(nomx, nomn);
516 noga = i8ssum(1, nomn, bord) + i8ssum(1, nomn, bkno) + nomn;
517 nohq = i8ssum(nobo + 1, noff, bord) + i8ssum(nobo + 1, noff, bkno) + noff;
518 nimf = i8ssum(nomn + 1, nsm2, bord) + i8ssum(nomn + 1, nsm2, bkno) + nsm2;
519 blsgen(nimf, nz, 3, &bmdl[3], &gamf[noga], &hq[nohq]);
520 if (Ctrl->CM4_DATA.coef) {
521 nygo = MAX(nygo,nsm1 * 5);
522 nygo = MAX(nygo,nsm2 * 5);
523 nout += nygo * MIN(1,nsm1);
524 nopo = i8ssum(1, nomn, bkno) + (nomn << 1);
525 getgmf(4, nsm2, &epch, &date, wb, &gamf[noga], &Ctrl->CM4_DATA.gmdl[nout-1],
526 &bkno[nomn], &bord[nomn], &bkpo[nopo]);
527 if (cerr > 49) {
528 clear_mem (mut, gpsq, gssq, gpmg, gsmg, hysq, epsq, essq, ecto, hyto, hq, ht, bkpo, ws, gamf, epmg,
529 esmg, hymg, f107x, pleg, rcur, gcto_or, gcto_mg);
530 return 1;
531 }
532 nout += nygo * MIN(1,nsm2);
533 }
534 }
535 if (Ctrl->CM4_DATA.pred[1] || Ctrl->CM4_DATA.pred[2] || Ctrl->CM4_DATA.pred[3]) {
536 if (!Ctrl->CM4_DATA.pred[0])
537 geocen(csys, re, rp, rm, alt, clat, &ro, &thetas, &sthe, &cthe);
538
539 psiz = thetas - clat;
540 sincos(psiz, &spsi, &cpsi);
541 rlgm[0] = -cpsi; rlgm[3] = 0.; rlgm[6] = spsi;
542 rlgm[1] = 0.; rlgm[4] = 1.; rlgm[7] = 0.;
543 rlgm[2] = -spsi; rlgm[5] = 0.; rlgm[8] = -cpsi;
544 sincos(thetas, &stgo, &ctgo);
545 sincos(elon, &sego, &cego);
546 rrgt[0] = ctgo * cego; rrgt[3] = -sego; rrgt[6] = stgo * cego;
547 rrgt[1] = ctgo * sego; rrgt[4] = cego; rrgt[7] = stgo * sego;
548 rrgt[2] = -stgo; rrgt[5] = 0.; rrgt[8] = ctgo;
549 rmerge_(rlgm, rrgt);
550 rrgt[0] = cemp * ctmp; rrgt[3] = semp * ctmp; rrgt[6] = -stmp;
551 rrgt[1] = -semp; rrgt[4] = cemp; rrgt[7] = 0.;
552 rrgt[2] = cemp * stmp; rrgt[5] = semp * stmp; rrgt[8] = ctmp;
553 rmerge_(rlgm, rrgt);
554 xg = stgo * cego;
555 yg = stgo * sego;
556 zg = ctgo;
557 xd = rrgt[0] * xg + rrgt[3] * yg + rrgt[6] * zg;
558 yd = rrgt[1] * xg + rrgt[4] * yg + rrgt[7] * zg;
559 zd = rrgt[2] * xg + rrgt[5] * yg + rrgt[8] * zg;
560 cdip = acos(zd);
561 edip = atan2(yd, xd);
562 ctmo = zd;
563 stmo = sin(cdip);
564 sincos(edip, &semo, &cemo);
565 rrgt[0] = -cemo * ctmo; rrgt[3] = -semo * ctmo; rrgt[6] = stmo;
566 rrgt[1] = -semo; rrgt[4] = cemo; rrgt[7] = 0.;
567 rrgt[2] = -cemo * stmo; rrgt[5] = -semo * stmo; rrgt[8] = -ctmo;
568 rmerge_(rlgm, rrgt);
569 taus = omgs * date;
570 taud = omgd * mut_now;
571 }
572 if (Ctrl->CM4_DATA.pred[1]) {
573 bfield(1, 11, 11, 1, 1, 6, 6, 0, 0, 0, 1, 3, 0, 0, epch, re, rp, rm,
574 date, cdip, edip, alt, dst, dstt, rse, &nu, &mu,
575 &ru, &thetas, us, us, us, us, ws, us, us, us, us, ws, us, gsmg, bc, gsmg, pleg, rcur,
576 trig, us, ws, ht, hq, hq, &cerr);
577 if (cerr > 49) {
578 clear_mem (mut, gpsq, gssq, gpmg, gsmg, hysq, epsq, essq, ecto, hyto, hq, ht, bkpo, ws, gamf, epmg,
579 esmg, hymg, f107x, pleg, rcur, gcto_or, gcto_mg);
580 return 1;
581 }
582 js = nu / 2;
583 jy = 1;
584 trigmp(2, taus, tsmg);
585 trigmp(5, taud, tdmg);
586 for (p = 0; p <= 5; ++p) {
587 cosp = tdmg[p];
588 sinp = tdmg[p + 6];
589 mpotent(11, 6, nu, c__1356, cosp, sinp, hq, & hymg[jy - 1]);
590 mpotent(11, 6, nu, c__1356, cosp, sinp, &hq[js], &hymg[jy + 4067]);
591 jy += 226;
592 }
593 bc[0] = bc[1] = bc[2] = 0.;
594 mseason(2, 5, c__1356, dst, tsmg, epmg, gpmg);
595 blsgen(c__1356, c__1356, 3, bc, epmg, &hymg[4068]);
596 ltrans(1, 1, bc, rlgm, &bmdl[6]);
597 bc[0] = bc[1] = bc[2] = 0.;
598 mseason(2, 5, c__1356, dst, tsmg, esmg, gsmg);
599 blsgen(c__1356, c__1356, 3, bc, esmg, hymg);
600 ltrans(1, 1, bc, rlgm, &bmdl[9]);
601 if (Ctrl->CM4_L.curr) {
602 bc[0] = bc[1] = bc[2] = 0.;
603 jtbelow(0, 5, 11, 6, ro, rm, c__1356, hymg);
604 blsgen(c__1356, c__1356, 3, bc, esmg, hymg);
605 ltrans(1, 1, bc, rlgm, &jmdl[0]);
606 }
607 if (Ctrl->CM4_DATA.coef) {
608 getgxf(0, 5, 11, 6, &js, epmg, &Ctrl->CM4_DATA.gmdl[nout-1], tdmg);
609 nout += nygo;
610 getgxf(0, 5, 11, 6, &js, esmg, &Ctrl->CM4_DATA.gmdl[nout-1], tdmg);
611 nout += nygo;
612 }
613 }
614 if (Ctrl->CM4_DATA.pred[2]) {
615 fsrf = Ctrl->CM4_I.F107 * .01485 + 1.;
616 if (ro < rion) {
617 bfield(1, 60, 60, 1, 1, 12, 12, 0, 0, 0, 1, 3, 0, 0, epch, re, rp, rm,
618 date, cdip, edip, alt, dst, dstt, rse, &nu,
619 &mu, &ru, &thetas, us, us, us, us, ws, us, us, us, us, ws, us, gssq,
620 bc, gssq, pleg, rcur, trig, us, ws, ht, hq, hq, &cerr);
621 if (cerr > 49) {
622 clear_mem (mut, gpsq, gssq, gpmg, gsmg, hysq, epsq, essq, ecto, hyto, hq, ht, bkpo, ws, gamf, epmg,
623 esmg, hymg, f107x, pleg, rcur, gcto_or, gcto_mg);
624 return 1;
625 }
626 js = nu / 2;
627 jy = 1;
628 trigmp(2, taus, tssq);
629 trigmp(4, taud, tdsq);
630 for (p = 0; p <= 4; ++p) {
631 cosp = tdsq[p];
632 sinp = tdsq[p + 5];
633 mpotent(60, 12, nu, c__13680, cosp, sinp, hq, &hysq[jy - 1]);
634 mpotent(60, 12, nu, c__13680, cosp, sinp, &hq[js], &hysq[jy + 41039]);
635 jy += 2736;
636 }
637 bc[0] = bc[1] = bc[2] = 0.;
638 iseason(2, 5, c__13680, fsrf, tssq, epsq, gpsq);
639 blsgen(c__13680, c__13680, 3, bc, epsq, &hysq[41040]);
640 ltrans(1, 1, bc, rlgm, &bmdl[12]);
641 bc[0] = bc[1] = bc[2] = 0.;
642 iseason(2, 5, c__13680, fsrf, tssq, essq, gssq);
643 blsgen(c__13680, c__13680, 3, bc, essq, hysq);
644 ltrans(1, 1, bc, rlgm, &bmdl[15]);
645 if (Ctrl->CM4_L.curr) {
646 bc[0] = bc[1] = bc[2] = 0.;
647 jtabove(0, 4, 60, 12, ro, rion, c__13680, &hysq[41040]);
648 blsgen(c__13680, c__13680, 3, bc, epsq, &hysq[41040]);
649 ltrans(1, 1, bc, rlgm, &jmdl[3]);
650 bc[0] = bc[1] = bc[2] = 0.;
651 jtbelow(0, 4, 60, 12, ro, rm, c__13680, hysq);
652 blsgen(c__13680, c__13680, 3, bc, essq, hysq);
653 ltrans(1, 1, bc, rlgm, &jmdl[6]);
654 }
655 if (Ctrl->CM4_DATA.coef) {
656 getgxf(0, 4, 60, 12, &js, epsq, &Ctrl->CM4_DATA.gmdl[nout-1], tdsq);
657 nout += nygo;
658 getgxf(0, 4, 60, 12, &js, essq, &Ctrl->CM4_DATA.gmdl[nout-1], tdsq);
659 nout += nygo;
660 }
661 }
662 else {
663 bfield(1, 60, 0, 1, 1, 12, 0, 0, 0, 0, 1, 3, 0, 0, epch, re, rp, rm,
664 date, cdip, edip, alt, dst, dstt,
665 rse, &nu, &mu, &ru, &thetas, us, us, us, us, ws, us, us, us, us, ws, us, gssq,
666 bc, gssq, pleg, rcur, trig, us, ws, ht, hq, hq, &cerr);
667 if (cerr > 49) {
668 clear_mem (mut, gpsq, gssq, gpmg, gsmg, hysq, epsq, essq, ecto, hyto, hq, ht, bkpo, ws, gamf, epmg,
669 esmg, hymg, f107x, pleg, rcur, gcto_or, gcto_mg);
670 return 1;
671 }
672 jy = 1;
673 trigmp(2, taus, tssq);
674 trigmp(4, taud, tdsq);
675 for (p = 0; p <= 4; ++p) {
676 cosp = tdsq[p];
677 sinp = tdsq[p + 5];
678 mpotent(60, 12, nu, c__13680, cosp, sinp, hq, &hysq[jy - 1]);
679 jy += 2736;
680 }
681 bc[0] = bc[1] = bc[2] = 0.;
682 iseason(2, 5, c__13680, fsrf, tssq, epsq, &gpsq[68400]);
683 blsgen(c__13680, c__13680, 3, bc, epsq, hysq);
684 ltrans(1, 1, bc, rlgm, &bmdl[12]);
685 bc[0] = bc[1] = bc[2] = 0.;
686 iseason(2, 5, c__13680, fsrf, tssq, essq, gssq);
687 blsgen(c__13680, c__13680, 3, bc, essq, hysq);
688 ltrans(1, 1, bc, rlgm, &bmdl[15]);
689 if (Ctrl->CM4_L.curr) {
690 bc[0] = bc[1] = bc[2] = 0.;
691 jtbelow(0, 4, 60, 12, ro, rion, c__13680, hysq);
692 blsgen(c__13680, c__13680, 3, bc, epsq, hysq);
693 ltrans(1, 1, bc, rlgm, &jmdl[3]);
694 bc[0] = bc[1] = bc[2] = 0.;
695 jtbcont(0, 4, 60, 12, rion, rm, c__13680, hysq);
696 blsgen(c__13680, c__13680, 3, bc, essq, hysq);
697 ltrans(1, 1, bc, rlgm, &jmdl[6]);
698 }
699 if (Ctrl->CM4_DATA.coef) {
700 getgxf(0, 4, 60, 12, &nu, epsq, &Ctrl->CM4_DATA.gmdl[nout-1], tdsq);
701 nout += nygo;
702 getgxf(0, 4, 60, 12, &nu, essq, &Ctrl->CM4_DATA.gmdl[nout-1], tdsq);
703 nout += nygo;
704 }
705 }
706 }
707 if (Ctrl->CM4_DATA.pred[3]) {
708 if (Ctrl->CM4_DATA.pred[4]) {
709 if (Ctrl->CM4_DATA.pred[5]) {
710 pbto = peto = 0;
711 nyto = 2736;
712 nsto = 3;
713 ntay = 1;
714 rtay = rtay_dw;
715 omdl = false;
716 mmdl = 1;
717 } else {
718 pbto = peto = 0;
719 nyto = 2736;
720 nsto = 3;
721 ntay = 1;
722 rtay = rtay_dk;
723 omdl = false;
724 mmdl = 2;
725 }
726 } else {
727 pbto = 0;
728 peto = 4;
729 nyto = 13680;
730 nsto = 5;
731 ntay = 1;
732 rtay = rtay_or;
733 omdl = true;
734 mmdl = 1;
735 }
736 bfield(1, 60, 0, 1, 1, 12, 0, 0, 0, 0, 1, 3, 0, 0, epch, re, rp, rm,
737 date, cdip, edip, 0., dst, dstt, rse, &nt, &mt,
738 &rt, &thetas, us, us, us, us, ws, us, us, us, us, ws, us, gcto_mg, bc, gcto_mg,
739 pleg, rcur, trig, us, ws, ht, hq, hq, &cerr);
740 if (cerr > 49) {
741 clear_mem (mut, gpsq, gssq, gpmg, gsmg, hysq, epsq, essq, ecto, hyto, hq, ht, bkpo, ws, gamf, epmg,
742 esmg, hymg, f107x, pleg, rcur, gcto_or, gcto_mg);
743 return 1;
744 }
745 frto = rm / ro;
746 frho = (ro - rtay) / rm;
747 jy = 1;
748 trigmp(2, taus, tsto);
749 trigmp(peto, taud, tdto);
750 for (p = pbto; p <= peto; ++p) {
751 cosp = tdto[p];
752 sinp = tdto[p + 1 + peto];
753 mstream(60, 12, nt, nyto, cosp, sinp, frto, hq, &hyto[jy - 1]);
754 jy += 2736;
755 }
756 bc[0] = bc[1] = bc[2] = 0.;
757 if (omdl)
758 tsearad(omdl, 2, ntay, nsto, nyto, frho, tsto, ecto, gcto_or);
759 else
760 tsearad(omdl, 2, ntay, nsto, nyto, frho, tsto, ecto, &gcto_mg[(((mmdl << 1) + 1) * 3 + 1) * 2736 - 27360]);
761
762 blsgen(nyto, nyto, 2, bc, ecto, hyto);
763 ltrans(1, 1, bc, rlgm, &bmdl[18]);
764 if (Ctrl->CM4_DATA.coef)
765 getgxf(pbto, peto, 60, 12, &nt, ecto, &Ctrl->CM4_DATA.gmdl[nout-1], tdto);
766
767 if (Ctrl->CM4_L.curr) {
768 bc[0] = bc[1] = bc[2] = 0.;
769 jpoloid(pbto, peto, 60, 12, ro, rm, nt, nyto, tdto, hq, hyto);
770 blsgen(nyto, nyto, 1, &bc[2], ecto, &hyto[nyto * 2]);
771 if (omdl)
772 tseardr(omdl, 2, ntay, nsto, nyto, frho, tsto, ecto, gcto_or);
773 else
774 tseardr(omdl, 2, ntay, nsto, nyto, frho, tsto, ecto, &gcto_mg[(((mmdl << 1) + 1) * 3 + 1) * 2736 - 27360]);
775
776 blsgen(nyto, nyto, 2, bc, ecto, hyto);
777 ltrans(1, 1, bc, rlgm, &jmdl[9]);
778 }
779 }
780
781 x = y = z = 0.0;
782 if (!Ctrl->CM4_L.curr) { /* Magnetic field */
783 for (k = 0; k < Ctrl->CM4_F.n_field_sources; k++) { /* Sum all field sources */
784 x += bmdl[Ctrl->CM4_F.field_sources[k]*3];
785 y += bmdl[Ctrl->CM4_F.field_sources[k]*3+1];
786 z += bmdl[Ctrl->CM4_F.field_sources[k]*3+2];
787 }
788 for (j = 0; j < Ctrl->CM4_F.n_field_components; j++) { /* Loop over vector field components */
789 h = 0.;
790 if (Ctrl->CM4_F.field_components[j] == 0) {
791 t = sqrt(x*x + y*y + z*z);
792 Ctrl->CM4_DATA.out_field[n*Ctrl->CM4_F.n_field_components+j] = t;
793 }
794 else if (Ctrl->CM4_F.field_components[j] == 1) {
795 h = sqrt(x*x + y*y);
796 Ctrl->CM4_DATA.out_field[n*Ctrl->CM4_F.n_field_components+j] = h;
797 }
798 else if (Ctrl->CM4_F.field_components[j] == 2)
799 Ctrl->CM4_DATA.out_field[n*Ctrl->CM4_F.n_field_components+j] = x;
800 else if (Ctrl->CM4_F.field_components[j] == 3)
801 Ctrl->CM4_DATA.out_field[n*Ctrl->CM4_F.n_field_components+j] = y;
802 else if (Ctrl->CM4_F.field_components[j] == 4)
803 Ctrl->CM4_DATA.out_field[n*Ctrl->CM4_F.n_field_components+j] = z;
804 else if (Ctrl->CM4_F.field_components[j] == 5)
805 Ctrl->CM4_DATA.out_field[n*Ctrl->CM4_F.n_field_components+j] = atan2(y,x) * R2D;
806 else if (Ctrl->CM4_F.field_components[j] == 6) {
807 if (!h) h = sqrt(x*x + y*y);
808 Ctrl->CM4_DATA.out_field[n*Ctrl->CM4_F.n_field_components+j] = atan2(z,h) * R2D;
809 }
810 }
811 }
812 else { /* Current density field (J) */
813 for (k = 0; k < Ctrl->CM4_L.n_curr_sources; k++) { /* Sum all current sources */
814 x += jmdl[Ctrl->CM4_L.curr_sources[k]*3];
815 y += jmdl[Ctrl->CM4_L.curr_sources[k]*3+1];
816 z += jmdl[Ctrl->CM4_L.curr_sources[k]*3+2];
817 }
818 for (j = 0; j < Ctrl->CM4_L.n_curr_components; j++) { /* Loop over current components */
819 if (Ctrl->CM4_L.curr_components[j] == 0) {
820 t = sqrt(x*x + y*y + z*z);
821 Ctrl->CM4_DATA.out_field[n*Ctrl->CM4_L.n_curr_components+j] = t;
822 }
823 else if (Ctrl->CM4_L.curr_components[j] == 1)
824 Ctrl->CM4_DATA.out_field[n*Ctrl->CM4_L.n_curr_components+j] = x;
825 else if (Ctrl->CM4_L.curr_components[j] == 2)
826 Ctrl->CM4_DATA.out_field[n*Ctrl->CM4_L.n_curr_components+j] = y;
827 else if (Ctrl->CM4_L.curr_components[j] == 3)
828 Ctrl->CM4_DATA.out_field[n*Ctrl->CM4_L.n_curr_components+j] = z;
829 }
830 }
831 }
832
833 clear_mem (mut, gpsq, gssq, gpmg, gsmg, hysq, epsq, essq, ecto, hyto, hq, ht, bkpo, ws, gamf, epmg,
834 esmg, hymg, f107x, pleg, rcur, gcto_or, gcto_mg);
835 return 0;
836 }
837
clear_mem(double * mut,double * gpsq,double * gssq,double * gpmg,double * gsmg,double * hysq,double * epsq,double * essq,double * ecto,double * hyto,double * hq,double * ht,double * bkpo,double * ws,double * gamf,double * epmg,double * esmg,double * hymg,double * f107x,double * pleg,double * rcur,double * gcto_or,double * gcto_mg)838 static void clear_mem (double *mut, double *gpsq, double *gssq, double *gpmg, double *gsmg, double *hysq, double *epsq, double *essq,
839 double *ecto, double *hyto, double *hq, double *ht, double *bkpo, double *ws, double *gamf, double *epmg,
840 double *esmg, double *hymg, double *f107x, double *pleg, double *rcur, double *gcto_or, double *gcto_mg) {
841 gmt_M_str_free (mut);
842 gmt_M_str_free (gpsq);
843 gmt_M_str_free (gssq);
844 gmt_M_str_free (gpmg);
845 gmt_M_str_free (gsmg);
846 gmt_M_str_free (hysq);
847 gmt_M_str_free (epsq);
848 gmt_M_str_free (essq);
849 gmt_M_str_free (ecto);
850 gmt_M_str_free (hyto);
851 gmt_M_str_free (hq);
852 gmt_M_str_free (ht);
853 gmt_M_str_free (bkpo);
854 gmt_M_str_free (ws);
855 gmt_M_str_free (gamf);
856 gmt_M_str_free (epmg);
857 gmt_M_str_free (esmg);
858 gmt_M_str_free (hymg);
859 gmt_M_str_free (f107x);
860 gmt_M_str_free (pleg);
861 gmt_M_str_free (rcur);
862 gmt_M_str_free (gcto_or);
863 gmt_M_str_free (gcto_mg);
864 }
865
ymdtomjd(int yearad,int month,int dayofmonth,int * mjd,int * dayofyear)866 void ymdtomjd(int yearad, int month, int dayofmonth, int *mjd, int *dayofyear) {
867 static int daysuptomonth[12] = { 0,31,59,90,120,151,181,212,243,273, 304,334 };
868 int remyear;
869
870 remyear = yearad - 1900;
871 if (remyear > 0) {
872 --remyear;
873 *mjd = remyear / 4 * 1461 + 15384;
874 remyear %= 4;
875 *dayofyear = daysuptomonth[month - 1] + dayofmonth;
876 if (month > 2) {
877 *dayofyear += I_DIM(remyear, 2);
878 }
879 *mjd = *mjd + remyear * 365 + *dayofyear;
880 } else {
881 *dayofyear = daysuptomonth[month - 1] + dayofmonth;
882 *mjd = *dayofyear + 15019;
883 }
884 }
885
ydtomjdx(int yearad,int dayofyear,int * mjd,int * month,int * dayofmonth,int * daysinmonth)886 void ydtomjdx(int yearad, int dayofyear, int *mjd, int *month, int *dayofmonth, int *daysinmonth) {
887 static int daysuptomonth[12] = { 0,31,59,90,120,151,181,212,243,273,304,334 };
888 int j, leapday, leapcum, remyear;
889
890 remyear = yearad - 1900;
891 if (remyear > 0) {
892 --remyear;
893 *mjd = remyear / 4 * 1461 + 15384;
894 remyear %= 4;
895 *mjd = *mjd + remyear * 365 + dayofyear;
896 leapday = I_DIM(remyear, 2);
897 } else {
898 *mjd = dayofyear + 15019;
899 leapday = 0;
900 }
901 for (j = 12; j >= 1; --j) {
902 leapcum = MIN(1,I_DIM(j, 2)) * leapday;
903 if (daysuptomonth[j - 1] + leapcum < dayofyear) {
904 *month = j;
905 *dayofmonth = dayofyear - daysuptomonth[j - 1] - leapcum;
906 break;
907 }
908 }
909
910 daysinmonth[0] = 31;
911 daysinmonth[1] = leapday + 28;
912 daysinmonth[2] = 31;
913 daysinmonth[3] = 30;
914 daysinmonth[4] = 31;
915 daysinmonth[5] = 30;
916 daysinmonth[6] = 31;
917 daysinmonth[7] = 31;
918 daysinmonth[8] = 30;
919 daysinmonth[9] = 31;
920 daysinmonth[10] = 30;
921 daysinmonth[11] = 31;
922 }
923
intdst(int mjdl,int mjdh,int mjdy,int msec,double * dstx,int * cerr)924 double intdst(int mjdl, int mjdh, int mjdy, int msec, double *dstx, int *cerr) {
925 int hbot, mjdt, jbot, mobs, htop, jtop, hour;
926 double ttop, dst;
927
928 hour = msec / 3600000;
929 mjdt = mjdy + hour / 24;
930 hour = hour % 24 + 1;
931 mobs = msec % 3600000;
932 if (mobs <= 1800000) {
933 ttop = (double) (mobs + 1800000) / 3.6e6;
934 jtop = mjdt;
935 if (hour > 1) {
936 jbot = mjdt;
937 htop = hour;
938 hbot = hour - 1;
939 } else {
940 jbot = mjdt - 1;
941 htop = 1;
942 hbot = 24;
943 }
944 } else {
945 ttop = (double) (mobs - 1800000) / 3.6e6;
946 if (hour < 24) {
947 jtop = jbot = mjdt;
948 htop = hour + 1;
949 hbot = hour;
950 } else {
951 jtop = mjdt + 1;
952 jbot = mjdt;
953 htop = 1;
954 hbot = 24;
955 }
956 }
957 if (jbot < mjdl || jtop > mjdh) {
958 *cerr = 50;
959 fprintf (stderr, "INTDST -- Error: T (%d; %d) LIES OUTSIDE OF DST TABLE TIME SPAN [%d; %d] -- ABORT\n", jbot, jtop, mjdl, mjdh);
960 dst = -1e12;
961 }
962 else
963 dst = ttop * dstx[(jtop - mjdl)*24 + (htop-1)] + (1. - ttop) * dstx[(jbot - mjdl)*24 + (hbot-1)];
964
965 return (dst);
966 }
967
intf107(int iyrl,int imol,int iyrh,int imoh,int iyr,int imon,int idom,int * idim,int msec,double * f107x,int * cerr)968 double intf107(int iyrl, int imol, int iyrh, int imoh, int iyr, int imon, int idom, int *idim, int msec, double *f107x, int *cerr) {
969
970 int mbot, ybot, mtop, ytop;
971 double tadd, thlf, tobs, tint, ttop, f107;
972
973 /* Parameter adjustments */
974 --idim;
975
976 /* Function Body */
977 thlf = (double) idim[imon] * .5;
978 tobs = (double) (idom - 1) + (double) (msec) / 8.64e7;
979 if (tobs <= thlf) {
980 if (imon > 1) {
981 tadd = (double) idim[imon - 1] * .5;
982 tobs += tadd;
983 tint = thlf + tadd;
984 ttop = tobs / tint;
985 ytop = iyr;
986 ybot = iyr;
987 mtop = imon;
988 mbot = imon - 1;
989 } else {
990 tobs += 15.5;
991 tint = thlf + 15.5;
992 ttop = tobs / tint;
993 ytop = iyr;
994 ybot = iyr - 1;
995 mtop = 1;
996 mbot = 12;
997 }
998 } else {
999 if (imon < 12) {
1000 tadd = (double) idim[imon + 1] * .5;
1001 tobs -= thlf;
1002 tint = thlf + tadd;
1003 ttop = tobs / tint;
1004 ytop = iyr;
1005 ybot = iyr;
1006 mtop = imon + 1;
1007 mbot = imon;
1008 } else {
1009 tobs += -15.5;
1010 tint = thlf + 15.5;
1011 ttop = tobs / tint;
1012 ytop = iyr + 1;
1013 ybot = iyr;
1014 mtop = 1;
1015 mbot = 12;
1016 }
1017 }
1018 if (ybot < iyrl || ytop > iyrh || (ybot == iyrl && mbot < imol) || (ytop == iyrh && mtop > imoh)) {
1019 fprintf (stderr, "SUBROUTINE INTF107 -- ERROR CODE 50 -- T LIES OUTSIDE OF F10.7 TABLE TIME SPAN -- ABORT\n");
1020 f107 = -1.;
1021 *cerr = 50;
1022 }
1023 else
1024 f107 = ttop * f107x[(ytop - iyrl)*12 + mtop-1] + (1. - ttop) * f107x[(ybot - iyrl)*12 + mbot-1];
1025
1026 return (f107);
1027 }
1028
getmut2(double thenmp,double phinmp,int iyear,int iday,int msec)1029 double getmut2(double thenmp, double phinmp, int iyear, int iday, int msec) {
1030 double gst, cth0, the0, phi0, sth0, eadj, sdec, cphd, sphd, cths, thes, phis, eopp, gmts, sths, slong, srasn;
1031
1032 the0 = thenmp * D2R;
1033 phi0 = phinmp * D2R;
1034 sincos(the0, &sth0, &cth0);
1035 gmts = (double) (msec) * 1e-3;
1036 sun2(iyear, iday, gmts, &gst, &slong, &srasn, &sdec);
1037 thes = (90. - sdec) * D2R;
1038 sincos(thes, &sths, &cths);
1039 phis = (srasn - gst) * D2R;
1040 sincos(phis - phi0, &sphd, &cphd);
1041 eopp = sths * sphd;
1042 eadj = cth0 * sths * cphd - sth0 * cths;
1043 return (12. - atan2(eopp, eadj) * R2D / 15.);
1044 }
1045
sun2(int iyr,int iday,double secs,double * gst,double * slong,double * srasn,double * sdec)1046 void sun2(int iyr, int iday, double secs, double *gst, double *slong, double *srasn, double *sdec) {
1047 double d__1, g, t, dj, vl, slp, fday, sined, obliq, cosined;
1048
1049 if (iyr < 1901 || iyr > 2099) {
1050 *gst = 0.;
1051 *slong = 0.;
1052 *srasn = 0.;
1053 *sdec = 0.;
1054 } else {
1055 fday = secs / 86400.;
1056 dj = (double) (iyr - 1900) * 365. + (double) ((iyr - 1901) / 4) + (double) (iday) + fday - .5;
1057 t = dj / 36525.;
1058 d__1 = dj * .9856473354 + 279.696678;
1059 vl = d_mod(d__1, 360);
1060 d__1 = dj * .9856473354 + 279.690983 + fday * 360. + 180.;
1061 *gst = d_mod(d__1, 360);
1062 d__1 = dj * .985600267 + 358.475845;
1063 g = d_mod(d__1, 360) * D2R;
1064 *slong = vl + (1.91946 - t * .004789) * sin(g) + sin(g * 2.) * .020094;
1065 obliq = (23.45229 - t * .0130125) * D2R;
1066 slp = (*slong - .005686) * D2R;
1067 sined = sin(obliq) * sin(slp);
1068 cosined = sqrt(1. - sined * sined);
1069 *sdec = R2D * atan(sined / cosined);
1070 *srasn = 180. - R2D * atan2(1. / tan(obliq) * sined / cosined, -cos(slp) / cosined);
1071 }
1072 }
1073
rmerge_(double * rmrg,double * rmlt)1074 void rmerge_(double *rmrg, double *rmlt) {
1075 double r1, r2, r3;
1076
1077 r1 = rmrg[0];
1078 r2 = rmrg[1];
1079 r3 = rmrg[2];
1080 rmrg[0] = r1 * rmlt[0] + r2 * rmlt[3] + r3 * rmlt[6];
1081 rmrg[1] = r1 * rmlt[1] + r2 * rmlt[4] + r3 * rmlt[7];
1082 rmrg[2] = r1 * rmlt[2] + r2 * rmlt[5] + r3 * rmlt[8];
1083 r1 = rmrg[3];
1084 r2 = rmrg[4];
1085 r3 = rmrg[5];
1086 rmrg[3] = r1 * rmlt[0] + r2 * rmlt[3] + r3 * rmlt[6];
1087 rmrg[4] = r1 * rmlt[1] + r2 * rmlt[4] + r3 * rmlt[7];
1088 rmrg[5] = r1 * rmlt[2] + r2 * rmlt[5] + r3 * rmlt[8];
1089 r1 = rmrg[6];
1090 r2 = rmrg[7];
1091 r3 = rmrg[8];
1092 rmrg[6] = r1 * rmlt[0] + r2 * rmlt[3] + r3 * rmlt[6];
1093 rmrg[7] = r1 * rmlt[1] + r2 * rmlt[4] + r3 * rmlt[7];
1094 rmrg[8] = r1 * rmlt[2] + r2 * rmlt[5] + r3 * rmlt[8];
1095 }
1096
tsearad(int full,int ks,int kr,int ns,int ng,double f,double * t,double * e,double * g)1097 void tsearad(int full, int ks, int kr, int ns, int ng, double f, double *t, double *e, double *g) {
1098 int i, j, k;
1099 double s, z;
1100
1101 /* Parameter adjustments */
1102 g -= (1 + ng * (1 + ns));
1103 --t;
1104
1105 /* Function Body */
1106 memset(e, 0, ng * sizeof(double));
1107 j = 1;
1108 s = 1.;
1109 r8vlinkt(1, 1, ng, s, &g[(j + ns) * ng + 1], &e[0]);
1110 for (i = 1; i <= ks; ++i) {
1111 ++j;
1112 s = t[i + 1];
1113 r8vlinkt(1, 1, ng, s, &g[(j + ns) * ng + 1], &e[0]);
1114 if (full) {
1115 ++j;
1116 s = t[i + ks + 2];
1117 r8vlinkt(1, 1, ng, s, &g[(j + ns) * ng + 1], &e[0]);
1118 }
1119 }
1120 z = 1.;
1121 for (k = 1; k <= kr; ++k) {
1122 j = 1;
1123 z = z * f / (double) k;
1124 r8vlinkt(1, 1, ng, z, &g[(j + (k + 1) * ns) * ng + 1], &e[0]);
1125 for (i = 1; i <= ks; ++i) {
1126 ++j;
1127 s = t[i + 1] * z;
1128 r8vlinkt(1, 1, ng, s, &g[(j + (k + 1) * ns) * ng + 1], &e[0]);
1129 if (full) {
1130 ++j;
1131 s = t[i + ks + 2] * z;
1132 r8vlinkt(1, 1, ng, s, &g[(j + (k + 1) * ns) * ng + 1], &e[0]);
1133 }
1134 }
1135 }
1136 }
1137
tseardr(int full,int ks,int kr,int ns,int ng,double f,double * t,double * e,double * g)1138 void tseardr(int full, int ks, int kr, int ns, int ng, double f, double *t, double *e, double *g) {
1139 int i, j, k;
1140 double s, z;
1141
1142 /* Parameter adjustments */
1143 g -= (1 + ng * (1 + ns));
1144 --t;
1145
1146 /* Function Body */
1147 memset(e, 0, ng * sizeof(double));
1148 z = 1.;
1149 for (k = 1; k <= kr; ++k) {
1150 j = 1;
1151 r8vlinkt(1, 1, ng, z, &g[(j + (k + 1) * ns) * ng + 1], &e[0]);
1152 for (i = 1; i <= ks; ++i) {
1153 ++j;
1154 s = t[i + 1] * z;
1155 r8vlinkt(1, 1, ng, s, &g[(j + (k + 1) * ns) * ng + 1], &e[0]);
1156 if (full) {
1157 ++j;
1158 s = t[i + ks + 2] * z;
1159 r8vlinkt(1, 1, ng, s, &g[(j + (k + 1) * ns) * ng + 1], &e[0]);
1160 }
1161 }
1162 z = z * f / (double) k;
1163 }
1164 }
1165
mseason(int ks,int ns,int ng,double d,double * t,double * e,double * g)1166 void mseason(int ks, int ns, int ng, double d, double *t, double *e, double *g) {
1167 int i, j;
1168 double s;
1169
1170 /* Parameter adjustments */
1171 g -= (1 + ng * (1 + ns));
1172
1173 /* Function Body */
1174 memset(e, 0, ng * sizeof(double));
1175 j = 1;
1176 s = 1.;
1177 r8vlinkt(1, 1, ng, s, &g[(j + ns) * ng + 1], &e[0]);
1178 r8vlinkt(1, 1, ng, d, &g[(j + (ns << 1)) * ng + 1], &e[0]);
1179 for (i = 1; i <= ks; ++i) {
1180 ++j;
1181 s = t[i];
1182 r8vlinkt(1, 1, ng, s, &g[(j + ns) * ng + 1], &e[0]);
1183 s *= d;
1184 r8vlinkt(1, 1, ng, s, &g[(j + (ns << 1)) * ng + 1], &e[0]);
1185 ++j;
1186 s = t[i + ks + 1];
1187 r8vlinkt(1, 1, ng, s, &g[(j + ns) * ng + 1], &e[0]);
1188 s *= d;
1189 r8vlinkt(1, 1, ng, s, &g[(j + (ns << 1)) * ng + 1], &e[0]);
1190 }
1191 }
1192
iseason(int ks,int ns,int ng,double f,double * t,double * e,double * g)1193 void iseason(int ks, int ns, int ng, double f, double *t, double *e, double *g) {
1194 int i, j;
1195 double s;
1196 gmt_M_unused(ns);
1197
1198 memset(e, 0, ng * sizeof(double));
1199 j = 0;
1200 r8vlinkt(1, 1, ng, f, &g[j * ng], &e[0]);
1201 for (i = 1; i <= ks; ++i) {
1202 ++j;
1203 s = f * t[i];
1204 r8vlinkt(1, 1, ng, s, &g[j * ng], &e[0]);
1205 ++j;
1206 s = f * t[i + ks + 1];
1207 r8vlinkt(1, 1, ng, s, &g[j * ng], &e[0]);
1208 }
1209 }
1210
mpotent(int nmax,int mmax,int nd,int nz,double cphi,double sphi,double * d,double * z)1211 void mpotent(int nmax, int mmax, int nd, int nz, double cphi, double sphi, double *d, double *z) {
1212 int m, n, id, iz, nd2;
1213
1214 /* Parameter adjustments */
1215 z -= (1 + nz);
1216 --d;
1217
1218 /* Function Body */
1219 nd2 = nd + nd;
1220 id = iz = 0;
1221 for (n = 1; n <= nmax; ++n) {
1222 ++id;
1223 ++iz;
1224 z[iz + nz] = d[id] * cphi;
1225 z[iz + (nz << 1)] = d[nd + id] * cphi;
1226 z[iz + nz * 3] = d[nd2 + id] * cphi;
1227 ++iz;
1228 z[iz + nz] = d[id] * sphi;
1229 z[iz + (nz << 1)] = d[nd + id] * sphi;
1230 z[iz + nz * 3] = d[nd2 + id] * sphi;
1231 for (m = 1; m <= MIN(n,mmax); ++m) {
1232 id += 2;
1233 ++iz;
1234 z[iz + nz] = d[id - 1] * cphi + d[id] * sphi;
1235 z[iz + (nz << 1)] = d[nd + id - 1] * cphi + d[nd + id] * sphi;
1236 z[iz + nz * 3] = d[nd2 + id - 1] * cphi + d[nd2 + id] * sphi;
1237 ++iz;
1238 z[iz + nz] = d[id] * cphi - d[id - 1] * sphi;
1239 z[iz + (nz << 1)] = d[nd + id] * cphi - d[nd + id - 1] * sphi;
1240 z[iz + nz * 3] = d[nd2 + id] * cphi - d[nd2 + id - 1] * sphi;
1241 ++iz;
1242 z[iz + nz] = d[id - 1] * cphi - d[id] * sphi;
1243 z[iz + (nz << 1)] = d[nd + id - 1] * cphi - d[nd + id] * sphi;
1244 z[iz + nz * 3] = d[nd2 + id - 1] * cphi - d[nd2 + id] * sphi;
1245 ++iz;
1246 z[iz + nz] = d[id] * cphi + d[id - 1] * sphi;
1247 z[iz + (nz << 1)] = d[nd + id] * cphi + d[nd + id - 1] * sphi;
1248 z[iz + nz * 3] = d[nd2 + id] * cphi + d[nd2 + id - 1] * sphi;
1249 }
1250 }
1251 }
1252
jtabove(int pmin,int pmax,int nmax,int mmax,double r,double rref,int nz,double * z)1253 void jtabove(int pmin, int pmax, int nmax, int mmax, double r, double rref, int nz, double *z) {
1254 /* Initialized data */
1255
1256 static double fgeo = 7.95774715459478e-4;
1257
1258 /* Local variables */
1259 int m, n, p, iz;
1260 double ffac, fcur, fpsi, frmu, frpw, ztemp;
1261
1262 /* Parameter adjustments */
1263 z -= (1 + nz);
1264
1265 /* Function Body */
1266 frmu = rref / r;
1267 iz = 0;
1268 for (p = pmin; p <= pmax; ++p) {
1269 frpw = fgeo;
1270 for (n = 1; n <= nmax; ++n) {
1271 ffac = frpw * (double) (n + n + 1);
1272 fcur = ffac / (double) (n + 1);
1273 fpsi = -(rref) * ffac / (double) (n * (n + 1));
1274 ++iz;
1275 ztemp = z[iz + nz];
1276 z[iz + nz] = -fcur * z[iz + (nz << 1)];
1277 z[iz + (nz << 1)] = fcur * ztemp;
1278 z[iz + nz * 3] = fpsi * z[iz + nz * 3];
1279 ++iz;
1280 ztemp = z[iz + nz];
1281 z[iz + nz] = -fcur * z[iz + (nz << 1)];
1282 z[iz + (nz << 1)] = fcur * ztemp;
1283 z[iz + nz * 3] = fpsi * z[iz + nz * 3];
1284 for (m = 1; m <= MIN(n,mmax); ++m) {
1285 ++iz;
1286 ztemp = z[iz + nz];
1287 z[iz + nz] = -fcur * z[iz + (nz << 1)];
1288 z[iz + (nz << 1)] = fcur * ztemp;
1289 z[iz + nz * 3] = fpsi * z[iz + nz * 3];
1290 ++iz;
1291 ztemp = z[iz + nz];
1292 z[iz + nz] = -fcur * z[iz + (nz << 1)];
1293 z[iz + (nz << 1)] = fcur * ztemp;
1294 z[iz + nz * 3] = fpsi * z[iz + nz * 3];
1295 ++iz;
1296 ztemp = z[iz + nz];
1297 z[iz + nz] = -fcur * z[iz + (nz << 1)];
1298 z[iz + (nz << 1)] = fcur * ztemp;
1299 z[iz + nz * 3] = fpsi * z[iz + nz * 3];
1300 ++iz;
1301 ztemp = z[iz + nz];
1302 z[iz + nz] = -fcur * z[iz + (nz << 1)];
1303 z[iz + (nz << 1)] = fcur * ztemp;
1304 z[iz + nz * 3] = fpsi * z[iz + nz * 3];
1305 }
1306 frpw *= frmu;
1307 }
1308 }
1309 }
1310
jtbelow(int pmin,int pmax,int nmax,int mmax,double r,double rref,int nz,double * z)1311 void jtbelow(int pmin, int pmax, int nmax, int mmax, double r, double rref, int nz, double *z) {
1312 /* Initialized data */
1313
1314 static double fgeo = 7.95774715459478e-4;
1315
1316 /* Local variables */
1317 int m, n, p, iz;
1318 double ffac, frbg, fcur, fpsi, frmu, frpw, ztemp;
1319
1320 /* Parameter adjustments */
1321 z -= (1 + nz);
1322
1323 /* Function Body */
1324 frmu = r / rref;
1325 frbg = fgeo * (frmu * frmu * frmu);
1326 iz = 0;
1327 for (p = pmin; p <= pmax; ++p) {
1328 frpw = frbg;
1329 for (n = 1; n <= nmax; ++n) {
1330 ffac = frpw * (double) (n + n + 1);
1331 fcur = ffac / (double) n;
1332 fpsi = -(rref) * ffac / (double) (n * (n + 1));
1333 ++iz;
1334 ztemp = z[iz + nz];
1335 z[iz + nz] = fcur * z[iz + (nz << 1)];
1336 z[iz + (nz << 1)] = -fcur * ztemp;
1337 z[iz + nz * 3] = fpsi * z[iz + nz * 3];
1338 ++iz;
1339 ztemp = z[iz + nz];
1340 z[iz + nz] = fcur * z[iz + (nz << 1)];
1341 z[iz + (nz << 1)] = -fcur * ztemp;
1342 z[iz + nz * 3] = fpsi * z[iz + nz * 3];
1343 for (m = 1; m <= MIN(n,mmax); ++m) {
1344 ++iz;
1345 ztemp = z[iz + nz];
1346 z[iz + nz] = fcur * z[iz + (nz << 1)];
1347 z[iz + (nz << 1)] = -fcur * ztemp;
1348 z[iz + nz * 3] = fpsi * z[iz + nz * 3];
1349 ++iz;
1350 ztemp = z[iz + nz];
1351 z[iz + nz] = fcur * z[iz + (nz << 1)];
1352 z[iz + (nz << 1)] = -fcur * ztemp;
1353 z[iz + nz * 3] = fpsi * z[iz + nz * 3];
1354 ++iz;
1355 ztemp = z[iz + nz];
1356 z[iz + nz] = fcur * z[iz + (nz << 1)];
1357 z[iz + (nz << 1)] = -fcur * ztemp;
1358 z[iz + nz * 3] = fpsi * z[iz + nz * 3];
1359 ++iz;
1360 ztemp = z[iz + nz];
1361 z[iz + nz] = fcur * z[iz + (nz << 1)];
1362 z[iz + (nz << 1)] = -fcur * ztemp;
1363 z[iz + nz * 3] = fpsi * z[iz + nz * 3];
1364 }
1365 frpw *= frmu;
1366 }
1367 }
1368 }
1369
jtbcont(int pmin,int pmax,int nmax,int mmax,double rold,double rnew,int nz,double * z)1370 void jtbcont(int pmin, int pmax, int nmax, int mmax, double rold, double rnew, int nz, double *z) {
1371 int m, n, p, iz, ii;
1372 double frbg, fcur, fpsi, frmu;
1373
1374 /* Parameter adjustments */
1375 z -= (1 + nz);
1376
1377 /* Function Body */
1378 frmu = rold / rnew;
1379 frbg = frmu * frmu;
1380 iz = 0;
1381 for (p = pmin; p <= pmax; ++p) {
1382 fpsi = frbg;
1383 for (n = 1; n <= nmax; ++n) {
1384 fcur = fpsi * frmu;
1385 ++iz;
1386 ii = iz + nz; z[ii] = fcur * z[ii];
1387 ii = iz + (nz << 1); z[ii] = fcur * z[ii];
1388 ii = iz + nz * 3; z[ii] = fpsi * z[ii];
1389 ++iz;
1390 ii = iz + nz; z[ii] = fcur * z[ii];
1391 ii = iz + (nz << 1); z[ii] = fcur * z[ii];
1392 ii = iz + nz * 3; z[ii] = fpsi * z[ii];
1393 for (m = 1; m <= MIN(n,mmax); ++m) {
1394 ++iz;
1395 ii = iz + nz; z[ii] = fcur * z[ii];
1396 ii = iz + (nz << 1); z[ii] = fcur * z[ii];
1397 ii = iz + nz * 3; z[ii] = fpsi * z[ii];
1398 ++iz;
1399 ii = iz + nz; z[ii] = fcur * z[ii];
1400 ii = iz + (nz << 1); z[ii] = fcur * z[ii];
1401 ii = iz + nz * 3; z[ii] = fpsi * z[ii];
1402 ++iz;
1403 ii = iz + nz; z[ii] = fcur * z[ii];
1404 ii = iz + (nz << 1); z[ii] = fcur * z[ii];
1405 ii = iz + nz * 3; z[ii] = fpsi * z[ii];
1406 ++iz;
1407 ii = iz + nz; z[ii] = fcur * z[ii];
1408 ii = iz + (nz << 1); z[ii] = fcur * z[ii];
1409 ii = iz + nz * 3; z[ii] = fpsi * z[ii];
1410 }
1411 fpsi *= frmu;
1412 }
1413 }
1414 }
1415
mstream(int nmax,int mmax,int nd,int nz,double cphi,double sphi,double faor,double * d,double * z)1416 void mstream(int nmax, int mmax, int nd, int nz, double cphi, double sphi, double faor, double *d, double *z) {
1417 int m, n, id, iz;
1418
1419 /* Parameter adjustments */
1420 z -= (1 + nz);
1421 --d;
1422
1423 /* Function Body */
1424 id = iz = 0;
1425 for (n = 1; n <= nmax; ++n) {
1426 ++id;
1427 ++iz;
1428 z[iz + nz] = faor * d[nd + id] * cphi;
1429 z[iz + (nz << 1)] = -(faor) * d[id] * cphi;
1430 ++iz;
1431 z[iz + nz] = faor * d[nd + id] * sphi;
1432 z[iz + (nz << 1)] = -(faor) * d[id] * sphi;
1433 for (m = 1; m <= MIN(n,mmax); ++m) {
1434 id += 2;
1435 ++iz;
1436 z[iz + nz] = faor * (d[nd + id - 1] * cphi + d[nd + id] * sphi);
1437 z[iz + (nz << 1)] = -(faor) * (d[id - 1] * cphi + d[id] * sphi);
1438 ++iz;
1439 z[iz + nz] = faor * (d[nd + id] * cphi - d[nd + id - 1] * sphi);
1440 z[iz + (nz << 1)] = -(faor) * (d[id] * cphi - d[id - 1] * sphi);
1441 ++iz;
1442 z[iz + nz] = faor * (d[nd + id - 1] * cphi - d[nd + id] * sphi);
1443 z[iz + (nz << 1)] = -(faor) * (d[id - 1] * cphi - d[id] * sphi);
1444 ++iz;
1445 z[iz + nz] = faor * (d[nd + id] * cphi + d[nd + id - 1] * sphi);
1446 z[iz + (nz << 1)] = -(faor) * (d[id] * cphi + d[id - 1] * sphi);
1447 }
1448 }
1449 }
1450
jpoloid(int pmin,int pmax,int nmax,int mmax,double r,double rm,int nd,int nz,double * t,double * d,double * z)1451 void jpoloid(int pmin, int pmax, int nmax, int mmax, double r, double rm, int nd, int nz, double *t, double *d, double *z) {
1452 /* Initialized data */
1453
1454 static double u0 = .0012566370614359158;
1455
1456 /* Local variables */
1457 int m, n, p, id, iz, nd2;
1458 double fdeg, fhtj, frtj, cosp, sinp, ztemp;
1459
1460 /* Parameter adjustments */
1461 z -= (1 + nz);
1462 --d;
1463
1464 /* Function Body */
1465 nd2 = nd + nd;
1466 fhtj = 1. / rm / u0;
1467 frtj = rm / (r * r) / u0;
1468 iz = 0;
1469 for (p = pmin; p <= pmax; ++p) {
1470 id = 0;
1471 cosp = t[p];
1472 sinp = t[p + 1 + pmax];
1473 for (n = 1; n <= nmax; ++n) {
1474 fdeg = frtj * (double) n;
1475 ++id;
1476 ++iz;
1477 ztemp = z[iz + nz];
1478 z[iz + nz] = fhtj * z[iz + (nz << 1)];
1479 z[iz + (nz << 1)] = -fhtj * ztemp;
1480 z[iz + nz * 3] = fdeg * d[nd2 + id] * cosp;
1481 ++iz;
1482 ztemp = z[iz + nz];
1483 z[iz + nz] = fhtj * z[iz + (nz << 1)];
1484 z[iz + (nz << 1)] = -fhtj * ztemp;
1485 z[iz + nz * 3] = fdeg * d[nd2 + id] * sinp;
1486 for (m = 1; m <= MIN(n,mmax); ++m) {
1487 id += 2;
1488 ++iz;
1489 ztemp = z[iz + nz];
1490 z[iz + nz] = fhtj * z[iz + (nz << 1)];
1491 z[iz + (nz << 1)] = -fhtj * ztemp;
1492 z[iz + nz * 3] = fdeg * (d[nd2 + id - 1] * cosp + d[nd2 + id] * sinp);
1493 ++iz;
1494 ztemp = z[iz + nz];
1495 z[iz + nz] = fhtj * z[iz + (nz << 1)];
1496 z[iz + (nz << 1)] = -fhtj * ztemp;
1497 z[iz + nz * 3] = fdeg * (d[nd2 + id] * cosp - d[nd2 + id - 1] * sinp);
1498 ++iz;
1499 ztemp = z[iz + nz];
1500 z[iz + nz] = fhtj * z[iz + (nz << 1)];
1501 z[iz + (nz << 1)] = -fhtj * ztemp;
1502 z[iz + nz * 3] = fdeg * (d[nd2 + id - 1] * cosp - d[nd2 + id] * sinp);
1503 ++iz;
1504 ztemp = z[iz + nz];
1505 z[iz + nz] = fhtj * z[iz + (nz << 1)];
1506 z[iz + (nz << 1)] = -fhtj * ztemp;
1507 z[iz + nz * 3] = fdeg * (d[nd2 + id] * cosp + d[nd2 + id - 1] * sinp);
1508 }
1509 }
1510 }
1511 }
1512
blsgen(int nc,int nd,int ni,double * b,double * c,double * dldc)1513 void blsgen(int nc, int nd, int ni, double *b, double *c, double *dldc) {
1514 int j;
1515
1516 /* Parameter adjustments */
1517 dldc -= (1 + nd);
1518 --b;
1519
1520 for (j = 1; j <= ni; ++j)
1521 b[j] += r8sdot(1, 1, nc, &dldc[j * nd + 1], &c[0]);
1522 }
1523
getgmf(int nder,int ns,double * ep,double * tm,double * b,double * c,double * g,int * h,int * o,double * p)1524 void getgmf(int nder, int ns, double *ep, double *tm, double *b, double *c, double *g, int *h, int *o, double *p) {
1525 /* Initialized data */
1526
1527 static int null = 0;
1528
1529 /* Local variables */
1530 int d, i, j, k, m, n, y, ic, la, lb, jc, ig, nb, ik, ni, im, nk, jm, no, js, nw;
1531 double gd, gm, bsum;
1532
1533 /* Parameter adjustments */
1534 --p;
1535 --o;
1536 --h;
1537 --g;
1538 --c;
1539 --b;
1540
1541 /* Function Body */
1542 ic = 1;
1543 ik = 1;
1544 for (i = 1; i <= ns; ++i) {
1545 g[i] = c[ic];
1546 ig = i;
1547 for (j = 1; j <= nder; ++j) {
1548 ig += ns;
1549 g[ig] = 0.;
1550 }
1551 n = o[i];
1552 if (n > 0) {
1553 k = h[i];
1554 if (*tm < p[ik] || *tm > p[ik + k + 1]) {
1555 fprintf(stderr, "SUBROUTINE GETGMF -- ERROR CODE 56 -- T LIES OUTSIDE OF KNOT DOMAIN -- ABORT\n");
1556 return;
1557 }
1558 m = n + 1;
1559 nk = k + 2;
1560 nb = n + k;
1561 ni = nb + 1;
1562 no = ni + 1;
1563 nw = no + no;
1564 r8slt(1, nk, *ep, &p[ik], &y);
1565 la = MIN(nk,y+1);
1566 r8slt(1, nk, *tm, &p[ik], &y);
1567 lb = MIN(nk,y+1);
1568 memset(&b[1], 0, nw * sizeof(double));
1569 dbspln_(&la, ep, &m, &null, &k, &p[ik], &b[la - 1], &b[nw + 1]);
1570 dbspln_(&lb, tm, &m, &null, &k, &p[ik], &b[no + lb - 1], &b[nw + 1]);
1571 r8vsub(1, no+1, 1, no, &b[1], &b[1], &b[1]);
1572 bsum = 0.;
1573 gm = 0.;
1574 jc = ic + nb;
1575 for (j = ni; j >= 2; --j) {
1576 bsum += b[j];
1577 im = MIN(j,nk);
1578 jm = MAX(j-n,1);
1579 gm += (p[ik + im - 1] - p[ik + jm - 1]) * bsum * c[jc];
1580 --jc;
1581 }
1582 g[i] += gm / (double) n;
1583 js = ic + lb - 1;
1584 ig = i;
1585 for (d = 0; d <= nder - 1; ++d) {
1586 ig += ns;
1587 dbspln_(&lb, tm, &n, &d, &k, &p[ik], &b[1], &b[nw + 1]);
1588 gd = 0.;
1589 jc = js;
1590 for (j = 1; j <= n; ++j) {
1591 gd += b[j] * c[jc];
1592 ++jc;
1593 }
1594 g[ig] = gd;
1595 }
1596 ic += nb;
1597 ik += nk;
1598 }
1599 ++ic;
1600 }
1601 }
1602
dbspln_(int * l,double * t,int * n,int * d__,int * k,double * x,double * b,double * w)1603 void dbspln_(int *l, double *t, int *n, int * d__, int *k, double *x, double *b, double *w) {
1604 /* Local variables */
1605 int i__, j, m, q, ik, jk, iw, lenb, posb, posw, posx;
1606 double fn, difi, difj, delx, temp;
1607
1608 /* Parameter adjustments */
1609 --w;
1610 --b;
1611 --x;
1612
1613 /* Function Body */
1614 q = *n - *d__;
1615 if (q == 1)
1616 b[1] = 1.;
1617 else {
1618 m = *l;
1619 ik = MIN(m,*k + 2);
1620 jk = MAX(m-1,1);
1621 difi = x[ik] - *t;
1622 delx = x[ik] - x[jk];
1623 if (delx == 0.)
1624 b[q] = 0.;
1625 else
1626 b[q] = 1. / delx;
1627
1628 posb = q - 1;
1629 for (j = 2; j <= q; ++j) {
1630 jk = MAX(m-j,1);
1631 temp = difi * b[posb + 1];
1632 delx = x[ik] - x[jk];
1633 if (delx == 0.)
1634 temp = 0.;
1635 else if (j < *n)
1636 temp /= delx;
1637
1638 b[posb] = temp;
1639 --posb;
1640 }
1641 b[q + 1] = 0.;
1642 ++m;
1643 for (i__ = 2; i__ <= q; ++i__) {
1644 ik = MIN(m,*k + 2);
1645 difi = x[ik] - *t;
1646 posb = q;
1647 for (j = i__; j <= q; ++j) {
1648 jk = MAX(m-j,1);
1649 difj = *t - x[jk];
1650 temp = difj * b[posb] + difi * b[posb + 1];
1651 delx = x[ik] - x[jk];
1652 if (delx == 0.)
1653 temp = 0.;
1654 else if (j < *n)
1655 temp /= delx;
1656
1657 b[posb] = temp;
1658 --posb;
1659 }
1660 ++m;
1661 }
1662 }
1663 posx = *l + *n - 1;
1664 posw = *d__ + *n;
1665 for (i__ = 1; i__ <= *n; ++i__) {
1666 lenb = MIN(q, posw - *d__);
1667 memset(&w[1], 0, posw * sizeof(double));
1668 r8vgathp(1, 1, *d__ + 1, lenb, &b[1], &w[1]);
1669 for (j = 1; j <= *d__; ++j) {
1670 ik = posx;
1671 jk = posx - q - j;
1672 iw = posw;
1673 fn = (double) (q + j - 1);
1674 for (m = j; m <= *d__; ++m) {
1675 if (j < *d__) {
1676 delx = x[MAX(1,MIN(ik,*k + 2))] - x[MAX(jk,1)];
1677 if (delx == 0.)
1678 w[iw] = 0.;
1679 else
1680 w[iw] = fn * (w[iw - 1] - w[iw]) / delx;
1681
1682 } else
1683 w[iw] = fn * (w[iw - 1] - w[iw]);
1684
1685 --ik;
1686 --jk;
1687 --iw;
1688 }
1689 }
1690 --posx;
1691 --posw;
1692 }
1693 r8vgathp(*d__ + 1, 1, 1, *n, &w[1], &b[1]);
1694 }
1695
getgxf(int pmin,int pmax,int nmax,int mmax,int * ng,double * e,double * g,double * t)1696 void getgxf(int pmin, int pmax, int nmax, int mmax, int *ng, double *e, double *g, double *t) {
1697 int m, n, p, ie, ig;
1698 double cosp, sinp;
1699
1700 /* Function Body */
1701 memset(g, 0, *ng * sizeof(double));
1702 ie = 0;
1703 for (p = pmin; p <= pmax; ++p) {
1704 ig = 0;
1705 cosp = t[p];
1706 sinp = t[p + 1 + pmax];
1707 for (n = 1; n <= nmax; ++n) {
1708 g[ig++] += (e[ie] * cosp + e[ie + 1] * sinp);
1709 ie += 2;
1710 for (m = 1; m <= MIN(n,mmax); ++m) {
1711 g[ig++] += ((e[ie] + e[ie + 2]) * cosp + (e[ie + 3] - e[ie + 1]) * sinp);
1712 g[ig++] += ((e[ie + 3] + e[ie + 1]) * cosp + (e[ie] - e[ie + 2]) * sinp);
1713 ie += 4;
1714 }
1715 }
1716 }
1717 }
1718
bfield(int rgen,int nmxi,int nmxe,int nmni,int nmne,int mmxi,int mmxe,int mmni,int mmne,int grad,int ctyp,int dtyp,int ityp,int etyp,double ep,double re,double rp,double rm,double tm,double clat,double elon,double h,double dst,double dstt,double * rse,int * nc,int * na,double * ro,double * theta,int * atyp,int * dsti,int * bori,int * bkni,double * bkpi,int * tdgi,int * dste,int * bore,int * bkne,double * bkpe,int * tdge,double * a,double * b,double * c,double * p,double * r,double * t,int * u,double * w,double * dsdc,double * dldc,double * dlda,int * cerr)1719 static void bfield(int rgen, int nmxi, int nmxe, int nmni, int nmne, int mmxi, int mmxe, int mmni,
1720 int mmne, int grad, int ctyp, int dtyp, int ityp, int etyp, double ep, double re,
1721 double rp, double rm, double tm, double clat, double elon, double h, double dst, double dstt,
1722 double *rse, int *nc, int *na, double *ro, double *theta, int *atyp, int *dsti, int *bori, int *bkni,
1723 double *bkpi, int *tdgi, int *dste, int *bore, int *bkne, double *bkpe, int *tdge, double *a,
1724 double *b, double *c, double *p, double *r, double *t, int *u, double *w, double *dsdc,
1725 double *dldc, double *dlda, int *cerr) {
1726
1727 int ia, ie, ii, np, ns, nce, nci, nse, nsi, nmax, mmin, mmax;
1728 double phi;
1729
1730 if (*cerr <= 49) {
1731 phi = elon;
1732 prebf_(&rgen, &ityp, &etyp, &dtyp, &grad, &nmni, &nmxi, &nmne, &nmxe, &mmni, &mmxi, &mmne, &mmxe, &nmax, &mmin, &mmax, &ns, &nsi,
1733 &nse, nc, &nci, &nce, na, &np, &ii, &ie, &atyp[0], &dsti[0], &bori[0], &bkni[0], &tdgi[0], &dste[0], &bore[0],
1734 &bkne[0], &tdge[0], &u[0], cerr);
1735 if (*cerr >= 50) return;
1736 if (*nc > 0) {
1737 fdsdc_(&rgen, &ityp, &etyp, &nsi, &nse, nc, &nci, &ep, &tm, &dst, &dstt, &dsti[0], &bori[0],
1738 &bkni[0], &bkpi[0], &tdgi[0], &dste[0], &bore[0], &bkne[0], &bkpe[0], &tdge[0], &u[0],
1739 &w[0], &dsdc[0], cerr);
1740 if (*cerr >= 50) return;
1741 fdlds_(&rgen, &grad, &ctyp, &clat, &phi, &h, &re, &rp, &rm, ro, &nsi, nc, &nci, &np, &ii, &ie,
1742 &nmni, &nmxi, &nmne, &nmxe, &nmax, &mmni, &mmxi, &mmne, &mmxe, &mmin, &mmax, theta, &p[0],
1743 &r[0], &t[0], &u[0], &w[0], &dldc[0], cerr);
1744 if (*cerr >= 50) return;
1745 if (rgen > 0) {
1746 rgen = 0;
1747 memset(b, 0, (grad * 36 + 28) * sizeof(double));
1748 fdldc(grad, *nc, &dsdc[0], &dldc[0]);
1749 blgen(grad, *nc, &b[0], &c[0], &dldc[0]);
1750 bngen_(&b[0]);
1751 }
1752 if (dtyp == 2) {
1753 tec(grad, atyp[0], *nc, theta, &phi, &b[0], &dldc[0], &w[0]);
1754 tse(grad, atyp[1], *nc, &rse[0], &b[0], &dldc[0], &w[0]);
1755 }
1756 }
1757 r8vgathp(1, 1, 15, 14, &b[0], &b[0]);
1758 if (grad == 1)
1759 r8vgathp(29, 1, 47, 18, &b[0], &b[0]);
1760
1761 ia = 0;
1762 if (*na > 0) {
1763 memset(dlda, 0, (*na * 6) * sizeof(double));
1764 if (dtyp == 1)
1765 tbi_(&atyp[0], na, &ia, &a[0], &b[0], &dlda[0]);
1766 else if (dtyp == 2) {
1767 tms(grad, atyp[2], *nc, *na, ia, &a[0], &b[0], &dldc[0], &dlda[0], &w[0]);
1768 tnm_(&grad, &atyp[3], nc, na, &ia, &a[0], &b[0], &dldc[0], &dlda[0], &w[0]);
1769 tvn_(&grad, &atyp[4], nc, na, &ia, &a[0], &b[0], &dldc[0], &dlda[0], &w[0]);
1770 tbi_(&atyp[5], na, &ia, &a[0], &b[0], &dlda[0]);
1771 }
1772 }
1773 }
1774 }
1775
prebf_(int * rgen,int * ityp,int * etyp,int * dtyp,int * grad,int * nmni,int * nmxi,int * nmne,int * nmxe,int * mmni,int * mmxi,int * mmne,int * mmxe,int * nmax,int * mmin,int * mmax,int * ns,int * nsi,int * nse,int * nc,int * nci,int * nce,int * na,int * np,int * ii,int * ie,int * atyp,int * dsti,int * bori,int * bkni,int * tdgi,int * dste,int * bore,int * bkne,int * tdge,int * u,int * cerr)1776 void prebf_(int *rgen, int *ityp, int *etyp, int *dtyp, int *grad, int *nmni, int *nmxi, int *
1777 nmne, int *nmxe, int *mmni, int *mmxi, int *mmne, int *mmxe, int *nmax, int *mmin, int *mmax, int *
1778 ns, int *nsi, int *nse, int *nc, int *nci, int *nce, int *na, int *np, int *ii, int *ie, int *
1779 atyp, int *dsti, int *bori, int *bkni, int *tdgi, int *dste, int *bore, int *bkne, int *tdge, int *u, int *cerr) {
1780
1781 static int nx = 0;
1782 int i__1, edst, esvr, idst, isvr;
1783 gmt_M_unused(grad);
1784
1785 /* Parameter adjustments */
1786 --u;
1787 --tdge;
1788 --bkne;
1789 --bore;
1790 --dste;
1791 --tdgi;
1792 --bkni;
1793 --bori;
1794 --dsti;
1795 --atyp;
1796
1797 /* Function Body */
1798 if (*rgen == 1) {
1799 i__1 = MIN(MIN(*nmni,*nmxi), *nmne);
1800 if (MIN(i__1,*nmxe) < 0) {
1801 fprintf(stderr, "SUBROUTINE BFIELD -- ERROR CODE 50 -- NMNI, NMXI, NMNE, OR NMXE < 0 -- ABORT\n");
1802 *cerr = 50;
1803 return;
1804 }
1805 i__1 = MIN(MIN(*mmni,*mmxi), *mmne);
1806 if (MIN(i__1,*mmxe) < 0) {
1807 fprintf(stderr, "SUBROUTINE BFIELD -- ERROR CODE 51 -- MMNI, MMXI, MMNE, OR MMXE < 0 -- ABORT\n");
1808 *cerr = 51;
1809 return;
1810 }
1811 if (*mmni > *mmxi || *mmne > *mmxe) {
1812 fprintf(stderr, "SUBROUTINE BFIELD -- ERROR CODE 52 -- EITHER MMNI > MMXI OR MMNE > MMXE -- ABORT\n");
1813 *cerr = 52;
1814 return;
1815 }
1816 if (*mmxi > *nmxi || *mmxe > *nmxe) {
1817 fprintf(stderr, "SUBROUTINE BFIELD -- ERROR CODE 53 -- EITHER MMXI > NMXI OR MMXE > NMXE -- ABORT\n");
1818 *cerr = 53;
1819 return;
1820 }
1821 isvr = *ityp % 3;
1822 idst = *ityp / 3;
1823 esvr = *etyp % 3;
1824 edst = *etyp / 3;
1825 *nmax = MAX(*nmxi,*nmxe);
1826 *mmin = MIN(*mmni,*mmne);
1827 *mmax = MAX(*mmxi,*mmxe);
1828 *nsi = nshx(*nmxi, *nmni, *mmxi, *mmni);
1829 *nse = nshx(*nmxe, *nmne, *mmxe, *mmne);
1830 *ns = *nsi + *nse;
1831 *np = nlpx(*nmax, *mmax, *mmin);
1832 *ii = nlpx(*nmni - 1, *mmax, *mmin);
1833 *ie = nlpx(*nmne - 1, *mmax, *mmin);
1834 *nci = 0;
1835 if (*nsi > 0) {
1836 i8vset(1, *nsi, 1, &u[1]);
1837 if (isvr == 1)
1838 i8vadd(1, 1, 1, *nsi, &tdgi[1], &u[1], &u[1]);
1839 else if (isvr == 2) {
1840 i8vadd(1, 1, 1, *nsi, &bori[1], &u[1], &u[1]);
1841 i8vadd(1, 1, 1, *nsi, &bkni[1], &u[1], &u[1]);
1842 }
1843 if (idst == 1)
1844 i8vadd(1, 1, 1, *nsi, &dsti[1], &u[1], &u[1]);
1845
1846 *nci = i8ssum(1, *nsi, &u[1]);
1847 }
1848 *nce = 0;
1849 if (*nse > 0) {
1850 i__1 = *nsi + 1;
1851 i8vset(i__1, *nse, 1, &u[1]);
1852 if (esvr == 1)
1853 i8vadd(1, i__1, i__1, *nse, &tdge[1], &u[1], &u[1]);
1854 else if (esvr == 2) {
1855 i8vadd(1, i__1, i__1, *nse, &bore[1], &u[1], &u[1]);
1856 i8vadd(1, i__1, i__1, *nse, &bkne[1], &u[1], &u[1]);
1857 }
1858 if (edst == 1)
1859 i8vadd(1, i__1, i__1, *nse, &dste[1], &u[1], &u[1]);
1860
1861 *nce = i8ssum(i__1, *nse, &u[1]);
1862 }
1863 *nc = *nci + *nce;
1864 *rgen = 7;
1865 }
1866 *rgen += nx;
1867 *na = 0;
1868 nx = 0;
1869 if (*dtyp == 1)
1870 *na += MIN(1,atyp[1]) * 3;
1871 else if (*dtyp == 2) {
1872 nx += atyp[1];
1873 nx += atyp[2];
1874 *na += MIN(1,atyp[3]) * 3;
1875 *na += MIN(1,atyp[4]) * 3;
1876 *na += MIN(1,atyp[5]) * 3;
1877 nx += *na;
1878 *na += MIN(1,atyp[6]) * 3;
1879 }
1880 nx = MIN(1,nx);
1881 }
1882
fdlds_(int * rgen,int * grad,int * ctyp,double * clat,double * phi,double * h__,double * re,double * rp,double * rm,double * ro,int * nsi,int * nc,int * nci,int * np,int * ii,int * ie,int * nmni,int * nmxi,int * nmne,int * nmxe,int * nmax,int * mmni,int * mmxi,int * mmne,int * mmxe,int * mmin,int * mmax,double * theta,double * p,double * r__,double * t,int * u,double * w,double * dldc,int * cerr)1883 void fdlds_(int *rgen, int *grad, int *ctyp, double *clat, double *phi, double *h__, double *re,
1884 double *rp, double *rm, double *ro, int *nsi, int *nc, int *nci, int *np, int *ii, int *ie, int *
1885 nmni, int *nmxi, int *nmne, int *nmxe, int *nmax, int *mmni, int *mmxi, int *mmne, int *mmxe, int *
1886 mmin, int *mmax, double *theta, double *p, double *r__, double *t, int *u, double *w, double *dldc, int *cerr) {
1887 /* Initialized data */
1888
1889 static double roo = 0.;
1890 static double phio = 0.;
1891 static double thetao = 0.;
1892 static double clato = 0.;
1893
1894 /* Local variables */
1895 int m, n, ic, id, ip, lend, pgen, tgen;
1896 double fa, fc, fd, ar, fm, ra, fp, fr, fs, fn, fnm1, fnp1, fnp2, fnfp, fprr, fdrr;
1897 double pbppp = 0, pbppr = 0, pbrpp = 0, pbtpp = 0, pbppt = 0, pbtpr = 0, pbrpt = 0, pbtpt = 0, pbrpr = 0, cscth2;
1898 double cscthe, costhe, cotthe, sinthe, fmfpst;
1899 gmt_M_unused(nsi); gmt_M_unused(nci); gmt_M_unused(cerr);
1900
1901 /* Parameter adjustments */
1902 --dldc;
1903 --w;
1904 --u;
1905 --t;
1906 --r__;
1907 --p;
1908
1909 /* Function Body */
1910 geocen(*ctyp, *re, *rp, *rm, *h__, *clat, ro, theta, &sinthe, &costhe);
1911 pgen = *rgen - 6;
1912 tgen = pgen;
1913 if (phio != *phi) {
1914 pgen = 1;
1915 ++(*rgen);
1916 phio = *phi;
1917 }
1918 if (thetao != *theta) {
1919 tgen = 1;
1920 ++(*rgen);
1921 thetao = *theta;
1922 }
1923 if (clato != *clat) {
1924 ++(*rgen);
1925 clato = *clat;
1926 }
1927 if (roo != *ro) {
1928 ++(*rgen);
1929 roo = *ro;
1930 }
1931 if (*rgen > 0) {
1932 if (sinthe == 0.) {
1933 if (*grad == 0)
1934 fprintf(stderr, "SUBROUTINE BFIELD -- ERROR CODE 1 -- GEOGRAPHIC POLAR POSITION DETECTED, B-PHI INDETERMINABLE -- WARNING\n");
1935 else
1936 fprintf(stderr, "SUBROUTINE BFIELD -- ERROR CODE 2 -- GEOGRAPHIC POLAR POSITION DETECTED, B-PHI AND\n\t\t\t-- PHI-DERIVATIVE GRADIENT COMPONENTS INDETERMINABLE -- WARNING\n");
1937
1938 *phi = 0.;
1939 cscthe = 0.;
1940 cscth2 = 0.;
1941 cotthe = 0.;
1942 } else {
1943 cscthe = 1. / sinthe;
1944 cscth2 = cscthe * cscthe;
1945 cotthe = costhe * cscthe;
1946 }
1947 if (tgen > 0)
1948 schmit_(grad, rgen, nmax, mmin, mmax, &sinthe, &costhe, &p[1], &r__[1]);
1949
1950 if (pgen > 0)
1951 trigmp(*mmax, *phi, &t[1]);
1952
1953 memset(&dldc[1], 0, ((*grad * 18 + 6) * *nc) * sizeof(double));
1954 ic = 0;
1955 id = 1;
1956 ip = *ii;
1957 ar = *rm / *ro;
1958 fa = pow_di(ar, *nmni + 1);
1959 for (n = *nmni; n <= *nmxi; ++n) {
1960 fnp1 = (double) (n + 1);
1961 fnp2 = (double) (n + 2);
1962 fa *= ar;
1963 for (m = *mmin; m <= MIN(n,*mmax); ++m) {
1964 ++ip;
1965 if (m >= *mmni && m <= *mmxi) {
1966 ++ic;
1967 fm = (double) m;
1968 fp = fa * p[ip];
1969 fd = -fa * p[*np + ip];
1970 fc = t[m + 1];
1971 fs = t[*mmax + m + 2];
1972 fnfp = fnp1 * fp;
1973 fmfpst = fm * fp * cscthe;
1974 lend = u[ic];
1975 r8vset(id, u[ic], fd * fc, &dldc[1]);
1976 r8vset(*nc + id, u[ic], fmfpst * fs, &dldc[1]);
1977 r8vset((*nc << 1) + id, u[ic], fnfp * fc, &dldc[1]);
1978 if (*grad == 1) {
1979 fprr = fp / *ro;
1980 fdrr = fd / *ro;
1981 pbtpt = -fa * p[(*np << 1) + ip] / *ro;
1982 pbtpp = -fm * fdrr * cscthe;
1983 pbtpr = -fnp2 * fdrr;
1984 pbppt = -fm * (fdrr + fprr * cotthe) * cscthe;
1985 pbppp = fm * fm * fprr * cscth2;
1986 pbppr = -fnp2 * fm * fprr * cscthe;
1987 pbrpt = -fnp1 * fdrr;
1988 pbrpp = -fnp1 * fm * fprr * cscthe;
1989 pbrpr = -fnp1 * fnp2 * fprr;
1990 r8vset(*nc * 6 + id, lend, pbtpt * fc, &dldc[1]);
1991 r8vset(*nc * 7 + id, lend, pbtpp * fs, &dldc[1]);
1992 r8vset((*nc << 3) + id, lend, pbtpr * fc, &dldc[1]);
1993 r8vset(*nc * 9 + id, lend, pbppt * fs, &dldc[1]);
1994 r8vset(*nc * 10 + id, lend, pbppp * fc, &dldc[1]);
1995 r8vset(*nc * 11 + id, lend, pbppr * fs, &dldc[1]);
1996 r8vset(*nc * 12 + id, lend, pbrpt * fc, &dldc[1]);
1997 r8vset(*nc * 13 + id, lend, pbrpp * fs, &dldc[1]);
1998 r8vset(*nc * 14 + id, lend, pbrpr * fc, &dldc[1]);
1999 }
2000 id += lend;
2001 if (m > 0) {
2002 ++ic;
2003 lend = u[ic];
2004 r8vset(id, lend, fd * fs, &dldc[1]);
2005 r8vset(*nc + id, lend, -fmfpst * fc, &dldc[1]);
2006 r8vset((*nc << 1) + id, lend, fnfp * fs, &dldc[1]);
2007 if (*grad == 1) {
2008 r8vset(*nc * 6 + id, lend, pbtpt * fs, &dldc[1]);
2009 r8vset(*nc * 7 + id, lend, -pbtpp * fc, &dldc[1]);
2010 r8vset((*nc << 3) + id, lend, pbtpr * fs, &dldc[1]);
2011 r8vset(*nc * 9 + id, lend, -pbppt * fc, &dldc[1]);
2012 r8vset(*nc * 10 + id, lend, pbppp * fs, &dldc[1]);
2013 r8vset(*nc * 11 + id, lend, -pbppr * fc, &dldc[1]);
2014 r8vset(*nc * 12 + id, lend, pbrpt * fs, &dldc[1]);
2015 r8vset(*nc * 13 + id, lend, -pbrpp * fc, &dldc[1]);
2016 r8vset(*nc * 14 + id, lend, pbrpr * fs, &dldc[1]);
2017 }
2018 id += lend;
2019 }
2020 }
2021 }
2022 }
2023 ip = *ie;
2024 ra = *ro / *rm;
2025 fr = pow_di(ra, *nmne - 2);
2026 for (n = *nmne; n <= *nmxe; ++n) {
2027 fnm1 = (double) (n - 1);
2028 fn = (double) n;
2029 fr *= ra;
2030 for (m = *mmin; m <= MIN(n,*mmax); ++m) {
2031 ++ip;
2032 if (m >= *mmne && m <= *mmxe) {
2033 ++ic;
2034 fm = (double) m;
2035 fp = fr * p[ip];
2036 fd = -fr * p[*np + ip];
2037 fc = t[m + 1];
2038 fs = t[*mmax + m + 2];
2039 fnfp = -fn * fp;
2040 fmfpst = fm * fp * cscthe;
2041 lend = u[ic];
2042 r8vset(id, lend, fd * fc, &dldc[1]);
2043 r8vset(*nc + id, lend, fmfpst * fs, &dldc[1]);
2044 r8vset((*nc << 1) + id, lend, fnfp * fc, &dldc[1]);
2045 if (*grad == 1) {
2046 fprr = fp / *ro;
2047 fdrr = fd / *ro;
2048 pbtpt = -fr * p[(*np << 1) + ip] / *ro;
2049 pbtpp = -fm * fdrr * cscthe;
2050 pbtpr = fnm1 * fdrr;
2051 pbppt = -fm * (fdrr + fprr * cotthe) * cscthe;
2052 pbppp = fm * fm * fprr * cscth2;
2053 pbppr = fnm1 * fm * fprr * cscthe;
2054 pbrpt = fn * fdrr;
2055 pbrpp = fn * fm * fprr * cscthe;
2056 pbrpr = -fn * fnm1 * fprr;
2057 r8vset(*nc * 6 + id, lend, pbtpt * fc, &dldc[1]);
2058 r8vset(*nc * 7 + id, lend, pbtpp * fs, &dldc[1]);
2059 r8vset((*nc << 3) + id, lend, pbtpr * fc, &dldc[1]);
2060 r8vset(*nc * 9 + id, lend, pbppt * fs, &dldc[1]);
2061 r8vset(*nc * 10 + id, lend, pbppp * fc, &dldc[1]);
2062 r8vset(*nc * 11 + id, lend, pbppr * fs, &dldc[1]);
2063 r8vset(*nc * 12 + id, lend, pbrpt * fc, &dldc[1]);
2064 r8vset(*nc * 13 + id, lend, pbrpp * fs, &dldc[1]);
2065 r8vset(*nc * 14 + id, lend, pbrpr * fc, &dldc[1]);
2066 }
2067 id += lend;
2068 if (m > 0) {
2069 ++ic;
2070 lend = u[ic];
2071 r8vset(id, lend, fd * fs, &dldc[1]);
2072 r8vset(*nc + id, lend, -fmfpst * fc, &dldc[1]);
2073 r8vset((*nc << 1) + id, lend, fnfp * fs, &dldc[1]);
2074 if (*grad == 1) {
2075 r8vset(*nc * 6 + id, lend, pbtpt * fs, &dldc[1]);
2076 r8vset(*nc * 7 + id, lend, -pbtpp * fc, &dldc[1]);
2077 r8vset((*nc << 3) + id, lend, pbtpr * fs, &dldc[1]);
2078 r8vset(*nc * 9 + id, lend, -pbppt * fc, &dldc[1]);
2079 r8vset(*nc * 10 + id, lend, pbppp * fs, &dldc[1]);
2080 r8vset(*nc * 11 + id, lend, -pbppr * fc, &dldc[1]);
2081 r8vset(*nc * 12 + id, lend, pbrpt * fs, &dldc[1]);
2082 r8vset(*nc * 13 + id, lend, -pbrpp * fc, &dldc[1]);
2083 r8vset(*nc * 14 + id, lend, pbrpr * fs, &dldc[1]);
2084 }
2085 id += lend;
2086 }
2087 }
2088 }
2089 }
2090 tdc(*grad, *nc, *clat, *theta, &dldc[1], &w[1]);
2091 }
2092 }
2093
geocen(int ctyp,double re,double rp,double rm,double h,double clat,double * r,double * theta,double * sinthe,double * costhe)2094 void geocen(int ctyp, double re, double rp, double rm, double h, double clat, double *r, double *theta, double *sinthe, double *costhe) {
2095 int ifirst = 1;
2096 double re2 = 0, rp2 = 0, kappa, rhoew, rhons, costh2, sinth2;
2097
2098 *theta = clat;
2099 *r = rm + h;
2100 /* *sinthe = sin(*theta);
2101 *costhe = cos(*theta);*/
2102 sincos(*theta, sinthe, costhe);
2103 if (ctyp == 0) {
2104 if (ifirst == 1) {
2105 re2 = re * re;
2106 rp2 = rp * rp;
2107 }
2108 sinth2 = *sinthe * *sinthe;
2109 costh2 = *costhe * *costhe;
2110 kappa = 1. / sqrt(re2 * sinth2 + rp2 * costh2);
2111 rhoew = re2 * kappa + h;
2112 rhons = rp2 * kappa + h;
2113 *theta = atan2(rhoew * *sinthe, rhons * *costhe);
2114 *r = sqrt(rhoew * rhoew * sinth2 + rhons * rhons * costh2);
2115 /* *sinthe = sin(*theta);
2116 *costhe = cos(*theta); */
2117 sincos(*theta, sinthe, costhe);
2118 }
2119 }
2120
schmit_(int * grad,int * rgen,int * nmax,int * mmin,int * mmax,double * sinthe,double * costhe,double * p,double * r)2121 void schmit_(int *grad, int *rgen, int *nmax, int *mmin, int *mmax, double *sinthe, double *costhe, double *p, double *r) {
2122
2123 int l, n, np = 0, kd0, kd1, kd2, kp0, kp1, kp2, kq2, kq1, kq0, ksm2 = 0, ktm2 = 0, ksec, ktes, knnp, nd0sec = 0;
2124 int np1sec = 0, nd0tes = 0, nd1tes = 0, np1tes = 0, np2tes = 0, nq0nnp = 0, nq0msq = 0;
2125 double cscth2, costh2, sinth2, cscthe, cotthe, cthsth;
2126 double root3 = 1.732050807568877;
2127
2128 /* Parameter adjustments */
2129 --r;
2130 --p;
2131
2132 /* Function Body */
2133 if (*rgen > 6)
2134 srecur_(grad, nmax, mmin, mmax, &ksm2, &ktm2, &np, &np1sec, &nd0sec, &np1tes, &np2tes, &nd0tes, &nd1tes, &nq0nnp, &nq0msq, &r[1]);
2135
2136 if (*nmax >= 1) {
2137 kp0 = 1;
2138 kp2 = kp0;
2139 if (*mmin <= 0) {
2140 p[kp0] = *costhe;
2141 p[np + kp0] = -(*sinthe);
2142 if (*grad == 1)
2143 p[(np << 1) + kp0] = -(*costhe);
2144
2145 ++kp0;
2146 }
2147 if (*mmax >= 1) {
2148 p[kp0] = *sinthe;
2149 p[np + kp0] = *costhe;
2150 if (*grad == 1)
2151 p[(np << 1) + kp0] = -(*sinthe);
2152
2153 ++kp0;
2154 }
2155 if (*nmax >= 2) {
2156 kp1 = kp0;
2157 sinth2 = *sinthe * *sinthe;
2158 costh2 = *costhe * *costhe;
2159 cthsth = *costhe * *sinthe;
2160 if (*mmin <= 0) {
2161 p[kp0] = (costh2 * 3. - 1.) / 2.;
2162 p[np + kp0] = cthsth * -3.;
2163 if (*grad == 1) {
2164 p[(np << 1) + kp0] = (costh2 * 2. - 1.) * -3.;
2165 }
2166 ++kp0;
2167 }
2168 if (*mmin <= 1 && *mmax >= 1) {
2169 p[kp0] = root3 * cthsth;
2170 p[np + kp0] = root3 * (costh2 * 2. - 1.);
2171 if (*grad == 1) {
2172 p[(np << 1) + kp0] = root3 * -4. * cthsth;
2173 }
2174 ++kp0;
2175 }
2176 if (*mmax >= 2) {
2177 p[kp0] = root3 * sinth2 / 2.;
2178 p[np + kp0] = root3 * cthsth;
2179 if (*grad == 1) {
2180 p[(np << 1) + kp0] = root3 * (costh2 * 2. - 1.);
2181 }
2182 ++kp0;
2183 }
2184 kd2 = np + kp2;
2185 kd1 = np + kp1;
2186 kd0 = np + kp0;
2187 kq2 = np + kd2;
2188 kq1 = np + kd1;
2189 kq0 = np + kd0;
2190 ksec = 1;
2191 ktes = 1;
2192 knnp = 1;
2193 if (*sinthe == 0.) {
2194 for (n = 3; n <= *nmax; ++n) {
2195 l = MAX(0, MIN(n-1,*mmax) - *mmin + 1);
2196 if (l > 0) {
2197 r8vset(kp0, l, 0., &p[1]);
2198 r8vlinkq(kp1, np1tes + ktes, kp0, l, *costhe, &p[1], &r[1], &p[1]);
2199 r8vlinkq(kp2, np2tes + ktes, kp0, l, -1., &p[1], &r[1], &p[1]);
2200 r8vset(kd0, l, 0., &p[1]);
2201 r8vlinkq(kd1, np1tes + ktes, kd0, l, *costhe, &p[1], &r[1], &p[1]);
2202 r8vlinkq(kd2, np2tes + ktes, kd0, l, -1., &p[1], &r[1], &p[1]);
2203 if (*grad == 1) {
2204 r8vset(kq0, l, 0., &p[1]);
2205 r8vlinkq(kq1, np1tes + ktes, kq0, l, *costhe, &p[1], &r[1], &p[1]);
2206 r8vlinkq(kp1, np1tes + ktes, kq0, l, -(*costhe), &p[1], &r[1], &p[1]);
2207 r8vlinkq(kq2, np2tes + ktes, kq0, l, -1., &p[1], &r[1], &p[1]);
2208 }
2209 ktes += l;
2210 }
2211 if (n <= *mmax) {
2212 p[kp0 + l] = 0.;
2213 p[kd0 + l] = 0.;
2214 if (*grad == 1)
2215 p[kq0 + l] = 0.;
2216 ++l;
2217 }
2218 kp2 = kp1;
2219 kp1 = kp0;
2220 kp0 += l;
2221 kd2 = kd1;
2222 kd1 = kd0;
2223 kd0 += l;
2224 kq2 = kq1;
2225 kq1 = kq0;
2226 kq0 += l;
2227 }
2228 } else {
2229 cotthe = *costhe / *sinthe;
2230 cscthe = 1. / *sinthe;
2231 cscth2 = cscthe * cscthe;
2232 for (n = 3; n <= *nmax; ++n) {
2233 l = MAX(0, MIN(n-1,*mmax) - *mmin + 1);
2234 if (l > 0) {
2235 r8vset(kp0, l, 0., &p[1]);
2236 r8vlinkq(kp1, np1tes + ktes, kp0, l, *costhe, &p[1], &r[1], &p[1]);
2237 r8vlinkq(kp2, np2tes + ktes, kp0, l, -1., &p[1], &r[1], &p[1]);
2238 r8vset(kd0, l, 0., &p[1]);
2239 r8vlinkq(kp0, nd0tes + ktes, kd0, l, cotthe, &p[1], &r[1], &p[1]);
2240 r8vlinkq(kp1, nd1tes + ktes, kd0, l, -cscthe, &p[1], &r[1], &p[1]);
2241 if (*grad == 1) {
2242 r8vset(kq0, l, 0., &p[1]);
2243 r8vlinkq(kp0, nq0msq + ktm2, kq0, l, cscth2, &p[1], &r[1], &p[1]);
2244 r8vlinkt(kp0, kq0, l, -r[nq0nnp + knnp], &p[1], &p[1]);
2245 r8vlinkt(kd0, kq0, l, -cotthe, &p[1], &p[1]);
2246 }
2247 ktes += l;
2248 }
2249 if (n <= *mmax) {
2250 p[kp0 + l] = r[np1sec + ksec] * *sinthe * p[kp0 - 1];
2251 p[kd0 + l] = r[nd0sec + ksec] * cotthe * p[kp0 + l];
2252 if (*grad == 1)
2253 p[kq0 + l] = (r[nq0msq - ksm2 + n + 1] * cscth2 -r[nq0nnp + knnp]) * p[kp0 + l] - cotthe * p[kd0 + l];
2254
2255 ++ksec;
2256 ++l;
2257 }
2258 ++knnp;
2259 kp2 = kp1;
2260 kp1 = kp0;
2261 kp0 += l;
2262 kd0 += l;
2263 kq0 += l;
2264 }
2265 }
2266 }
2267 }
2268 }
2269
srecur_(int * grad,int * nmax,int * mmin,int * mmax,int * ksm2,int * ktm2,int * npall,int * nad1,int * nad2,int * nad3,int * nad4,int * nad5,int * nad6,int * nad7,int * nad8,double * r)2270 void srecur_(int *grad, int *nmax, int *mmin, int *mmax, int *ksm2, int *ktm2, int *npall, int *
2271 nad1, int *nad2, int *nad3, int *nad4, int *nad5, int *nad6, int *nad7, int *nad8, double *r) {
2272 int m, n, ksec, kmin, lmax, kmax, ktes, knnp, kmsq, nnnp1, nsect, ntess;
2273 double fn, f2n, fnl1, fnp1, f2nl1, fnsq, fmsq, fnl1sq, fdsqrt;
2274
2275 /* Parameter adjustments */
2276 --r;
2277
2278 /* Function Body */
2279 lmax = MIN(*nmax,2);
2280 kmax = MIN(*mmax,2);
2281 kmin = MIN(*mmin,2);
2282 *ksm2 = MIN(*mmin,3);
2283 *ktm2 = I_DIM(*mmin, 3) + 1;
2284 *npall = nlpx(*nmax, *mmax, *mmin);
2285 ntess = *npall - nlpx(lmax, kmax, kmin) + kmax - *mmax;
2286 nsect = MAX(0, *mmax - 2);
2287 nnnp1 = MAX(0, *nmax - 2);
2288 *nad1 = 0;
2289 *nad2 = nsect;
2290 *nad3 = nsect + *nad2;
2291 *nad4 = ntess + *nad3;
2292 *nad5 = ntess + *nad4;
2293 *nad6 = ntess + *nad5;
2294 *nad7 = ntess + *nad6;
2295 *nad8 = nnnp1 + *nad7;
2296 ksec = ktes = knnp = 0;
2297 for (n = 3; n <= *nmax; ++n) {
2298 fnp1 = (double) (n + 1);
2299 fn = (double) n;
2300 fnl1 = (double) (n - 1);
2301 f2n = fn * 2.;
2302 f2nl1 = f2n - 1.;
2303 fnsq = fn * fn;
2304 fnl1sq = fnl1 * fnl1;
2305 if (n <= *mmax) {
2306 ++ksec;
2307 r[*nad1 + ksec] = sqrt(f2nl1 / f2n);
2308 r[*nad2 + ksec] = fn;
2309 }
2310 if (*grad == 1) {
2311 ++knnp;
2312 r[*nad7 + knnp] = fn * fnp1;
2313 }
2314 for (m = *mmin; m <= MIN(n-1,*mmax); ++m) {
2315 ++ktes;
2316 fmsq = (double)(m * m);
2317 fdsqrt = sqrt(fnsq - fmsq);
2318 r[*nad3 + ktes] = f2nl1 / fdsqrt;
2319 r[*nad4 + ktes] = sqrt(fnl1sq - fmsq) / fdsqrt;
2320 r[*nad5 + ktes] = fn;
2321 r[*nad6 + ktes] = fdsqrt;
2322 }
2323 }
2324 if (*grad == 1) {
2325 kmsq = 0;
2326 for (m = *ksm2; m <= *mmax; ++m) {
2327 ++kmsq;
2328 r[*nad8 + kmsq] = (double)(m * m);
2329 }
2330 }
2331 }
2332
trigmp(int mmax,double phi,double * t)2333 void trigmp(int mmax, double phi, double *t) {
2334 int m;
2335
2336 /* Parameter adjustments */
2337 --t;
2338
2339 t[1] = 1.;
2340 t[mmax + 2] = 0.;
2341 if (mmax > 0) {
2342 t[2] = cos(phi);
2343 t[mmax + 3] = sin(phi);
2344 for (m = 2; m <= mmax; ++m) {
2345 t[m + 1] = t[2] * 2. * t[m] - t[m - 1];
2346 t[mmax + m + 2] = t[2] * 2. * t[mmax + m + 1] - t[mmax + m];
2347 }
2348 }
2349 }
2350
tdc(int grad,int nc,double clat,double theta,double * dldc,double * r)2351 void tdc(int grad, int nc, double clat, double theta, double *dldc, double *r) {
2352 double cpsi, spsi;
2353
2354 sincos(theta - clat, &spsi, &cpsi);
2355 r[0] = -cpsi;
2356 r[1] = 0.;
2357 r[2] = -spsi;
2358 r[3] = 0.;
2359 r[4] = 1.;
2360 r[5] = 0.;
2361 r[6] = spsi;
2362 r[7] = 0.;
2363 r[8] = -cpsi;
2364 ltranv(1, nc, nc, &r[0], &dldc[0]);
2365 if (grad == 1) {
2366 ltranv(0, nc * 3, nc * 3, &r[0], &dldc[nc * 6]);
2367 ltranv(0, nc, nc, &r[0], &dldc[nc * 6]);
2368 ltranv(0, nc, nc, &r[0], &dldc[nc * 9]);
2369 ltranv(0, nc, nc, &r[0], &dldc[nc * 12]);
2370 }
2371 }
2372
fdsdc_(int * rgen,int * ityp,int * etyp,int * nsi,int * nse,int * nc,int * nci,double * ta,double * tb,double * dst,double * dstt,int * dsti,int * bori,int * bkni,double * bkpi,int * tdgi,int * dste,int * bore,int * bkne,double * bkpe,int * tdge,int * u,double * w,double * dsdc,int * cerr)2373 void fdsdc_(int *rgen, int *ityp, int *etyp, int *nsi, int *nse, int *nc, int *nci, double *ta,
2374 double *tb, double *dst, double *dstt, int *dsti, int *bori, int *bkni, double *bkpi, int *tdgi,
2375 int *dste, int *bore, int *bkne, double *bkpe, int *tdge, int *u, double *w, double *dsdc, int *cerr) {
2376
2377 static double tbo = 0.;
2378
2379 /* Local variables */
2380 int i__1, tgen, edst, idst, esvr, isvr;
2381
2382 /* Parameter adjustments */
2383 --dsdc;
2384 --w;
2385 --u;
2386 --tdge;
2387 --bkpe;
2388 --bkne;
2389 --bore;
2390 --dste;
2391 --tdgi;
2392 --bkpi;
2393 --bkni;
2394 --bori;
2395 --dsti;
2396
2397 /* Function Body */
2398 tgen = MAX(0,*rgen - 6);
2399 if (tbo != *tb) {
2400 tgen = MIN(1,tgen + *ityp + *etyp);
2401 *rgen += tgen;
2402 tbo = *tb;
2403 }
2404 if (tgen > 0) {
2405 r8vset(1, *nc << 1, 0., &dsdc[1]);
2406 if (*nsi > 0) {
2407 isvr = *ityp % 3;
2408 idst = *ityp / 3;
2409 i8vcum(1, 1, *nsi, &u[1]);
2410 r8vscats(1, *nsi, 1., &u[1], &dsdc[1]);
2411 r8vscats(1, *nsi, 0., &u[1], &dsdc[*nc + 1]);
2412 i8vadds(1, 1, *nsi, 1, &u[1], &u[1]);
2413 if (isvr == 1)
2414 taylor(*nc, *nsi, *ta, *tb, &tdgi[1], &u[1], &w[1], &dsdc[1]);
2415 else if (isvr == 2) {
2416 bsplyn(*nc, *nsi, ta, tb, &bori[1], &bkni[1], &bkpi[1], &u[1], &w[1], &dsdc[1], cerr);
2417 if (*cerr >= 50) goto L10;
2418 }
2419 if (idst == 1)
2420 dstorm(*nc, *nsi, dst, dstt, &dsti[1], &u[1], &dsdc[1]);
2421 L10:
2422 i8vdel(1, 1, *nsi, &u[1]);
2423 }
2424 if (*nse > 0) {
2425 esvr = *etyp % 3;
2426 edst = *etyp / 3;
2427 i__1 = *nsi + 1;
2428 i8vcum(1, i__1, *nse, &u[1]);
2429 r8vscats(i__1, *nse, 1., &u[1], &dsdc[*nci + 1]);
2430 r8vscats(i__1, *nse, 0., &u[1], &dsdc[*nc + *nci + 1]);
2431 i8vadds(i__1, i__1, *nse, 1, &u[1], &u[1]);
2432 if (esvr == 1)
2433 taylor(*nc, *nse, *ta, *tb, &tdge[1], &u[*nsi + 1], &w[1], &dsdc[*nci + 1]);
2434 else if (esvr == 2) {
2435 bsplyn(*nc, *nse, ta, tb, &bore[1], &bkne[1], &bkpe[1], &u[*nsi + 1], &w[1], &dsdc[*nci + 1], cerr);
2436 if (*cerr >= 50) goto L20;
2437 }
2438 if (edst == 1)
2439 dstorm(*nc, *nse, dst, dstt, &dste[1], &u[*nsi + 1], &dsdc[*nci + 1]);
2440 L20:
2441 i8vdel(1, *nsi + 1, *nse, &u[1]);
2442 }
2443 }
2444 }
2445
taylor(int nc,int ns,double ta,double tb,int * tdeg,int * u,double * dsdt,double * dsdc)2446 void taylor(int nc, int ns, double ta, double tb, int *tdeg, int *u, double *dsdt, double *dsdc) {
2447 int i, j, n, iu;
2448 double dt;
2449
2450 /* Parameter adjustments */
2451 --dsdt;
2452 --u;
2453 --tdeg;
2454
2455 /* Function Body */
2456 dt = tb - ta;
2457 for (i = 1; i <= ns; ++i) {
2458 n = tdeg[i];
2459 if (n > 0) {
2460 iu = u[i];
2461 dsdt[1] = 1.;
2462 for (j = 1; j <= n; ++j) {
2463 dsdt[j + 1] = dsdt[j] * dt / (double) j;
2464 }
2465 r8vgathp(2, 1, iu, n, &dsdt[1], &dsdc[0]);
2466 r8vgathp(1, 1, nc + iu, n, &dsdt[1], &dsdc[0]);
2467 u[i] += n;
2468 }
2469 }
2470 }
2471
bsplyn(int nc,int ns,double * ta,double * tb,int * bord,int * bkno,double * bkpo,int * u,double * dtdb,double * dsdc,int * cerr)2472 void bsplyn(int nc, int ns, double *ta, double *tb, int *bord, int *bkno, double *bkpo, int *u, double *dtdb, double *dsdc, int *cerr) {
2473 int i, k, l, m, n, ik, iu;
2474
2475 /* Parameter adjustments */
2476 --dsdc;
2477 --dtdb;
2478 --u;
2479 --bkpo;
2480 --bkno;
2481 --bord;
2482
2483 /* Function Body */
2484 ik = 1;
2485 for (i = 1; i <= ns; ++i) {
2486 n = bord[i];
2487 k = bkno[i];
2488 l = n + k;
2489 iu = u[i];
2490 if (n > 0) {
2491 m = n + 1;
2492 sbspln_(ta, tb, &m, &k, &bkpo[ik], &dtdb[1], &dsdc[iu], cerr);
2493 if (*cerr >= 50) return;
2494 r8vset(1, l + 1, 0., &dtdb[1]);
2495 tbspln_(tb, &n, &k, &bkpo[ik], &dtdb[1], cerr);
2496 if (*cerr >= 50) return;
2497 r8vgathp(1, 1, nc + iu, l, &dtdb[1], &dsdc[1]);
2498 } else if (k > 0) {
2499 r8vset(iu, l, 0., &dsdc[1]);
2500 r8vset(nc + iu, l, 0., &dsdc[1]);
2501 }
2502 u[i] += l;
2503 ik = ik + k + 2;
2504 }
2505 }
2506
sbspln_(double * ta,double * tb,int * n,int * k,double * bkpo,double * dtdb,double * dsdc,int * cerr)2507 void sbspln_(double *ta, double *tb, int *n, int *k, double *bkpo, double *dtdb, double *dsdc, int *cerr) {
2508 int i, ik, jk, ip, nl, np, ns;
2509 double rn;
2510
2511 /* Parameter adjustments */
2512 --dsdc;
2513 --dtdb;
2514 --bkpo;
2515
2516 /* Function Body */
2517 if (*n > 1) {
2518 nl = *n + *k + 1;
2519 r8vset(1, nl << 1, 0., &dtdb[1]);
2520 tbspln_(tb, n, k, &bkpo[1], &dtdb[1], cerr);
2521 if (*cerr >= 50) return;
2522 tbspln_(ta, n, k, &bkpo[1], &dtdb[nl + 1], cerr);
2523 if (*cerr >= 50) return;
2524 r8vsub(nl+1, 1, 1, nl, &dtdb[1], &dtdb[1], &dtdb[1]);
2525 rn = 1. / (double) (*n - 1);
2526 np = *n + *k - 1;
2527 ns = np;
2528 for (i = 1; i <= np; ++i) {
2529 ip = i + 1;
2530 ik = MIN(ip,*k + 2);
2531 jk = MAX(ip - *n + 1,1);
2532 dsdc[i] = (bkpo[ik] - bkpo[jk]) * r8ssum_(&ip, &ns, &dtdb[1]);
2533 --ns;
2534 }
2535 if (np > 0)
2536 r8vscale(1, np, rn, &dsdc[1]);
2537 }
2538 }
2539
tbspln_(double * t,int * n,int * k,double * bkpo,double * dtdb,int * cerr)2540 void tbspln_(double *t, int *n, int *k, double *bkpo, double *dtdb, int *cerr) {
2541 int i, j, l, m, p, y, ik, jk, posb;
2542 double difi, difj, delk, temp;
2543
2544 /* Parameter adjustments */
2545 --dtdb;
2546 --bkpo;
2547
2548 /* Function Body */
2549 if (*t >= bkpo[1] && *t <= bkpo[*k + 2]) {
2550 r8slt(1, *k + 2, *t, &bkpo[1], &y);
2551 l = MIN(*k+2,y+1);
2552 p = l + *n - 2;
2553 if (*n == 1)
2554 dtdb[p] = 1.;
2555 else if (*n > 1) {
2556 m = l;
2557 ik = MIN(m,*k+2);
2558 jk = MAX(m-1,1);
2559 difi = bkpo[ik] - *t;
2560 delk = bkpo[ik] - bkpo[jk];
2561 if (delk == 0.)
2562 dtdb[p] = 0.;
2563 else
2564 dtdb[p] = 1. / delk;
2565
2566 posb = p - 1;
2567 for (j = 2; j <= *n; ++j) {
2568 jk = MAX(m-j,1);
2569 temp = difi * dtdb[posb + 1];
2570 delk = bkpo[ik] - bkpo[jk];
2571 if (delk == 0.)
2572 temp = 0.;
2573 else if (j < *n)
2574 temp /= delk;
2575
2576 dtdb[posb] = temp;
2577 --posb;
2578 }
2579 dtdb[p + 1] = 0.;
2580 ++m;
2581 for (i = 2; i <= *n; ++i) {
2582 ik = MIN(m,*k+2);
2583 difi = bkpo[ik] - *t;
2584 posb = p;
2585 for (j = i; j <= *n; ++j) {
2586 jk = MAX(m-j,1);
2587 difj = *t - bkpo[jk];
2588 temp = difj * dtdb[posb] + difi * dtdb[posb + 1];
2589 delk = bkpo[ik] - bkpo[jk];
2590 if (delk == 0.)
2591 temp = 0.;
2592 else if (j < *n)
2593 temp /= delk;
2594
2595 dtdb[posb] = temp;
2596 --posb;
2597 }
2598 ++m;
2599 }
2600 }
2601 } else {
2602 fprintf (stderr, "TBSPLN -- Error: T (%f) LIES OUTSIDE OF KNOT DOMAIN [%f; %f] -- ABORT\n", *t, bkpo[1], bkpo[*k + 2]);
2603 *cerr = 50;
2604 }
2605 }
2606
dstorm(int nc,int ns,double * dst,double * dstt,int * dstm,int * u,double * dsdc)2607 void dstorm(int nc, int ns, double *dst, double *dstt, int *dstm, int *u, double *dsdc) {
2608 int i, n, iu;
2609
2610 /* Parameter adjustments */
2611 --u;
2612 --dstm;
2613
2614 for (i = 1; i <= ns; ++i) {
2615 n = dstm[i];
2616 if (n > 0) {
2617 iu = u[i];
2618 r8vset(iu, n, *dst, &dsdc[0]);
2619 r8vset(nc + iu, n, *dstt, &dsdc[0]);
2620 u[i] += n;
2621 }
2622 }
2623 }
2624
fdldc(int grad,int nc,double * dsdc,double * dldc)2625 void fdldc(int grad, int nc, double *dsdc, double *dldc) {
2626 int i__1, i, j;
2627
2628 i = 1;
2629 for (j = 0; j < 3; ++j) {
2630 r8vmul(nc+1, i, nc * 3 + i, nc, &dsdc[0], &dldc[0], &dldc[0]);
2631 i += nc;
2632 }
2633 i = 1;
2634 for (j = 0; j < 3; ++j) {
2635 r8vmul(1, i, i, nc, &dsdc[0], &dldc[0], &dldc[0]);
2636 i += nc;
2637 }
2638 if (grad == 1) {
2639 i = 1;
2640 for (j = 0; j < 9; ++j) {
2641 r8vmul(nc+1, nc * 6 + i, nc * 15 + i, nc, &dsdc[0], &dldc[0], &dldc[0]);
2642 i += nc;
2643 }
2644 i = 1;
2645 for (j = 0; j < 9; ++j) {
2646 i__1 = nc * 6 + i;
2647 r8vmul(1, i__1, i__1, nc, &dsdc[0], &dldc[0], &dldc[0]);
2648 i += nc;
2649 }
2650 }
2651 }
2652
blgen(int grad,int nc,double * b,double * c,double * dldc)2653 void blgen(int grad, int nc, double *b, double *c, double *dldc) {
2654 int i, j;
2655
2656 i = 1;
2657 for (j = 0; j < 6; ++j) {
2658 b[j] += r8sdot(i, 1, nc, &dldc[0], &c[0]);
2659 i += nc;
2660 }
2661 if (grad == 1) {
2662 i = 1;
2663 for (j = 28; j < 46; ++j) {
2664 b[j] += r8sdot(nc * 6 + i, 1, nc, &dldc[0], &c[0]);
2665 i += nc;
2666 }
2667 }
2668 }
2669
bngen_(double * b)2670 void bngen_(double *b) {
2671 double cd, cf, ch, ci, cx, cy, cz, cdt, cft, cht, cit, cxt, cyt, czt;
2672
2673 cx = b[0];
2674 cy = b[1];
2675 cz = b[2];
2676 cxt = b[3];
2677 cyt = b[4];
2678 czt = b[5];
2679 ch = sqrt(cx * cx + cy * cy);
2680 cf = sqrt(ch * ch + cz * cz);
2681 if (ch != 0.) {
2682 cd = atan(cy / (ch + cx)) * 2.;
2683 cdt = (cx * cyt - cy * cxt) / ch / ch;
2684 cht = (cx * cxt + cy * cyt) / ch;
2685 } else {
2686 cd = cdt = cht = 0.;
2687 }
2688 if (cf != 0.) {
2689 ci = atan(cz / (cf + ch)) * 2.;
2690 cit = (ch * czt - cz * cht) / cf / cf;
2691 cft = (ch * cht + cz * czt) / cf;
2692 } else {
2693 ci = cit = cft = 0.;
2694 }
2695 b[6] = cd;
2696 b[7] = ci;
2697 b[8] = ch;
2698 b[9] = cf;
2699 b[10] = cdt;
2700 b[11] = cit;
2701 b[12] = cht;
2702 b[13] = cft;
2703 }
2704
tec(int grad,int k,int nc,double * theta,double * phi,double * b,double * dldc,double * r)2705 void tec(int grad, int k, int nc, double *theta, double *phi, double *b, double *dldc, double *r) {
2706 int i__1;
2707
2708 /* Local variables */
2709 double cphi, sphi, costhe, sinthe;
2710
2711 /* Parameter adjustments */
2712 --b;
2713
2714 if (k >= 1) {
2715 sinthe = sin(*theta);
2716 costhe = cos(*theta);
2717 sincos(*theta, &sinthe, &costhe);
2718 sincos(*phi, &sphi, &cphi);
2719 r[0] = -costhe * cphi;
2720 r[1] = -sphi;
2721 r[2] = -sinthe * cphi;
2722 r[3] = -costhe * sphi;
2723 r[4] = cphi;
2724 r[5] = -sinthe * sphi;
2725 r[6] = sinthe;
2726 r[7] = 0.;
2727 r[8] = -costhe;
2728 ltrans(1, 1, &b[1], &r[0], &b[1]);
2729 ltrans(1, 1, &b[4], &r[0], &b[4]);
2730 ltranv(1, nc, nc, &r[0], &dldc[0]);
2731 ltranv(0, nc, nc, &r[0], &dldc[nc * 3]);
2732 bngen_(&b[1]);
2733 if (grad == 1) {
2734 ltranv(0, 3, 3, &r[0], &b[29]);
2735 ltranv(0, 3, 3, &r[0], &b[38]);
2736 ltrans(1, 1, &b[29], &r[0], &b[29]);
2737 ltrans(1, 1, &b[32], &r[0], &b[32]);
2738 ltrans(1, 1, &b[35], &r[0], &b[35]);
2739 ltrans(1, 1, &b[38], &r[0], &b[38]);
2740 ltrans(1, 1, &b[41], &r[0], &b[41]);
2741 ltrans(1, 1, &b[44], &r[0], &b[44]);
2742 i__1 = nc * 3;
2743 ltranv(0, i__1, i__1, &r[0], &dldc[nc * 6]);
2744 ltranv(0, i__1, i__1, &r[0], &dldc[nc * 15]);
2745 ltranv(0, nc, nc, &r[0], &dldc[nc * 6]);
2746 ltranv(0, nc, nc, &r[0], &dldc[nc * 9]);
2747 ltranv(0, nc, nc, &r[0], &dldc[nc * 12]);
2748 ltranv(0, nc, nc, &r[0], &dldc[nc * 15]);
2749 ltranv(0, nc, nc, &r[0], &dldc[nc * 18]);
2750 ltranv(0, nc, nc, &r[0], &dldc[nc * 21]);
2751 }
2752 }
2753 }
2754
tse(int grad,int k,int nc,double * rse,double * b,double * dldc,double * r)2755 void tse(int grad, int k, int nc, double *rse, double *b, double *dldc, double *r) {
2756 int i__1;
2757
2758 if (k >= 1) {
2759 r8vgathp(1, 1, 1, 9, &rse[0], &r[0]);
2760 ltrans(1, 1, &b[0], &r[0], &b[0]);
2761 ltrans(1, 1, &b[3], &r[0], &b[3]);
2762 ltranv(1, nc, nc, &r[0], &dldc[0]);
2763 ltranv(0, nc, nc, &r[0], &dldc[nc * 3]);
2764 bngen_(&b[0]);
2765 if (grad == 1) {
2766 ltranv(0, 3, 3, &r[0], &b[28]);
2767 ltranv(0, 3, 3, &r[0], &b[37]);
2768 ltrans(1, 1, &b[28], &r[0], &b[28]);
2769 ltrans(1, 1, &b[31], &r[0], &b[31]);
2770 ltrans(1, 1, &b[34], &r[0], &b[34]);
2771 ltrans(1, 1, &b[37], &r[0], &b[37]);
2772 ltrans(1, 1, &b[40], &r[0], &b[40]);
2773 ltrans(1, 1, &b[43], &r[0], &b[43]);
2774 i__1 = nc * 3;
2775 ltranv(0, i__1, i__1, &r[0], &dldc[nc * 6]);
2776 ltranv(0, i__1, i__1, &r[0], &dldc[nc * 15]);
2777 ltranv(0, nc, nc, &r[0], &dldc[nc * 6]);
2778 ltranv(0, nc, nc, &r[0], &dldc[nc * 9]);
2779 ltranv(0, nc, nc, &r[0], &dldc[nc * 12]);
2780 ltranv(0, nc, nc, &r[0], &dldc[nc * 15]);
2781 ltranv(0, nc, nc, &r[0], &dldc[nc * 18]);
2782 ltranv(0, nc, nc, &r[0], &dldc[nc * 21]);
2783 }
2784 }
2785 }
2786
tms(int grad,int k,int nc,int na,int ia,double * a,double * b,double * dldc,double * dlda,double * r)2787 void tms(int grad, int k, int nc, int na, int ia, double *a, double *b, double *dldc, double *dlda, double *r) {
2788 int i__1;
2789 double eulx, euly, eulz, ceulx, ceuly, ceulz, seulx, seuly, seulz;
2790
2791 if (k >= 1) {
2792 eulx = a[ia];
2793 euly = a[ia + 1];
2794 eulz = a[ia + 2];
2795 sincos(eulx, &seulx, &ceulx);
2796 sincos(euly, &seuly, &ceuly);
2797 sincos(eulz, &seulz, &ceulz);
2798 fdldeu_(&k, &na, &ia, &seulx, &ceulx, &seuly, &ceuly, &seulz, &ceulz, &r[0], &b[1], &dlda[0]);
2799 r[0] = ceuly * ceulz;
2800 r[1] = -seulz;
2801 r[2] = -seuly * ceulz;
2802 r[3] = ceuly * ceulx * seulz + seuly * seulx;
2803 r[4] = ceulx * ceulz;
2804 r[5] = ceuly * seulx - seuly * seulz * ceulx;
2805 r[6] = -ceuly * seulx * seulz + seuly * ceulx;
2806 r[7] = -seulx * ceulz;
2807 r[8] = ceuly * ceulx + seuly * seulx * seulz;
2808 ltrans(1, 1, &b[0], &r[0], &b[0]);
2809 ltrans(1, 1, &b[3], &r[0], &b[3]);
2810 ltranv(1, nc, nc, &r[0], &dldc[0]);
2811 ltranv(0, nc, nc, &r[0], &dldc[nc * 3]);
2812 ltranv(0, na, ia, &r[0], &dlda[0]);
2813 ltranv(0, na, ia, &r[0], &dlda[na * 3]);
2814 bngen_(&b[1]);
2815 if (grad == 1) {
2816 ltranv(0, 3, 3, &r[0], &b[28]);
2817 ltranv(0, 3, 3, &r[0], &b[37]);
2818 ltrans(1, 1, &b[28], &r[0], &b[28]);
2819 ltrans(1, 1, &b[31], &r[0], &b[31]);
2820 ltrans(1, 1, &b[34], &r[0], &b[34]);
2821 ltrans(1, 1, &b[37], &r[0], &b[37]);
2822 ltrans(1, 1, &b[40], &r[0], &b[40]);
2823 ltrans(1, 1, &b[43], &r[0], &b[43]);
2824 i__1 = nc * 3;
2825 ltranv(0, i__1, i__1, &r[0], &dldc[nc * 6]);
2826 ltranv(0, i__1, i__1, &r[0], &dldc[nc * 15]);
2827 ltranv(0, nc, nc, &r[0], &dldc[nc * 6]);
2828 ltranv(0, nc, nc, &r[0], &dldc[nc * 9]);
2829 ltranv(0, nc, nc, &r[0], &dldc[nc * 12]);
2830 ltranv(0, nc, nc, &r[0], &dldc[nc * 15]);
2831 ltranv(0, nc, nc, &r[0], &dldc[nc * 18]);
2832 ltranv(0, nc, nc, &r[0], &dldc[nc * 21]);
2833 }
2834 ia += 3;
2835 }
2836 }
2837
fdldeu_(int * k,int * na,int * ia,double * seulx,double * ceulx,double * seuly,double * ceuly,double * seulz,double * ceulz,double * r,double * b,double * dlda)2838 void fdldeu_(int *k, int *na, int *ia, double *seulx, double *ceulx, double *seuly, double *ceuly,
2839 double *seulz, double *ceulz, double *r, double *b, double *dlda) {
2840
2841 int i, j;
2842
2843 if (*k == 1) {
2844 i = *ia;
2845 for (j = 1; j <= 6; ++j) {
2846 dlda[i] = 0.;
2847 dlda[i + 1] = 0.;
2848 dlda[i + 2] = 0.;
2849 i += *na;
2850 }
2851 } else {
2852 r[0] = 0.;
2853 r[1] = 0.;
2854 r[2] = 0.;
2855 r[3] = -(*ceuly) * *seulx * *seulz + *seuly * *ceulx;
2856 r[4] = -(*seulx) * *ceulz;
2857 r[5] = *ceuly * *ceulx + *seuly * *seulz * *seulx;
2858 r[6] = -(*ceuly) * *ceulx * *seulz - *seuly * *seulx;
2859 r[7] = -(*ceulx) * *ceulz;
2860 r[8] = -(*ceuly) * *seulx + *seuly * *ceulx * *seulz;
2861 ltrans(*na, 1, &b[0], &r[0], &dlda[*ia]);
2862 ltrans(*na, 1, &b[3], &r[0], &dlda[*na * 3 + *ia]);
2863 r[0] = -(*seuly) * *ceulz;
2864 r[1] = 0.;
2865 r[2] = -(*ceuly) * *ceulz;
2866 r[3] = -(*seuly) * *ceulx * *seulz + *ceuly * *seulx;
2867 r[4] = 0.;
2868 r[5] = -(*seuly) * *seulx - *ceuly * *seulz * *ceulx;
2869 r[6] = *seuly * *seulx * *seulz + *ceuly * *ceulx;
2870 r[7] = 0.;
2871 r[8] = -(*seuly) * *ceulx + *ceuly * *seulx * *seulz;
2872 ltrans(*na, 1, &b[0], &r[0], &dlda[*ia + 1]);
2873 ltrans(*na, 1, &b[3], &r[0], &dlda[*na * 3 + *ia + 1]);
2874 r[0] = -(*ceuly) * *seulz;
2875 r[1] = -(*ceulz);
2876 r[2] = *seuly * *seulz;
2877 r[3] = *ceuly * *ceulx * *ceulz;
2878 r[4] = -(*ceulx) * *seulz;
2879 r[5] = -(*seuly) * *ceulz * *ceulx;
2880 r[6] = -(*ceuly) * *seulx * *ceulz;
2881 r[7] = *seulx * *seulz;
2882 r[8] = *seuly * *seulx * *ceulz;
2883 ltrans(*na, 1, &b[0], &r[0], &dlda[*ia + 2]);
2884 ltrans(*na, 1, &b[3], &r[0], &dlda[*na * 3 + *ia + 2]);
2885 }
2886 }
2887
tnm_(int * grad,int * k,int * nc,int * na,int * ia,double * a,double * b,double * dldc,double * dlda,double * r)2888 void tnm_(int *grad, int *k, int *nc, int *na, int *ia, double *a, double *b, double *dldc, double *dlda, double *r) {
2889 int i__1;
2890 double chix, chiy, chiz, cchix, cchiy, cchiz, schix, schiy, schiz;
2891
2892 if (*k >= 1) {
2893 chix = a[*ia];
2894 chiy = a[*ia + 1];
2895 chiz = a[*ia + 2];
2896 sincos(chix, &schix, &cchix);
2897 sincos(chiy, &schiy, &cchiy);
2898 sincos(chiz, &schiz, &cchiz);
2899 fdldno_(k, na, ia, &schix, &cchix, &schiy, &cchiy, &schiz, &cchiz, &r[0], &b[0], &dlda[0]);
2900 r[0] = 1.;
2901 r[1] = 0.;
2902 r[2] = 0.;
2903 r[3] = schix;
2904 r[4] = cchix;
2905 r[5] = 0.;
2906 r[6] = schiy * cchiz;
2907 r[7] = schiy * schiz;
2908 r[8] = cchiy;
2909 ltrans(1, 1, &b[0], &r[0], &b[0]);
2910 ltrans(1, 1, &b[3], &r[0], &b[3]);
2911 ltranv(1, *nc, *nc, &r[0], &dldc[0]);
2912 ltranv(0, *nc, *nc, &r[0], &dldc[*nc * 3]);
2913 ltranv(0, *na, *ia, &r[0], &dlda[0]);
2914 ltranv(0, *na, *ia, &r[0], &dlda[*na * 3]);
2915 bngen_(&b[0]);
2916 if (*grad == 1) {
2917 ltranv(0, 3, 3, &r[0], &b[28]);
2918 ltranv(0, 3, 3, &r[0], &b[37]);
2919 ltrans(1, 1, &b[28], &r[0], &b[28]);
2920 ltrans(1, 1, &b[31], &r[0], &b[31]);
2921 ltrans(1, 1, &b[34], &r[0], &b[34]);
2922 ltrans(1, 1, &b[37], &r[0], &b[37]);
2923 ltrans(1, 1, &b[40], &r[0], &b[40]);
2924 ltrans(1, 1, &b[43], &r[0], &b[43]);
2925 i__1 = *nc * 3;
2926 ltranv(0, i__1, i__1, &r[0], &dldc[*nc * 6]);
2927 ltranv(0, i__1, i__1, &r[0], &dldc[*nc * 15]);
2928 ltranv(0, *nc, *nc, &r[0], &dldc[*nc * 6]);
2929 ltranv(0, *nc, *nc, &r[0], &dldc[*nc * 9]);
2930 ltranv(0, *nc, *nc, &r[0], &dldc[*nc * 12]);
2931 ltranv(0, *nc, *nc, &r[0], &dldc[*nc * 15]);
2932 ltranv(0, *nc, *nc, &r[0], &dldc[*nc * 18]);
2933 ltranv(0, *nc, *nc, &r[0], &dldc[*nc * 21]);
2934 }
2935 *ia += 3;
2936 }
2937 }
2938
fdldno_(int * k,int * na,int * ia,double * schix,double * cchix,double * schiy,double * cchiy,double * schiz,double * cchiz,double * r,double * b,double * dlda)2939 void fdldno_(int *k, int *na, int *ia, double *schix, double *cchix, double *schiy, double *cchiy,
2940 double *schiz, double *cchiz, double *r, double *b, double *dlda) {
2941
2942 int i, j;
2943
2944 if (*k == 1) {
2945 i = *ia;
2946 for (j = 1; j <= 6; ++j) {
2947 dlda[i] = 0.;
2948 dlda[i + 1] = 0.;
2949 dlda[i + 2] = 0.;
2950 i += *na;
2951 }
2952 } else {
2953 r[0] = r[1] = r[2] = 0.;
2954 r[3] = *cchix;
2955 r[4] = -(*schix);
2956 r[5] = r[6] = r[7] = r[8] = 0.;
2957 ltrans(*na, 1, &b[0], &r[0], &dlda[*ia]);
2958 ltrans(*na, 1, &b[3], &r[0], &dlda[*na * 3 + *ia]);
2959 r[0] = r[1] = r[2] = r[3] = r[4] = r[5] = 0.;
2960 r[6] = *cchiy * *cchiz;
2961 r[7] = *cchiy * *schiz;
2962 r[8] = -(*schiy);
2963 ltrans(*na, 1, &b[0], &r[0], &dlda[*ia + 1]);
2964 ltrans(*na, 1, &b[3], &r[0], &dlda[*na * 3 + *ia + 1]);
2965 r[0] = r[1] = r[2] = r[3] = r[4] = r[5] = 0.;
2966 r[6] = -(*schiy) * *schiz;
2967 r[7] = *schiy * *cchiz;
2968 r[8] = 0.;
2969 ltrans(*na, 1, &b[0], &r[0], &dlda[*ia + 2]);
2970 ltrans(*na, 1, &b[3], &r[0], &dlda[*na * 3 + *ia + 2]);
2971 }
2972 }
2973
tvn_(int * grad,int * k,int * nc,int * na,int * ia,double * a,double * b,double * dldc,double * dlda,double * r)2974 void tvn_(int *grad, int *k, int *nc, int *na, int *ia, double *a, double *b, double *dldc, double *dlda, double *r) {
2975 int i__1;
2976
2977 /* Local variables */
2978 double slpx, slpy, slpz;
2979
2980 if (*k >= 1) {
2981 slpx = a[*ia];
2982 slpy = a[*ia + 1];
2983 slpz = a[*ia + 2];
2984 fdldsl_(k, na, ia, &b[0], &dlda[0]);
2985 r[0] = slpx;
2986 r[1] = 0.;
2987 r[2] = 0.;
2988 r[3] = 0.;
2989 r[4] = slpy;
2990 r[5] = 0.;
2991 r[6] = 0.;
2992 r[7] = 0.;
2993 r[8] = slpz;
2994 ltrans(1, 1, &b[0], &r[0], &b[0]);
2995 ltrans(1, 1, &b[3], &r[0], &b[3]);
2996 ltranv(1, *nc, *nc, &r[0], &dldc[0]);
2997 ltranv(0, *nc, *nc, &r[0], &dldc[*nc * 3]);
2998 ltranv(0, *na, *ia, &r[0], &dlda[0]);
2999 ltranv(0, *na, *ia, &r[0], &dlda[*na * 3]);
3000 bngen_(&b[0]);
3001 if (*grad == 1) {
3002 ltranv(0, 3, 3, &r[0], &b[28]);
3003 ltranv(0, 3, 3, &r[0], &b[37]);
3004 ltrans(1, 1, &b[28], &r[0], &b[28]);
3005 ltrans(1, 1, &b[31], &r[0], &b[31]);
3006 ltrans(1, 1, &b[34], &r[0], &b[34]);
3007 ltrans(1, 1, &b[37], &r[0], &b[37]);
3008 ltrans(1, 1, &b[40], &r[0], &b[40]);
3009 ltrans(1, 1, &b[43], &r[0], &b[43]);
3010 i__1 = *nc * 3;
3011 ltranv(0, i__1, i__1, &r[0], &dldc[*nc * 6]);
3012 ltranv(0, i__1, i__1, &r[0], &dldc[*nc * 15]);
3013 ltranv(0, *nc, *nc, &r[0], &dldc[*nc * 6]);
3014 ltranv(0, *nc, *nc, &r[0], &dldc[*nc * 9]);
3015 ltranv(0, *nc, *nc, &r[0], &dldc[*nc * 12]);
3016 ltranv(0, *nc, *nc, &r[0], &dldc[*nc * 15]);
3017 ltranv(0, *nc, *nc, &r[0], &dldc[*nc * 18]);
3018 ltranv(0, *nc, *nc, &r[0], &dldc[*nc * 21]);
3019 }
3020 *ia += 3;
3021 }
3022 }
3023
fdldsl_(int * k,int * na,int * ia,double * b,double * dlda)3024 void fdldsl_(int *k, int *na, int *ia, double *b, double *dlda) {
3025 int i, j;
3026
3027 i = *ia;
3028 for (j = 0; j < 6; ++j) {
3029 dlda[i] = 0.;
3030 dlda[i + 1] = 0.;
3031 dlda[i + 2] = 0.;
3032 i += *na;
3033 }
3034 if (*k > 1) {
3035 i = *ia;
3036 dlda[i] = b[0];
3037 i += *na;
3038 dlda[i + 1] = b[1];
3039 i += *na;
3040 dlda[i + 2] = b[2];
3041 i += *na;
3042 dlda[i] = b[3];
3043 i += *na;
3044 dlda[i + 1] = b[4];
3045 i += *na;
3046 dlda[i + 2] = b[5];
3047 }
3048 }
3049
tbi_(int * k,int * na,int * ia,double * a,double * b,double * dlda)3050 void tbi_(int *k, int *na, int *ia, double *a, double *b, double *dlda) {
3051 double biax, biay, biaz;
3052
3053 if (*k >= 1) {
3054 biax = a[*ia];
3055 biay = a[*ia + 1];
3056 biaz = a[*ia + 2];
3057 fdldbi_(k, na, ia, &dlda[0]);
3058 b[0] += biax;
3059 b[1] += biay;
3060 b[2] += biaz;
3061 bngen_(&b[0]);
3062 *ia += 3;
3063 }
3064 }
3065
fdldbi_(int * k,int * na,int * ia,double * dlda)3066 void fdldbi_(int *k, int *na, int *ia, double *dlda) {
3067 int i, j;
3068
3069 i = *ia;
3070 for (j = 0; j < 6; ++j) {
3071 dlda[i] = 0.;
3072 dlda[i + 1] = 0.;
3073 dlda[i + 2] = 0.;
3074 i += *na;
3075 }
3076 if (*k > 1) {
3077 i = *ia;
3078 dlda[i] = 1.;
3079 i += *na;
3080 dlda[i + 1] = 1.;
3081 i += *na;
3082 dlda[i + 2] = 1.;
3083 }
3084 }
3085
ltrans(int n,int m,double * q,double * r,double * s)3086 void ltrans(int n, int m, double *q, double *r, double *s) {
3087 double q3;
3088
3089 q3 = q[m + m];
3090 s[0] = q[0] * r[0] + q[m] * r[1] + q3 * r[2];
3091 s[n] = q[0] * r[3] + q[m] * r[4] + q3 * r[5];
3092 s[n + n] = q[0] * r[6] + q[m] * r[7] + q3 * r[8];
3093 }
3094
ltranv(int rfac,int n,int m,double * r,double * v)3095 void ltranv(int rfac, int n, int m, double *r, double *v) {
3096 /* System generated locals */
3097 int i__1;
3098
3099 /* Parameter adjustments */
3100 --r;
3101
3102 /* Function Body */
3103 if (m > 0) {
3104 if (rfac == 1) {
3105 r[10] = r[4] / r[1];
3106 r[11] = r[5] - r[2] * r[10];
3107 r[12] = r[6] - r[3] * r[10];
3108 r[13] = r[7] / r[1];
3109 r[14] = (r[8] - r[2] * r[13]) / r[11];
3110 r[15] = r[9] - r[3] * r[13] - r[12] * r[14];
3111 r[13] -= r[10] * r[14];
3112 }
3113 r8vscale(1, m, r[1], &v[0]);
3114 r8vlinkt(n+1, 1, m, r[2], &v[0], &v[0]);
3115 r8vlinkt(n + n + 1, 1, m, r[3], &v[0], &v[0]);
3116 i__1 = n + 1;
3117 r8vscale(i__1, m, r[11], &v[0]);
3118 r8vlinkt(1, i__1, m, r[10], &v[0], &v[0]);
3119 i__1 = n + n + 1;
3120 r8vlinkt(i__1, n+1, m, r[12], &v[0], &v[0]);
3121 r8vscale(i__1, m, r[15], &v[0]);
3122 r8vlinkt(1, i__1, m, r[13], &v[0], &v[0]);
3123 r8vlinkt(n+1, i__1, m, r[14], &v[0], &v[0]);
3124 }
3125 }
3126
nshx(int nmax,int nmin,int mmax,int mmin)3127 int nshx(int nmax, int nmin, int mmax, int mmin) {
3128 int ret_val, i__2, i__5, i__6, kmax;
3129
3130 kmax = mmax + 1;
3131 i__5 = MIN(nmin,mmin);
3132 i__6 = MIN(nmin,kmax);
3133 i__2 = kmax * kmax - mmin * mmin + i__5 * i__5 - i__6 * i__6 + (nmax - mmax - I_DIM(nmin, kmax)) *
3134 ((mmax << 1) + 1) + (I_DIM(nmin, mmin) - nmax + mmin - 1) * MAX(0, (mmin << 1) - 1);
3135 ret_val = MAX(0,i__2);
3136 return ret_val;
3137 }
3138
nlpx(int nmax,int mmax,int mmin)3139 int nlpx(int nmax, int mmax, int mmin) {
3140 int mdif;
3141
3142 mdif = MAX(0, MIN(nmax,mmax) - mmin + 1);
3143 return(mdif * (mdif + 1) / 2 + mdif * I_DIM(nmax, mmax) + mmin - 1);
3144 }
3145
i8ssum(int abeg,int alen,int * a)3146 int i8ssum(int abeg, int alen, int *a) {
3147 int i, aadr, ret_val;
3148
3149 /* Parameter adjustments */
3150 --a;
3151
3152 ret_val = 0;
3153 aadr = abeg;
3154 for (i = 0; i < alen; ++i)
3155 ret_val += a[aadr++];
3156
3157 return ret_val;
3158 }
3159
i8vset(int abeg,int alen,int s,int * a)3160 void i8vset(int abeg, int alen, int s, int *a) {
3161 int i, aadr;
3162
3163 /* Parameter adjustments */
3164 --a;
3165
3166 /* Function Body */
3167 aadr = abeg;
3168 for (i = 0; i < alen; ++i)
3169 a[aadr++] = s;
3170 }
3171
i8vadd(int abeg,int bbeg,int cbeg,int vlen,int * a,int * b,int * c)3172 void i8vadd(int abeg, int bbeg, int cbeg, int vlen, int *a, int *b, int *c) {
3173 int i, aadr, badr, cadr;
3174
3175 /* Parameter adjustments */
3176 --c;
3177 --b;
3178 --a;
3179
3180 /* Function Body */
3181 aadr = abeg;
3182 badr = bbeg;
3183 cadr = cbeg;
3184 for (i = 0; i < vlen; ++i)
3185 c[cadr++] = b[badr++] + a[aadr++];
3186 }
3187
i8vadds(int abeg,int bbeg,int vlen,int s,int * a,int * b)3188 void i8vadds(int abeg, int bbeg, int vlen, int s, int *a, int *b) {
3189 int i, aadr, badr;
3190
3191 /* Parameter adjustments */
3192 --b;
3193 --a;
3194
3195 /* Function Body */
3196 aadr = abeg;
3197 badr = bbeg;
3198 for (i = 0; i < vlen; ++i)
3199 b[badr++] = a[aadr++] + s;
3200 }
3201
i8vcum(int abas,int abeg,int alen,int * a)3202 void i8vcum(int abas, int abeg, int alen, int *a) {
3203 int i, aadr, acur, aprv;
3204
3205 /* Parameter adjustments */
3206 --a;
3207
3208 aprv = a[abeg];
3209 a[abeg] = abas;
3210 aadr = abeg + 1;
3211 for (i = 0; i < alen - 1; ++i) {
3212 acur = a[aadr];
3213 a[aadr] = a[aadr - 1] + aprv;
3214 aprv = acur;
3215 ++aadr;
3216 }
3217 }
3218
i8vdel(int abas,int abeg,int alen,int * a)3219 void i8vdel(int abas, int abeg, int alen, int *a) {
3220 int i, aadr, acur, aprv;
3221
3222 /* Parameter adjustments */
3223 --a;
3224
3225 aprv = abas;
3226 aadr = abeg;
3227 for (i = 0; i < alen; ++i) {
3228 acur = a[aadr];
3229 a[aadr] = acur - aprv;
3230 aprv = acur;
3231 ++aadr;
3232 }
3233 }
3234
r8vset(int abeg,int alen,double s,double * a)3235 void r8vset(int abeg, int alen, double s, double *a) {
3236 int i, aadr;
3237
3238 /* Parameter adjustments */
3239 --a;
3240
3241 aadr = abeg;
3242 for (i = 0; i < alen; ++i)
3243 a[aadr++] = s;
3244 }
3245
r8sdot(int abeg,int bbeg,int vlen,double * a,double * b)3246 double r8sdot(int abeg, int bbeg, int vlen, double *a, double *b) {
3247 double ret_val;
3248 int i, aadr, badr;
3249
3250 /* Parameter adjustments */
3251 --b;
3252 --a;
3253
3254 ret_val = 0.;
3255 aadr = abeg;
3256 badr = bbeg;
3257 for (i = 0; i < vlen; ++i)
3258 ret_val += a[aadr++] * b[badr++];
3259
3260 return ret_val;
3261 }
3262
r8ssum_(int * abeg,int * alen,double * a)3263 double r8ssum_(int *abeg, int *alen, double *a) {
3264 double ret_val;
3265 int i, aadr;
3266
3267 /* Parameter adjustments */
3268 --a;
3269
3270 /* Function Body */
3271 ret_val = 0.;
3272 aadr = *abeg;
3273 for (i = 0; i < *alen; ++i) {
3274 ret_val += a[aadr];
3275 ++aadr;
3276 }
3277 return ret_val;
3278 }
3279
r8slt(int abeg,int alen,double s,double * a,int * j)3280 void r8slt(int abeg, int alen, double s, double *a, int *j) {
3281 int i, aadr;
3282
3283 /* Parameter adjustments */
3284 --a;
3285
3286 /* Function Body */
3287 aadr = abeg;
3288 for (i = 1; i <= alen; ++i) {
3289 if (s < a[aadr]) {
3290 *j = i - 1;
3291 return;
3292 }
3293 ++aadr;
3294 }
3295 *j = alen;
3296 }
3297
r8vsub(int abeg,int bbeg,int cbeg,int vlen,double * a,double * b,double * c)3298 void r8vsub(int abeg, int bbeg, int cbeg, int vlen, double *a, double *b, double *c) {
3299 int i, aadr, badr, cadr;
3300
3301 /* Parameter adjustments */
3302 --c;
3303 --b;
3304 --a;
3305
3306 /* Function Body */
3307 aadr = abeg;
3308 badr = bbeg;
3309 cadr = cbeg;
3310 for (i = 0; i < vlen; ++i)
3311 c[cadr++] = b[badr++] - a[aadr++];
3312 }
3313
r8vmul(int abeg,int bbeg,int cbeg,int vlen,double * a,double * b,double * c)3314 void r8vmul(int abeg, int bbeg, int cbeg, int vlen, double *a, double *b, double *c) {
3315 int i, aadr, badr, cadr;
3316
3317 /* Parameter adjustments */
3318 --c;
3319 --b;
3320 --a;
3321
3322 aadr = abeg;
3323 badr = bbeg;
3324 cadr = cbeg;
3325 for (i = 0; i < vlen; ++i)
3326 c[cadr++] = b[badr++] * a[aadr++];
3327 }
3328
r8vscale(int abeg,int alen,double s,double * a)3329 void r8vscale(int abeg, int alen, double s, double *a) {
3330 int i, aadr;
3331
3332 /* Parameter adjustments */
3333 --a;
3334
3335 aadr = abeg;
3336 for (i = 0; i < alen; ++i) {
3337 a[aadr] *= s;
3338 ++aadr;
3339 }
3340 }
3341
r8vscats(int qbeg,int qlen,double s,int * q,double * a)3342 void r8vscats(int qbeg, int qlen, double s, int *q, double *a) {
3343 int i, qadr;
3344
3345 /* Parameter adjustments */
3346 --a;
3347 --q;
3348
3349 qadr = qbeg;
3350 for (i = 0; i < qlen; ++i)
3351 a[q[qadr++]] = s;
3352
3353 }
3354
r8vlinkt(int abeg,int bbeg,int vlen,double s,double * a,double * b)3355 void r8vlinkt(int abeg, int bbeg, int vlen, double s, double *a, double *b) {
3356 int i, aadr, badr;
3357
3358 /* Parameter adjustments */
3359 --b;
3360 --a;
3361
3362 aadr = abeg;
3363 badr = bbeg;
3364 for (i = 0; i < vlen; ++i)
3365 b[badr++] += s * a[aadr++];
3366 }
3367
r8vlinkq(int abeg,int bbeg,int cbeg,int vlen,double s,double * a,double * b,double * c)3368 void r8vlinkq(int abeg, int bbeg, int cbeg, int vlen, double s, double *a, double *b, double *c) {
3369 int i, aadr, badr, cadr;
3370
3371 /* Parameter adjustments */
3372 --c;
3373 --b;
3374 --a;
3375
3376 aadr = abeg;
3377 badr = bbeg;
3378 cadr = cbeg;
3379 for (i = 0; i < vlen; ++i)
3380 c[cadr++] += s * a[aadr++] * b[badr++];
3381 }
3382
r8vgathp(int abeg,int ainc,int bbeg,int blen,double * a,double * b)3383 void r8vgathp(int abeg, int ainc, int bbeg, int blen, double *a, double *b) {
3384 int i, aadr, badr;
3385
3386 /* Parameter adjustments */
3387 --b; --a;
3388
3389 aadr = abeg;
3390 badr = bbeg;
3391 for (i = 0; i < blen; ++i) {
3392 b[badr++] = a[aadr];
3393 aadr += ainc;
3394 }
3395 }
3396
d_mod(double x,double y)3397 double d_mod(double x, double y) {
3398 double quotient;
3399 if( (quotient = x / y) >= 0)
3400 quotient = floor(quotient);
3401 else
3402 quotient = -floor(-quotient);
3403 return(x - y * quotient );
3404 }
3405
pow_di(double ap,int bp)3406 double pow_di(double ap, int bp) {
3407 double pow, x;
3408 int n;
3409 unsigned long u;
3410
3411 pow = 1;
3412 x = ap;
3413 n = bp;
3414
3415 if(n != 0) {
3416 if(n < 0) {
3417 n = -n;
3418 x = 1/x;
3419 }
3420 for(u = n; ; ) {
3421 if(u & 01)
3422 pow *= x;
3423 if(u >>= 1)
3424 x *= x;
3425 else
3426 break;
3427 }
3428 }
3429 return(pow);
3430 }
3431
i_dnnt(double x)3432 int i_dnnt(double x) {
3433 return (int)(x >= 0. ? floor(x + 0.5) : -floor(0.5 - x));
3434 }
3435