1 #include "splines/oskar_dierckx_sphere.h"
2 #include "splines/oskar_dierckx_fpsphe.h"
3 #include <math.h>
4
5 #ifdef __cplusplus
6 extern "C" {
7 #endif
8
oskar_dierckx_sphere(int iopt,int m,const double * theta,const double * phi,const double * r,const double * w,double s,int ntest,int npest,double eps,int * nt,double * tt,int * np,double * tp,double * c,double * fp,double * wrk1,int lwrk1,double * wrk2,int lwrk2,int * iwrk,int kwrk,int * ier)9 void oskar_dierckx_sphere(int iopt, int m, const double* theta,
10 const double* phi, const double* r, const double* w, double s,
11 int ntest, int npest, double eps, int* nt, double* tt, int* np,
12 double* tp, double* c, double* fp, double* wrk1, int lwrk1,
13 double* wrk2, int lwrk2, int* iwrk, int kwrk, int* ier)
14 {
15 /* Local variables */
16 int i = 0, j = 0, la = 0, lf = 0, ki = 0, lh = 0, kn = 0;
17 int lq = 0, ib1 = 0, ib3 = 0;
18 int np4 = 0, nt4 = 0, lcc = 0, ncc = 0, lff = 0, lco = 0;
19 int lbp = 0, lbt = 0, lcs = 0, lfp = 0, lro = 0, npp = 0;
20 int lsp = 0, lst = 0, ntt = 0, ncof = 0, nreg = 0, ncest = 0;
21 int maxit = 0, nrint = 0, kwest = 0, lwest = 0;
22 double tol = 0.0, pi = 0.0, pi2 = 0.0;
23
24 /* Parameter adjustments */
25 --tt;
26 --tp;
27 --wrk1;
28 --wrk2;
29 --iwrk;
30
31 /* Function Body */
32 /* we set up the parameters tol and maxit. */
33 maxit = 20;
34 tol = 0.001;
35 /* before starting computations a data check is made. if the input data
36 * are invalid, control is immediately repassed to the calling program. */
37 *ier = 10;
38 if (eps <= 0.0 || eps >= 1.0) return;
39 if (iopt < -1 || iopt > 1) return;
40 if (m < 2) return;
41 if (ntest < 8 || npest < 8) return;
42 nt4 = ntest - 4;
43 np4 = npest - 4;
44 ncest = nt4 * np4;
45 ntt = ntest - 7;
46 npp = npest - 7;
47 ncc = npp * (ntt - 1) + 6;
48 nrint = ntt + npp;
49 nreg = ntt * npp;
50 ncof = npp * 3 + 6;
51 ib1 = npp << 2;
52 ib3 = ib1 + 3;
53 if (ncof > ib1) ib1 = ncof;
54 if (ncof > ib3) ib3 = ncof;
55 lwest = 185 + 52 * npp + 10 * ntt + 14 * ntt * npp + 8 * (m + (ntt - 1) *
56 (npp * npp));
57 kwest = m + nreg;
58 if (lwrk1 < lwest || kwrk < kwest) return;
59 pi = 4.0 * atan(1.0);
60 pi2 = pi + pi;
61 if (iopt >= 0 && s < 0.0) return;
62 if (iopt <= 0)
63 {
64 for (i = 0; i < m; ++i)
65 {
66 if (w[i] <= 0.0) return;
67 if (theta[i] < 0.0 || theta[i] > pi) return;
68 if (phi[i] < 0.0 || phi[i] > pi2) return;
69 }
70 if (iopt < 0)
71 {
72 ntt = *nt - 8;
73 if (ntt < 0 || *nt > ntest) return;
74 if (ntt != 0)
75 {
76 tt[4] = 0.0;
77 for (i = 1; i <= ntt; ++i)
78 {
79 j = i + 4;
80 if (tt[j] <= tt[j - 1] || tt[j] >= pi) return;
81 }
82 }
83 npp = *np - 8;
84 if (npp < 1 || *np > npest) return;
85 tp[4] = 0.0;
86 for (i = 1; i <= npp; ++i)
87 {
88 j = i + 4;
89 if (tp[j] <= tp[j - 1] || tp[j] >= pi2) return;
90 }
91 }
92 }
93 *ier = 0;
94 /* we partition the working space and determine the spline approximation */
95 kn = 1;
96 ki = kn + m;
97 lq = 2;
98 la = lq + ncc * ib3;
99 lf = la + ncc * ib1;
100 lff = lf + ncc;
101 lfp = lff + ncest;
102 lco = lfp + nrint;
103 lh = lco + nrint;
104 lbt = lh + ib3;
105 lbp = lbt + ntest * 5;
106 lro = lbp + npest * 5;
107 lcc = lro + npest;
108 lcs = lcc + npest;
109 lst = lcs + npest;
110 lsp = lst + (m << 2);
111 oskar_dierckx_fpsphe(iopt, m, theta, phi, r, w, s, ntest, npest, eps,
112 tol, maxit, ncc, nt, &tt[1], np, &tp[1], c, fp, &wrk1[1],
113 &wrk1[lfp], &wrk1[lco], &wrk1[lf], &wrk1[lff], &wrk1[lro],
114 &wrk1[lcc], &wrk1[lcs], &wrk1[la], &wrk1[lq], &wrk1[lbt],
115 &wrk1[lbp], &wrk1[lst], &wrk1[lsp], &wrk1[lh], &iwrk[ki],
116 &iwrk[kn], &wrk2[1], lwrk2, ier);
117 }
118
119 #ifdef __cplusplus
120 }
121 #endif
122