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