1 /*
2 Copyright (C) 2016 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira, X. Andrade
3
4 // Spline routines produced at the Lawrence Livermore National
5 // Laboratory. Written by Xavier Andrade (xavier@llnl.gov), Erik
6 // Draeger (draeger1@llnl.gov) and Francois Gygi (fgygi@ucdavis.edu).
7
8 This program is free software; you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation; either version 2, or (at your option)
11 any later version.
12
13 This program is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with this program; if not, write to the Free Software
20 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 02110-1301, USA.
22
23 */
24
25 #include <config.h>
26
27 #include <stdio.h>
28 #include <math.h>
29 #include <stdlib.h>
30
31 #include <gsl/gsl_spline.h>
32
33 #include "string_f.h"
34
35 #include <fortran_types.h>
36
37 #include <assert.h>
38 #include <iostream>
39 #include <vector>
40
41 /* Interpolation */
FC_FUNC_(oct_spline_end,OCT_SPLINE_END)42 extern "C" void FC_FUNC_(oct_spline_end, OCT_SPLINE_END)(void **spl, void **acc) {
43 gsl_spline_free((gsl_spline *)(*spl));
44 gsl_interp_accel_free((gsl_interp_accel *)(*acc));
45 }
46
FC_FUNC_(oct_spline_fit,OCT_SPLINE_FIT)47 extern "C" void FC_FUNC_(oct_spline_fit, OCT_SPLINE_FIT)
48 (const fint *nrc, const double *x, const double *y, void **spl, void **acc){
49
50 /* the GSL headers actually specify size_t instead of const int for nrc */
51 *acc = (void *)gsl_interp_accel_alloc();
52 *spl = (void *)gsl_spline_alloc(gsl_interp_cspline, *nrc);
53 gsl_spline_init((gsl_spline *)(*spl), x, y, *nrc);
54 fflush(stdout);
55 }
56
FC_FUNC_(oct_spline_eval,OCT_SPLINE_EVAL)57 extern "C" double FC_FUNC_(oct_spline_eval, OCT_SPLINE_EVAL)
58 (const double *x, const void **spl, void **acc){
59 /* the GSL headers specify double x instead of const double x */
60 return gsl_spline_eval((gsl_spline *)(*spl), *x, (gsl_interp_accel *)(*acc));
61 }
62
63
64 template <typename Type, int stride>
oct_spline_eval_array(fint nn,Type * xf,const void ** spl,void ** acc)65 void oct_spline_eval_array(fint nn, Type * xf, const void **spl, void **acc){
66 for(fint ii = 0; ii < nn; ii++){
67 xf[stride*ii] = Type(gsl_spline_eval((gsl_spline *)(*spl), xf[stride*ii], (gsl_interp_accel *)(*acc)));
68 }
69 }
70
FC_FUNC_(oct_spline_eval_array,OCT_SPLINE_EVAL_ARRAY)71 extern "C" void FC_FUNC_(oct_spline_eval_array, OCT_SPLINE_EVAL_ARRAY)
72 (const fint * nn, double *xf, const void **spl, void **acc){
73 oct_spline_eval_array<double, 1>(*nn, xf, spl, acc);
74 }
75
76 /* use a stride of 2 to store into just the real part of a Fortran complex array */
FC_FUNC_(oct_spline_eval_arrayz,OCT_SPLINE_EVAL_ARRAYZ)77 extern "C" void FC_FUNC_(oct_spline_eval_arrayz, OCT_SPLINE_EVAL_ARRAYZ)
78 (const fint * nn, double *xf, const void **spl, void **acc){
79 oct_spline_eval_array<double, 2>(*nn, xf, spl, acc);
80 }
81
82 /* This function returns the number of points with which a spline
83 was constructed (the size component of the gsl_spline struct). */
FC_FUNC_(oct_spline_npoints,OCT_SPLINE_NPOINTS)84 extern "C" fint FC_FUNC_(oct_spline_npoints, OCT_SPLINE_NPOINTS)(const void **spl, void **acc){
85 return (fint)((gsl_spline *)(*spl))->size;
86 }
87
88 /* This function places in the x array the x values of a given spline spl*/
FC_FUNC_(oct_spline_x,OCT_SPLINE_X)89 extern "C" void FC_FUNC_(oct_spline_x, OCT_SPLINE_X)(const void **spl, void **acc, double *x){
90
91 int size, i;
92 size = (int)((gsl_spline *)(*spl))->size;
93 for(i=0; i<size; i++) x[i] = ((gsl_spline *)(*spl))->x[i];
94 }
95
96 /* This function places in the y array the y values of a given spline spl*/
FC_FUNC_(oct_spline_y,OCT_SPLINE_Y)97 extern "C" void FC_FUNC_(oct_spline_y, OCT_SPLINE_Y)(const void **spl, void **acc, double *y){
98 int size, i;
99
100 size = (int)((gsl_spline *)(*spl))->size;
101 for(i=0; i<size; i++) y[i] = ((gsl_spline *)(*spl))->y[i];
102 }
103
104 /* Returns the integral of the spline stored in spl, between a and b */
FC_FUNC_(oct_spline_eval_integ,OCT_SPLINE_EVAL_INTEG)105 extern "C" double FC_FUNC_(oct_spline_eval_integ, OCT_SPLINE_EVAL_INTEG)
106 (const void **spl, const double *a, const double *b, void **acc){
107 /* the GSL headers specify double a, double b */
108 return gsl_spline_eval_integ((gsl_spline *)(*spl), *a, *b, (gsl_interp_accel *)(* acc));
109 }
110
FC_FUNC_(oct_spline_eval_integ_full,OCT_SPLINE_EVAL_INTEG_FULL)111 extern "C" double FC_FUNC_(oct_spline_eval_integ_full, OCT_SPLINE_EVAL_INTEG_FULL)
112 (const void **spl, void **acc) {
113 /* the GSL headers specify double a, double b */
114 const int size = (int)((gsl_spline *)(*spl))->size;
115 const double a = ((gsl_spline *)(*spl))->x[0];
116 const double b = ((gsl_spline *)(*spl))->x[size - 1];
117 return gsl_spline_eval_integ((gsl_spline *)(*spl), a, b, (gsl_interp_accel *)(* acc));
118 }
119
120 /* Performs the derivative of a spline */
FC_FUNC_(oct_spline_eval_der,OCT_SPLINE_EVAL_DER)121 extern "C" double FC_FUNC_(oct_spline_eval_der, OCT_SPLINE_EVAL_DER)
122 (const double *x, const void **spl, void **acc){
123 /* the GSL headers specify double x */
124 return gsl_spline_eval_deriv((gsl_spline *)(*spl), *x, (gsl_interp_accel *)(*acc));
125 }
126
127 /* Performs the second derivative of a spline */
FC_FUNC_(oct_spline_eval_der2,OCT_SPLINE_EVAL_DER2)128 extern "C" double FC_FUNC_(oct_spline_eval_der2, OCT_SPLINE_EVAL_DER2)
129 (const double *x, const void **spl, void **acc){
130 /* the GSL headers specify double x */
131 return gsl_spline_eval_deriv2((gsl_spline *)(*spl), *x, (gsl_interp_accel *)(*acc));
132 }
133