home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / lsqrft15.zip / funclib2.c < prev    next >
Text File  |  1994-01-05  |  13KB  |  517 lines
  1. #include <math.h>
  2. #include <string.h>
  3. #include "fit.h"
  4. #define NUM_FCNS 20   /* CHANGE if you define more than 20 functions */
  5.  
  6. extern int debug;
  7. /* a few function declarations */
  8. int fgauss();
  9. int fgaussc();
  10. int fgaussn();
  11. int florenz();
  12. int florenz2();
  13. int fline();
  14. int fpoly();
  15. int fexpn();
  16. int fxyquad();
  17. int fxygauss();
  18. int fconic();
  19. int fsincos();
  20.  
  21. /* The structure fcn holds the information about each function */
  22. /* Specifically, it is here that the name used by the command line, */
  23. /* the name of the function in c, the number of parameters, */
  24. /* and a comment are defined.  If gnuplot is used for plotting, */
  25. /* then the comment must be a valid definition of the function */
  26. /* on the gnuplot command line. */
  27. /* to add a function, you must add a line to the fcn initialization */
  28. /* For a function with a variable number of parameters */
  29. /* you must initialize na to be -1.  You have to choose the number */
  30. /* of parameters when you choose the function at the fit> command */
  31. /* prompt. */
  32. /* If a comment begins with "FILE", plotting is done through a */
  33. /* tempory file.  For functions with variable number of parameters, */
  34. /* you must begin your comment with "FILE" */
  35. /* Functions of two independent variables which are not plotted */
  36. /* through a file should be expressed in terms of u and v instead */
  37. /* of x and y.  This is because gnuplot wants it that way. */
  38.  
  39. /* The elements of the struct routine are as follows: */
  40. /* char name[20]    the name of the function at the fit> command line */
  41. /* int num_indep   the number of independent variables */
  42. /* int linflag   Is the function linear in the parameters? */
  43. /* linflag = 0 for a non-linear function.  linflag = 1 for */
  44. /* a function linear in the parameters */
  45. /* int (*func)()   the name of the function as known to the */
  46. /*                         compiler */
  47. /* int na        the number of parameters, -1 for variable number */
  48. /* char comment[160] a comment about the function */
  49.  
  50. struct routine{
  51.     char name[20];
  52.     int num_indep;
  53.     int linflag;
  54.     int (*func)(double *x, double *a, double *y, double *dyda,
  55.            int na, int dyda_flag, int *fita, int dydx_flag, 
  56.            double *dydx, double ydat);
  57.     int na;
  58.     char comment[160];
  59. } fcn[NUM_FCNS]={
  60.  
  61.    {"gauss",1,0, &fgauss, 3, "f(x) = a2*exp(-((x-a0)/a1)**2)" },
  62.    {"gaussc",1,0, &fgaussc, 4,"f(x) = a2*exp(-((x-a0)/a1)**2) + a3" },
  63.    {"ngauss",1,0, &fgaussn, -1,"FILE f(x) = sum( a[i+2]*exp(-((x-ai)/a[i+1])**2) ) + a[n-1]"},
  64.    {"lorenz", 1,0, &florenz,3,"f(x) = a1*a2/(4*(x-a0)**2 + a1)" },
  65.    {"2lorenz",1,0,&florenz2,7,"f(x) = a1*a2/(4*(x-a0)**2 + a1) + a4*a5/(4*(x-a3)**2 + a4) + a6" },
  66.    {"line",1,1,&fline, 2, "f(x) = a0 + a1*x"},
  67.    {"poly", 1,1,&fpoly, -1, "FILE f(x) = sum( ai * pow(x,i) )"},
  68.    {"nexp", 1,0,&fexpn, -1, "FILE f(x) = sum( a[i+2]*exp((x-a[i])/a[i+1]) ) + a[n-1]"},
  69.    {"xyquad", 2,1,&fxyquad, 6, "f(u,v) = a0 + a1*u + a2*v + a3*u*u + a4*v*v + a5*u*v"},
  70.    {"xygauss",2,0,&fxygauss, 6, "f(u,v) = a4*exp(-((u-a0)/a2)**2 - ((v-a1)/a3)**2 ) + a5"},
  71.    {"conic",1,0,&fconic, 6, "FILE a0*x*x + a1*x*y + a2*y*y + a3*x + a4*y +a5 = 0"},
  72.    {"sincos",1,0,&fsincos, -1, "FILE a0 + a[i]*sin(a[i+1]*x) + a[i+2]*cos(a[i+3]*x)"},
  73. };
  74.  
  75. /* The function getfcnptr() looks through the struct fcn to find a */
  76. /* pointer corresponding to the function specified at the fit> */
  77. /* prompt with the fn command */
  78.  
  79. int (*getfcnptr(char *name, int *num_indep, int *linflag, int *na, char *comment))
  80. {
  81. int i=0;
  82. while(fcn[i].na != 0 ){
  83.     if(strcmp(name, fcn[i].name) == 0){
  84.         *na = fcn[i].na;
  85.         *linflag = fcn[i].linflag;
  86.         *num_indep = fcn[i].num_indep;
  87.         strcpy(comment,fcn[i].comment);
  88.         return fcn[i].func;
  89.     }
  90.     i++;
  91. }
  92.  
  93. /* If function name not found return pointer to NULL */
  94. *na = 0;
  95. return (int *)NULL;
  96. }
  97.  
  98. /* listfcns() lists the functions available */
  99. int listfcns(){
  100. int i=0 ;
  101. while(fcn[i].na != 0 ){
  102.     printf("%s %s\n", fcn[i].name, fcn[i].comment);
  103.     i++;
  104.     if(i == 20){
  105.         i = 0;
  106.         gets();
  107.     }
  108. }
  109. return 0;
  110. }
  111.  
  112. /* the definition of a function must be of the form: */
  113. /* int function(double *x, double *a, double *y, double *dyda, */
  114. /*                      int na, int dyda_flag, int *fita, int dydx_flag, */
  115. /*                      double *dydx, double ydat); */
  116.  
  117. /*
  118. /* The x[i]'s are the independent variables passed to the function */
  119. /* *y should equal the value of the function at x. */
  120. /* The a[i]'s are the parameters. */
  121. /* the dyda[i]'s should be set equal to the first derivative */
  122. /* of the fitting function with respect to each parameter */
  123. /* dyda_flag, fita[], and dydx_flag[] are flags which are */
  124. /* passed to the function and used for optimization */
  125. /* purposes only.  These need not be used at all, */
  126. /* but may result in a significant performance boost. */
  127. /* if calculating the derivative is slow and it is not */
  128. /* needed by the calling function. */
  129. /* dydx[i] should be set equal to the derivative of the fitting */
  130. /* function with respect to x[i].  It is used only if */
  131. /* the fitting algorithm is told to consider errors in the x[i] */
  132. /* ydat is the data value that you are fitting to. */
  133. /* It might be useful in multivalued functions */
  134. /* Using it might destroy the statistical relevance of your fit. */
  135.  
  136. /* The linear fitting routine makes different use of user defined */
  137. /* fitting functions.  the dyda[i]'s and dydx[i]'s are not needed. */
  138. /* Linear fitting functions which are set up properly for non-linear fitting */
  139. /* work fine with the linear fitting routine.  Using dydx_flag, */
  140. /* dyda_flag, and fita[i] will speed up linear fitting. */
  141.  
  142. /* a conic section, doesn't work very well */
  143. int fconic(x, a,y,dyda,na,dyda_flag, fita, dydx_flag, dydx, ydat)
  144. double x[];
  145. double a[],*y,dyda[];
  146. int na;
  147. int dyda_flag;
  148. int fita[];
  149. int dydx_flag[];
  150. double dydx[], ydat;
  151. {
  152. double A, B, C, D, E, F, Q, R, S;
  153. double yp, ym, x2;
  154. int pm;
  155.  
  156. A = a[0];
  157. B = a[1];
  158. C = a[2];
  159. D = a[3];
  160. E = a[4];
  161. F = a[5];
  162.  
  163. x2 = x[0]*x[0];
  164. R = B*x[0] + E;
  165. S = A*x2 + D*x[0] + F;
  166. if( (R*R - 4*C*S) < 0 ){
  167.     printf("Function conic failed, wanted to take sqrt of negative #.\n");
  168.     return 1;
  169. }
  170. Q = sqrt(R*R - 4*C*S);
  171.  
  172. yp = (-R + Q)/2*C;
  173. ym = (-R - Q)/2*C;
  174.  
  175. if( fabs(yp - ydat) < fabs(ym - ydat) ){
  176.     pm = +1;
  177.     *y = yp;
  178. }
  179. else{
  180.     pm = -1;
  181.     *y = ym;
  182. }
  183.  
  184. if(dyda_flag){
  185.     if(fita[0]) dyda[0] = -pm*x2/Q;
  186.     if(fita[1]) dyda[1] = -x[0]/(2*C) + pm*R*x[0]/(2*Q*C);
  187.     if(fita[2]) dyda[2] = (R/(2*C) + pm*S/Q)/C;
  188.     if(fita[3]) dyda[3] = -pm*x[0]/Q;
  189.     if(fita[4]) dyda[4] = -1/(2*C) + pm/(2*Q*C);
  190.     if(fita[5]) dyda[5] = -pm/Q;
  191. }
  192. if(dydx_flag[0])
  193.     dydx[0] = -B/(2*C) + pm*(2*B*(B*x[0]+E) - 
  194.                                     4*C*(2*x[0]*A + D))/(2*C*Q);
  195.  
  196. return 0;
  197. }
  198.  
  199. /*f(u,v) = a4*exp(-((u-a0)/a2)**2 - ((v-a1)/a3)**2 ) + a5 */
  200. /* a two dimensional gaussian */
  201. int fxygauss(x, a,y,dyda,na,dyda_flag, fita, dydx_flag, dydx, ydat)
  202. double x[];
  203. double a[],*y,dyda[];
  204. int na;
  205. int dyda_flag;
  206. int fita[];
  207. int dydx_flag[];
  208. double dydx[], ydat;
  209. {
  210. double fac0, ex0, arg0,fac1, ex1, arg1;
  211. arg0 = (x[0]-a[0])/a[2];
  212. ex0 = exp(-arg0*arg0);
  213. fac0 = a[4]*ex0*2.0*arg0;
  214. arg1 = (x[1]-a[1])/a[3];
  215. ex1 = exp(-arg1*arg1);
  216. fac1 = a[4]*ex1*2.0*arg1;
  217.  
  218. *y = a[4]*ex0*ex1 + a[5];
  219.  
  220. if(dyda_flag){
  221.     if(fita[0]) dyda[0] = ex1*fac0/a[2];
  222.     if(fita[1]) dyda[1] = ex0*fac1/a[3];
  223.     if(fita[2]) dyda[2] = ex1*fac0*arg0/a[2];
  224.     if(fita[3]) dyda[3] = ex0*fac1*arg1/a[3];
  225.     if(fita[4]) dyda[4] = ex0*ex1;
  226.     if(fita[5]) dyda[5] = 1;
  227. }
  228. if(dydx_flag[0]) dydx[0] = -ex1*fac0/a[2];
  229. if(dydx_flag[1]) dydx[1] = -ex0*fac1/a[3];
  230. return 0;
  231. }
  232.  
  233.  
  234. /* a sum of sines and cosines */
  235. int fsincos(x, a,y,dyda,na,dyda_flag, fita, dydx_flag, dydx, ydat)
  236. double x[];
  237. double a[],*y,dyda[];
  238. int na;
  239. int dyda_flag;
  240. int fita[];
  241. int dydx_flag[];
  242. double dydx[], ydat;
  243. {
  244.  
  245. int i;
  246. double temp;
  247.  
  248. *y = a[0];
  249. dydx[0] = 0;
  250.  
  251. dyda[0] = 1;
  252. for(i = 1; i < na-3; i+=4){
  253.     dyda[i] = sin(a[i+1]*x[0]);
  254.     *y += a[i]*dyda[i];
  255.     if(dydx_flag[0] || fita[i+1]) temp = a[i]*cos(a[i+1]*x[0]);
  256.     dydx[0] += a[i+1]*temp;
  257.     dyda[i+1] = x[0]*temp;
  258. }
  259.  
  260. for(i = 3; i < na-1; i+=4){
  261.     dyda[i] = cos(a[i+1]*x[0]);
  262.     *y += a[i]*dyda[i];
  263.     if(dydx_flag[0] || fita[i+1]) temp = -a[i]*sin(a[i+1]*x[0]);
  264.     dydx[0] += a[i+1]*temp;
  265.     dyda[i+1] = x[0]*temp;
  266. }
  267.  
  268. return 0;
  269. }
  270.  
  271.  
  272.  
  273. /* a quadradic in x and y */
  274. int fxyquad(x, a,y,dyda,na,dyda_flag, fita, dydx_flag, dydx, ydat)
  275. double x[];
  276. double a[],*y,dyda[];
  277. int na;
  278. int dyda_flag;
  279. int fita[];
  280. int dydx_flag[];
  281. double dydx[], ydat;
  282. {
  283. double x0, x1, x02, x12;
  284.  
  285. x0 = x[0];
  286. x1 = x[1];
  287. x02 = x0*x0;
  288. x12 = x1*x1;
  289.  
  290. *y = a[0] + a[1]*x0 + a[2]*x1 + a[3]*x02 + a[4]*x12 + a[5]*x1*x0;
  291. if(dyda_flag){
  292.     if(fita[0]) dyda[0] = 1;
  293.     if(fita[1]) dyda[1] = x0;
  294.     if(fita[2]) dyda[2] = x1;
  295.     if(fita[3]) dyda[3] = x02;
  296.     if(fita[4]) dyda[4] = x12;
  297.     if(fita[5]) dyda[5] = x1*x0;
  298. }
  299. if(dydx_flag[0]) dydx[0] = a[1] + 2*x0*a[3] + x1*a[5];
  300. if(dydx_flag[1]) dydx[1] = a[2] + 2*x1*a[4] + x0*a[5];
  301. return 0;
  302. }
  303.  
  304. /* a single lorenzian */
  305. int florenz(x, a,y,dyda,na,dyda_flag, fita, dydx_flag, dydx, ydat)
  306. double x[];
  307. double a[],*y,dyda[];
  308. int na;
  309. int dyda_flag;
  310. int fita[];
  311. int dydx_flag[];
  312. double dydx[], ydat;
  313. {
  314. double denom, del, del2;
  315.  
  316. del = x[0] - a[0];
  317. del2 = del * del;
  318. denom = 4.0*del2 + a[1];
  319. *y = a[2]*a[1]/denom;
  320. if(dyda_flag){
  321.     if(fita[0]) dyda[0] = 8.0*del*(*y)/denom;
  322.     if(fita[1]) dyda[1] = (*y)/a[1] - (*y)/denom;
  323.     if(fita[2]) dyda[2] = (*y)/a[2];
  324. }
  325. if(dydx_flag[0]) dydx[0] = -8.0*del*(*y)/denom;
  326. return 0;
  327. }
  328.  
  329. /* sum of two lorenzians */
  330. int florenz2(x,a,y,dyda,na,dyda_flag, fita, dydx_flag, dydx, ydat)
  331. double x[],a[],*y,dyda[];
  332. int na;
  333. int dyda_flag;
  334. int fita[];
  335. int dydx_flag[];
  336. double dydx[], ydat;
  337. {
  338. double denom, del, del2,y1, y2;
  339.  
  340. del = x[0] - a[0];
  341. del2 = del * del;
  342. denom = 4.0*del2 + a[1];
  343.  y1 = a[2]*a[1]/denom;
  344. if(dyda_flag){
  345.     if(fita[0]) dyda[0] = 8.0*del*y1/denom;
  346.     if(fita[1]) dyda[1] = y1/a[1] - y1/denom;
  347.     if(fita[1]) dyda[2] = y1/a[2];
  348. }
  349. if(dydx_flag[0]) dydx[0] =  -8.0*del*y1/denom;
  350.  
  351. del = x[0] - a[3];
  352. del2 = del * del;
  353. denom = 4.0*del2 + a[4];
  354. y2 = a[4]*a[5]/denom;
  355. if(dyda_flag){
  356.     if(fita[3]) dyda[3] = 8.0*del*y2/denom;
  357.     if(fita[4]) dyda[4] = y2/a[4] - y2/denom;
  358.     if(fita[5]) dyda[5] = y2/a[5];
  359.     if(fita[6]) dyda[6] = 1;
  360. }
  361. if(dydx_flag[0]) dydx[0] +=  -8.0*del*y2/denom;
  362. *y = y1+y2 + a[6];
  363. return 0;
  364. }
  365.  
  366. /* a gaussian */
  367. int fgauss(x,a,y,dyda,na,dyda_flag, fita, dydx_flag, dydx, ydat)
  368. double x[],a[],*y,dyda[];
  369. int na;
  370. int dyda_flag;
  371. int fita[];
  372. int dydx_flag[];
  373. double dydx[], ydat;
  374. {
  375. double fac, ex, arg;
  376. arg = (x[0]-a[0])/a[1];
  377. ex = exp (-arg*arg);
  378. fac = a[2]*ex*2.0*arg;
  379. *y = a[2]*ex;
  380.  
  381. if(dyda_flag){
  382.     if(fita[0]) dyda[0] = fac/a[1];
  383.     if(fita[1]) dyda[1] = fac*arg/a[1];
  384.     if(fita[2]) dyda[2] = ex;
  385. }
  386. if(dydx_flag[0]) dydx[0] = -fac/a[1];
  387. return 0;
  388. }
  389.  
  390. /* a gaussian plus a constant */
  391. int fgaussc(x,a,y,dyda,na,dyda_flag, fita, dydx_flag, dydx, ydat)
  392. double x[],a[],*y,dyda[];
  393. int na;
  394. int dyda_flag;
  395. int fita[];
  396. int dydx_flag[];
  397. double dydx[], ydat;
  398. {
  399. double fac, ex, arg;
  400. arg = (x[0]-a[0])/a[1];
  401. ex = exp (-arg*arg);
  402. fac = a[2]*ex*2.0*arg;
  403. *y = a[2]*ex + a[3];
  404. if(dyda_flag){
  405.     if(fita[0]) dyda[0] = fac/a[1];
  406.     if(fita[1]) dyda[1] = fac*arg/a[1];
  407.     if(fita[2]) dyda[2] = ex;
  408.     if(fita[3]) dyda[3] = 1;
  409. }
  410. if(dydx_flag[0]) dydx[0] = -fac/a[1];
  411. return 0;
  412. }
  413.  
  414. /* sum of gaussians plus a constant */
  415. int fgaussn(x,a,y,dyda,na,dyda_flag, fita, dydx_flag, dydx, ydat)
  416. double x[],a[],*y,dyda[];
  417. int na;
  418. int dyda_flag;
  419. int fita[];
  420. int dydx_flag[];
  421. double dydx[], ydat;
  422. {
  423. double fac, ex, arg;
  424. int i;
  425. dydx[0] = 0;
  426. *y = 0;
  427. for(i = 0; i < na-1; i+=3){
  428.     arg = (x[0]-a[i])/a[i+1];
  429.     ex = exp (-arg*arg);
  430.     fac = a[i+2]*ex*2.0*arg;
  431.     *y += a[i+2]*ex;
  432.     if(dyda_flag){
  433.         if(fita[i]) dyda[i] = fac/a[i+1];
  434.         if(fita[i+1]) dyda[i+1] = fac*arg/a[i+1];
  435.         if(fita[i+2]) dyda[i+2] = ex;
  436.     }
  437.     if(dydx_flag[0]) dydx[0] += -fac/a[i+1];
  438. }
  439. dyda[na-1] = 1;
  440. *y += a[na-1];
  441. return 0;
  442. }
  443.  
  444.  
  445. int fline(x,a,y,dyda,na,dyda_flag, fita, dydx_flag, dydx, ydat)
  446. double x[],a[],*y,dyda[];
  447. int na;
  448. int dyda_flag;
  449. int fita[];
  450. int dydx_flag[];
  451. double dydx[], ydat;
  452. {
  453. dyda[0] = 1;
  454. dyda[1] = x[0];
  455. *y = a[0] + a[1]*x[0];
  456. dydx[0] = a[1];
  457. return 0;
  458. }
  459.  
  460. int fpoly(x,a,y,dyda,na,dyda_flag, fita, dydx_flag, dydx, ydat)
  461. double x[],a[],*y,dyda[];
  462. int na;
  463. int dyda_flag;
  464. int fita[];
  465. int dydx_flag[];
  466. double dydx[], ydat;
  467. {
  468. int i;
  469. double xn, ynum;
  470. xn = 1;
  471. ynum = 0;
  472. dydx[0] = 0;
  473. for( i = 0; i < na; i++){
  474.     ynum += xn*a[i];
  475.     if(dydx_flag[0] && i > 0) dydx[0] += a[i-1] * xn;
  476.     if(fita[i]) dyda[i] = xn;
  477.     xn *= x[0];
  478. }
  479. *y = ynum;
  480. return 0;
  481. }
  482.  
  483. /* sum of exponentials plus a constant */
  484. int fexpn(x,a,y,dyda,na,dyda_flag, fita, dydx_flag, dydx, ydat)
  485. double x[],a[],*y,dyda[];
  486. int na;
  487. int dyda_flag;
  488. int fita[];
  489. int dydx_flag[];
  490. double dydx[], ydat;
  491. {
  492. double fac, ex, arg;
  493. int i;
  494. double s;
  495.  
  496. *y = 0;
  497. dydx[0] = 0;
  498. for(i = 0; i < na-1; i+=3){
  499.      if(x[0] < a[i]) s = 0;
  500.         else s = 1;
  501.     arg = (x[0]-a[i])/a[i+1];
  502.     ex = exp (-arg);
  503.     fac = a[i+2]*ex;
  504.     *y += s*fac;
  505.     if(dyda_flag){
  506.     if(fita[i]) dyda[i] = s*fac/a[i+1];
  507.     if(fita[i+1]) dyda[i+1] = s*fac*arg/a[i+1];
  508.     if(fita[i+2]) dyda[i+2] = s*ex;
  509.     }
  510.     if(dydx_flag[0]) dydx[0] += -s*fac/a[i+1];
  511. }
  512. dyda[na-1] = 1;
  513. *y += a[na-1];
  514. return 0;
  515. }
  516.  
  517.