1 /*
2 * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
3 * Applied Mathematics, Norway.
4 *
5 * Contact information: E-mail: tor.dokken@sintef.no
6 * SINTEF ICT, Department of Applied Mathematics,
7 * P.O. Box 124 Blindern,
8 * 0314 Oslo, Norway.
9 *
10 * This file is part of SISL.
11 *
12 * SISL is free software: you can redistribute it and/or modify
13 * it under the terms of the GNU Affero General Public License as
14 * published by the Free Software Foundation, either version 3 of the
15 * License, or (at your option) any later version.
16 *
17 * SISL is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Affero General Public License for more details.
21 *
22 * You should have received a copy of the GNU Affero General Public
23 * License along with SISL. If not, see
24 * <http://www.gnu.org/licenses/>.
25 *
26 * In accordance with Section 7(b) of the GNU Affero General Public
27 * License, a covered work must retain the producer line in every data
28 * file that is created or manipulated using SISL.
29 *
30 * Other Usage
31 * You can be released from the requirements of the license by purchasing
32 * a commercial license. Buying such a license is mandatory as soon as you
33 * develop commercial activities involving the SISL library without
34 * disclosing the source code of your own applications.
35 *
36 * This file may be used in accordance with the terms contained in a
37 * written agreement between you and SINTEF ICT.
38 */
39
40 #include "sisl-copyright.h"
41
42 /*
43 *
44 *
45 *
46 */
47
48
49 #define S1516
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
s1516(double ep[],double epar[],int im,int idim,double ** ev,int * jstat)54 void s1516(double ep[],double epar[],int im,int idim,
55 double **ev,int *jstat)
56 #else
57 void s1516(ep,epar,im,idim,ev,jstat)
58 double ep[];
59 double epar[];
60 int im;
61 int idim;
62 double **ev;
63 int *jstat;
64 #endif
65 /*
66 ************************************************************************
67 *
68 * Purpose: To estimate the first derivative at each point in a sequence.
69 *
70 * Input:
71 * ep - Array containing the point in sequence
72 * (x,y,..,x,y,..), length idim * im.
73 * epar - Parametrization array. The array should be increasing
74 * in value.
75 * im - Number of point and derivatives
76 * idim - The dimension of the space the points and derivatives
77 * lie in
78 * Output:
79 * ev - Pointer to array containing the derivatives in sequence
80 * (x,y,..,x,y,..), length idim * im.
81 * jstat - Status variable
82 * < 0 - Error.
83 * Method:
84 *
85 * REFERENCES :
86 *
87 *-
88 * CALLS :
89 *
90 * WRITTEN BY : Michael Floater, SI 1993-10
91 *
92 *********************************************************************
93 */
94 {
95 int ki,kj; /* Loop variables */
96 int kk; /* Polynomial order */
97 int kpos=0; /* Position of error */
98 int kstat=0; /* Status variable */
99 double *gpar;
100 int kcnsta;
101 int kcnend;
102 int iopen;
103 int iorder;
104 int ileft;
105 double *ntype;
106 SISLCurve *qc;
107 int knbpar;
108 double *evtemp;
109 double cendpar;
110 double *eder;
111
112
113
114
115 /* Check input */
116
117 if (idim < 1 || im < 2) goto err102;
118
119
120 /* Allocate array for derivatives */
121
122 evtemp = newarray(idim*im,DOUBLE);
123 if (evtemp == SISL_NULL) goto err101;
124
125 ntype = newarray(im,DOUBLE);
126 if (ntype == SISL_NULL) goto err101;
127
128 for(ki=0; ki<im; ki++)
129 {
130 ntype[ki] = 1.0;
131 }
132
133 eder = newarray(2 * idim,DOUBLE);
134 if (eder == SISL_NULL) goto err101;
135
136
137
138 kcnsta = 1;
139 kcnend = 1;
140 iopen = 1;
141 iorder = 4;
142
143 s1358(ep, im, idim, ntype, epar, kcnsta, kcnend, iopen, iorder,
144 epar[0],&cendpar, &qc, &gpar, &knbpar, &kstat);
145 if(kstat < 0) goto error;
146
147 for(ki=0,kk=0; ki<im; ki++,kk+=idim)
148 {
149 s1221(qc,1,epar[ki],&ileft,eder,&kstat);
150 if(kstat < 0) goto error;
151
152 for(kj=0; kj<idim; kj++)
153 {
154 evtemp[kk+kj] = eder[idim+kj];
155 }
156 }
157
158
159 /* Calculation completed */
160
161 /* Set result. */
162
163 (*ev) = evtemp;
164
165 *jstat = 0;
166 goto out;
167
168
169
170 /* Error in space allocation */
171
172 err101: *jstat = -101;
173 s6err("s1516",*jstat,kpos);
174 goto out;
175
176
177 /* Error in input. */
178
179 err102: *jstat = -102;
180 s6err("s1516",*jstat,kpos);
181 goto out;
182
183 /* Error in lower level routine. */
184
185 error: *jstat =kstat;
186 s6err("s1516",*jstat,kpos);
187 goto out;
188
189 out:
190 if (ntype != SISL_NULL) freearray(ntype);
191 if (gpar != SISL_NULL) freearray(gpar);
192 if (eder != SISL_NULL) freearray(eder);
193
194 return;
195 }
196