API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Jun 29th, 2014
336
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 14.02 KB | None | 0 0
raw download clone embed print report
  1.  
  2.  
  3. #define FALSE 0
  4. #define TRUE 1
  5.  
  6. //Define this if using with UFRaw
  7. //#define __WITH_UFRAW__
  8.  
  9.  
  10.  
  11. #include "nikon_curve.h"
  12. #include <stdio.h>
  13. #define MAX(a,b) ((a) > (b) ? (a) : (b))
  14. #define MIN(a,b) ((a) < (b) ? (a) : (b))
  15. #define g_fopen fopen
  16.  
  17. #include <stdlib.h>
  18. #include <string.h>
  19. #include <math.h>
  20. #include <errno.h>
  21. #include <stdarg.h> /* For variable argument lists */
  22.  
  23.  
  24. #include "nikon_curve.h"
  25.  
  26.  
  27. static void nc_merror(void *ptr, char *where)
  28. {
  29. if (ptr) return;
  30.  
  31. fprintf(stderr, "Out of memory in %s\n", where);
  32. exit(1);
  33.  
  34. }
  35.  
  36.  
  37.  
  38.  
  39. static void nc_message(int code, char *format, ...)
  40. {
  41. char message[256];
  42. va_list ap;
  43. va_start(ap, format);
  44.  
  45. vsnprintf(message, 255, format, ap);
  46. message[255] = '\0';
  47. va_end(ap);
  48.  
  49.  
  50. code = code;
  51. fprintf(stderr, "%s", message);
  52. fflush(stderr);
  53.  
  54.  
  55. }
  56.  
  57. //**********************************************************************
  58. //
  59. // Purpose:
  60. //
  61. // D3_NP_FS factors and solves a D3 system.
  62. //
  63. // Discussion:
  64. //
  65. // The D3 storage format is used for a tridiagonal matrix.
  66. // The superdiagonal is stored in entries (1,2:N), the diagonal in
  67. // entries (2,1:N), and the subdiagonal in (3,1:N-1). Thus, the
  68. // original matrix is "collapsed" vertically into the array.
  69. //
  70. // This algorithm requires that each diagonal entry be nonzero.
  71. // It does not use pivoting, and so can fail on systems that
  72. // are actually nonsingular.
  73. //
  74. // Example:
  75. //
  76. // Here is how a D3 matrix of order 5 would be stored:
  77. //
  78. // * A12 A23 A34 A45
  79. // A11 A22 A33 A44 A55
  80. // A21 A32 A43 A54 *
  81. //
  82. // Modified:
  83. //
  84. // 07 January 2005 Shawn Freeman (pure C modifications)
  85. // 15 November 2003 John Burkardt
  86. //
  87. // Author:
  88. //
  89. // John Burkardt
  90. //
  91. // Parameters:
  92. //
  93. // Input, int N, the order of the linear system.
  94. //
  95. // Input/output, double A[3*N].
  96. // On input, the nonzero diagonals of the linear system.
  97. // On output, the data in these vectors has been overwritten
  98. // by factorization information.
  99. //
  100. // Input, double B[N], the right hand side.
  101. //
  102. // Output, double D3_NP_FS[N], the solution of the linear system.
  103. // This is NULL if there was an error because one of the diagonal
  104. // entries was zero.
  105. //
  106. static double *d3_np_fs(int n, double a[], double b[])
  107. {
  108. int i;
  109. double *x;
  110. double xmult;
  111. //
  112. // Check.
  113. //
  114. for (i = 0; i < n; i++) {
  115. if (a[1 + i * 3] == 0.0E+00) {
  116. return NULL;
  117. }
  118. }
  119. x = (double *)calloc(n, sizeof(double));
  120. nc_merror(x, "d3_np_fs");
  121.  
  122. for (i = 0; i < n; i++) {
  123. x[i] = b[i];
  124. }
  125.  
  126. for (i = 1; i < n; i++) {
  127. xmult = a[2 + (i - 1) * 3] / a[1 + (i - 1) * 3];
  128. a[1 + i * 3] = a[1 + i * 3] - xmult * a[0 + i * 3];
  129. x[i] = x[i] - xmult * x[i - 1];
  130. }
  131.  
  132. x[n - 1] = x[n - 1] / a[1 + (n - 1) * 3];
  133. for (i = n - 2; 0 <= i; i--) {
  134. x[i] = (x[i] - a[0 + (i + 1) * 3] * x[i + 1]) / a[1 + i * 3];
  135. }
  136.  
  137. return x;
  138. }
  139.  
  140.  
  141. //**********************************************************************
  142. //
  143. // Purpose:
  144. //
  145. // SPLINE_CUBIC_SET computes the second derivatives of a piecewise cubic spline.
  146. //
  147. // Discussion:
  148. //
  149. // For data interpolation, the user must call SPLINE_SET to determine
  150. // the second derivative data, passing in the data to be interpolated,
  151. // and the desired boundary conditions.
  152. //
  153. // The data to be interpolated, plus the SPLINE_SET output, defines
  154. // the spline. The user may then call SPLINE_VAL to evaluate the
  155. // spline at any point.
  156. //
  157. // The cubic spline is a piecewise cubic polynomial. The intervals
  158. // are determined by the "knots" or abscissas of the data to be
  159. // interpolated. The cubic spline has continous first and second
  160. // derivatives over the entire interval of interpolation.
  161. //
  162. // For any point T in the interval T(IVAL), T(IVAL+1), the form of
  163. // the spline is
  164. //
  165. // SPL(T) = A(IVAL)
  166. // + B(IVAL) * ( T - T(IVAL) )
  167. // + C(IVAL) * ( T - T(IVAL) )**2
  168. // + D(IVAL) * ( T - T(IVAL) )**3
  169. //
  170. // If we assume that we know the values Y(*) and YPP(*), which represent
  171. // the values and second derivatives of the spline at each knot, then
  172. // the coefficients can be computed as:
  173. //
  174. // A(IVAL) = Y(IVAL)
  175. // B(IVAL) = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
  176. // - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
  177. // C(IVAL) = YPP(IVAL) / 2
  178. // D(IVAL) = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
  179. //
  180. // Since the first derivative of the spline is
  181. //
  182. // SPL'(T) = B(IVAL)
  183. // + 2 * C(IVAL) * ( T - T(IVAL) )
  184. // + 3 * D(IVAL) * ( T - T(IVAL) )**2,
  185. //
  186. // the requirement that the first derivative be continuous at interior
  187. // knot I results in a total of N-2 equations, of the form:
  188. //
  189. // B(IVAL-1) + 2 C(IVAL-1) * (T(IVAL)-T(IVAL-1))
  190. // + 3 * D(IVAL-1) * (T(IVAL) - T(IVAL-1))**2 = B(IVAL)
  191. //
  192. // or, setting H(IVAL) = T(IVAL+1) - T(IVAL)
  193. //
  194. // ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
  195. // - ( YPP(IVAL) + 2 * YPP(IVAL-1) ) * H(IVAL-1) / 6
  196. // + YPP(IVAL-1) * H(IVAL-1)
  197. // + ( YPP(IVAL) - YPP(IVAL-1) ) * H(IVAL-1) / 2
  198. // =
  199. // ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
  200. // - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * H(IVAL) / 6
  201. //
  202. // or
  203. //
  204. // YPP(IVAL-1) * H(IVAL-1) + 2 * YPP(IVAL) * ( H(IVAL-1) + H(IVAL) )
  205. // + YPP(IVAL) * H(IVAL)
  206. // =
  207. // 6 * ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
  208. // - 6 * ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
  209. //
  210. // Boundary conditions must be applied at the first and last knots.
  211. // The resulting tridiagonal system can be solved for the YPP values.
  212. //
  213. // Modified:
  214. //
  215. // 07 January 2005 Shawn Freeman (pure C modifications)
  216. // 06 February 2004 John Burkardt
  217. //
  218. //
  219. // Author:
  220. //
  221. // John Burkardt
  222. //
  223. // Parameters:
  224. //
  225. // Input, int N, the number of data points. N must be at least 2.
  226. // In the special case where N = 2 and IBCBEG = IBCEND = 0, the
  227. // spline will actually be linear.
  228. //
  229. // Input, double T[N], the knot values, that is, the points were data is
  230. // specified. The knot values should be distinct, and increasing.
  231. //
  232. // Input, double Y[N], the data values to be interpolated.
  233. //
  234. // Input, int IBCBEG, left boundary condition flag:
  235. // 0: the cubic spline should be a quadratic over the first interval;
  236. // 1: the first derivative at the left endpoint should be YBCBEG;
  237. // 2: the second derivative at the left endpoint should be YBCBEG.
  238. //
  239. // Input, double YBCBEG, the values to be used in the boundary
  240. // conditions if IBCBEG is equal to 1 or 2.
  241. //
  242. // Input, int IBCEND, right boundary condition flag:
  243. // 0: the cubic spline should be a quadratic over the last interval;
  244. // 1: the first derivative at the right endpoint should be YBCEND;
  245. // 2: the second derivative at the right endpoint should be YBCEND.
  246. //
  247. // Input, double YBCEND, the values to be used in the boundary
  248. // conditions if IBCEND is equal to 1 or 2.
  249. //
  250. // Output, double SPLINE_CUBIC_SET[N], the second derivatives of the cubic spline.
  251. //
  252. static double *spline_cubic_set(int n, double t[], double y[], int ibcbeg,
  253. double ybcbeg, int ibcend, double ybcend)
  254. {
  255. double *a;
  256. double *b;
  257. int i;
  258. double *ypp;
  259. //
  260. // Check.
  261. //
  262. if (n <= 1) {
  263. nc_message(NC_SET_ERROR, "spline_cubic_set() error: "
  264. "The number of data points must be at least 2.\n");
  265. return NULL;
  266. }
  267.  
  268. for (i = 0; i < n - 1; i++) {
  269. if (t[i + 1] <= t[i]) {
  270. nc_message(NC_SET_ERROR, "spline_cubic_set() error: "
  271. "The knots must be strictly increasing, but "
  272. "T(%u) = %e, T(%u) = %e\n", i, t[i], i + 1, t[i + 1]);
  273. return NULL;
  274. }
  275. }
  276. a = (double *)calloc(3 * n, sizeof(double));
  277. nc_merror(a, "spline_cubic_set");
  278. b = (double *)calloc(n, sizeof(double));
  279. nc_merror(b, "spline_cubic_set");
  280. //
  281. // Set up the first equation.
  282. //
  283. if (ibcbeg == 0) {
  284. b[0] = 0.0E+00;
  285. a[1 + 0 * 3] = 1.0E+00;
  286. a[0 + 1 * 3] = -1.0E+00;
  287. } else if (ibcbeg == 1) {
  288. b[0] = (y[1] - y[0]) / (t[1] - t[0]) - ybcbeg;
  289. a[1 + 0 * 3] = (t[1] - t[0]) / 3.0E+00;
  290. a[0 + 1 * 3] = (t[1] - t[0]) / 6.0E+00;
  291. } else if (ibcbeg == 2) {
  292. b[0] = ybcbeg;
  293. a[1 + 0 * 3] = 1.0E+00;
  294. a[0 + 1 * 3] = 0.0E+00;
  295. } else {
  296. nc_message(NC_SET_ERROR, "spline_cubic_set() error: "
  297. "IBCBEG must be 0, 1 or 2. The input value is %u.\n", ibcbeg);
  298. free(a);
  299. free(b);
  300. return NULL;
  301. }
  302. //
  303. // Set up the intermediate equations.
  304. //
  305. for (i = 1; i < n - 1; i++) {
  306. b[i] = (y[i + 1] - y[i]) / (t[i + 1] - t[i])
  307. - (y[i] - y[i - 1]) / (t[i] - t[i - 1]);
  308. a[2 + (i - 1) * 3] = (t[i] - t[i - 1]) / 6.0E+00;
  309. a[1 + i * 3] = (t[i + 1] - t[i - 1]) / 3.0E+00;
  310. a[0 + (i + 1) * 3] = (t[i + 1] - t[i]) / 6.0E+00;
  311. }
  312. //
  313. // Set up the last equation.
  314. //
  315. if (ibcend == 0) {
  316. b[n - 1] = 0.0E+00;
  317. a[2 + (n - 2) * 3] = -1.0E+00;
  318. a[1 + (n - 1) * 3] = 1.0E+00;
  319. } else if (ibcend == 1) {
  320. b[n - 1] = ybcend - (y[n - 1] - y[n - 2]) / (t[n - 1] - t[n - 2]);
  321. a[2 + (n - 2) * 3] = (t[n - 1] - t[n - 2]) / 6.0E+00;
  322. a[1 + (n - 1) * 3] = (t[n - 1] - t[n - 2]) / 3.0E+00;
  323. } else if (ibcend == 2) {
  324. b[n - 1] = ybcend;
  325. a[2 + (n - 2) * 3] = 0.0E+00;
  326. a[1 + (n - 1) * 3] = 1.0E+00;
  327. } else {
  328. nc_message(NC_SET_ERROR, "spline_cubic_set() error: "
  329. "IBCEND must be 0, 1 or 2. The input value is %u", ibcend);
  330. free(a);
  331. free(b);
  332. return NULL;
  333. }
  334. //
  335. // Solve the linear system.
  336. //
  337. if (n == 2 && ibcbeg == 0 && ibcend == 0) {
  338. ypp = (double *)calloc(2, sizeof(double));
  339. nc_merror(ypp, "spline_cubic_set");
  340.  
  341. ypp[0] = 0.0E+00;
  342. ypp[1] = 0.0E+00;
  343. } else {
  344. ypp = d3_np_fs(n, a, b);
  345.  
  346. if (!ypp) {
  347. nc_message(NC_SET_ERROR, "spline_cubic_set() error: "
  348. "The linear system could not be solved.\n");
  349. free(a);
  350. free(b);
  351. return NULL;
  352. }
  353.  
  354. }
  355.  
  356. free(a);
  357. free(b);
  358. return ypp;
  359. }
  360.  
  361. //**********************************************************************
  362. //
  363. // Purpose:
  364. //
  365. // SPLINE_CUBIC_VAL evaluates a piecewise cubic spline at a point.
  366. //
  367. // Discussion:
  368. //
  369. // SPLINE_CUBIC_SET must have already been called to define the values of YPP.
  370. //
  371. // For any point T in the interval T(IVAL), T(IVAL+1), the form of
  372. // the spline is
  373. //
  374. // SPL(T) = A
  375. // + B * ( T - T(IVAL) )
  376. // + C * ( T - T(IVAL) )**2
  377. // + D * ( T - T(IVAL) )**3
  378. //
  379. // Here:
  380. // A = Y(IVAL)
  381. // B = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
  382. // - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
  383. // C = YPP(IVAL) / 2
  384. // D = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
  385. //
  386. // Modified:
  387. //
  388. // 07 January 2005 Shawn Freeman (pure C modifications)
  389. // 04 February 1999 John Burkardt
  390. //
  391. // Author:
  392. //
  393. // John Burkardt
  394. //
  395. // Parameters:
  396. //
  397. // Input, int n, the number of knots.
  398. //
  399. // Input, double T[N], the knot values.
  400. //
  401. // Input, double TVAL, a point, typically between T[0] and T[N-1], at
  402. // which the spline is to be evalulated. If TVAL lies outside
  403. // this range, extrapolation is used.
  404. //
  405. // Input, double Y[N], the data values at the knots.
  406. //
  407. // Input, double YPP[N], the second derivatives of the spline at
  408. // the knots.
  409. //
  410. // Output, double *YPVAL, the derivative of the spline at TVAL.
  411. //
  412. // Output, double *YPPVAL, the second derivative of the spline at TVAL.
  413. //
  414. // Output, double SPLINE_CUBIC_VAL, the value of the spline at TVAL.
  415. //
  416. static double spline_cubic_val(int n, double t[], double tval, double y[],
  417. double ypp[], double *ypval, double *yppval)
  418. {
  419. double dt;
  420. double h;
  421. int i;
  422. int ival;
  423. double yval;
  424. //
  425. // Determine the interval [ T(I), T(I+1) ] that contains TVAL.
  426. // Values below T[0] or above T[N-1] use extrapolation.
  427. //
  428. ival = n - 2;
  429.  
  430. for (i = 0; i < n - 1; i++) {
  431. if (tval < t[i + 1]) {
  432. ival = i;
  433. break;
  434. }
  435. }
  436. //
  437. // In the interval I, the polynomial is in terms of a normalized
  438. // coordinate between 0 and 1.
  439. //
  440. dt = tval - t[ival];63
  441. h = t[ival + 1] - t[ival];
  442.  
  443. yval = y[ival]
  444. + dt * ((y[ival + 1] - y[ival]) / h
  445. - (ypp[ival + 1] / 6.0E+00 + ypp[ival] / 3.0E+00) * h
  446. + dt * (0.5E+00 * ypp[ival]
  447. + dt * ((ypp[ival + 1] - ypp[ival]) / (6.0E+00 * h))));
  448.  
  449. *ypval = (y[ival + 1] - y[ival]) / h
  450. - (ypp[ival + 1] / 6.0E+00 + ypp[ival] / 3.0E+00) * h
  451. + dt * (ypp[ival]
  452. + dt * (0.5E+00 * (ypp[ival + 1] - ypp[ival]) / h));
  453.  
  454. *yppval = ypp[ival] + dt * (ypp[ival + 1] - ypp[ival]) / h;
  455.  
  456. return yval;
  457. }
  458.  
  459.  
  460. int main(int argc, char* argv[])
  461. {
  462.  
  463. double x[]={0,1,2,3,4}, y[]={0, 3, 3.5, 4, 4.05};
  464. double *setres = spline_cubic_set(5, x, y, 2, 0, 2, 0);
  465. int i; double tv;
  466.  
  467. double ypval, yppval;
  468. double valres, tval;
  469.  
  470.  
  471. for (i=0; i!=5; i++)
  472. {printf("AAAAA %lf ", setres[i]); printf("%lf\n", (spline_cubic_val(5, x, x[i], y, setres, &ypval, &yppval), ypval) ); }
  473.  
  474. printf("end of derivatives\n");
  475.  
  476. for (tv=0; tv!=4.125; tv+=0.125)
  477. {printf("%lf %lf\n", tv , spline_cubic_val(5, x, tv, y, setres, &ypval, &yppval)); }
  478.  
  479.  
  480.  
  481. while(scanf("%lf",&tval)==1)
  482. { valres = spline_cubic_val(5, x, tval, y, setres, &ypval, &yppval);
  483.  
  484. printf("%lf\n", valres);}
  485. scanf("\n");
  486.  
  487. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!