home *** CD-ROM | disk | FTP | other *** search
/ Between Heaven & Hell 2 / BetweenHeavenHell.cdr / 300 / 249 / spline.arc / SPLINE.C < prev    next >
C/C++ Source or Header  |  1987-06-16  |  17KB  |  600 lines
  1. /*    spline - interpolate points in a tabulated function
  2.  
  3.     Usage...
  4.         spline [file][options]
  5.         options are:
  6.           -a  [step [start]] automatic abscissas 
  7.           -b                 break interpolation after each label 
  8.           -c                 general curve 
  9.           -i  file           interpolate at x values given in file.  Only
  10.                                1st number on each line of file is used -
  11.                                rest of line is ignored.  Blank lines and
  12.                                lines starting with ';' are ignored.
  13.           -n  num            interpolate over num intervals 
  14.           -q                 increase intervals fourfold 
  15.           -t  num            tension in interpolating curve 
  16.           -x  [min [max]]    interpolate from min to max 
  17.           -xl                take logs of x values before interpolating 
  18.           -yl                take logs of y values before interpolating 
  19.           -zl                take logs of z values before interpolating 
  20.           -3                 3D function or curve
  21.  
  22.     Notes...
  23.         This program fits a smooth curve through a given set of points,
  24.         using the splines under tension introduced by Schweikert [J. 
  25.         Math.  and Physics v 45 pp 312-317 (1966)].  They are similar
  26.         to cubic splines, but are less likely to introduce spurious
  27.         inflection points.
  28.  
  29.     History...
  30.         15 Jun 87    echoes command line to output file.
  31.         -- ver 2.1 --
  32.         25 Nov 86    3D
  33.         -- ver 2.0 --
  34.         16 Oct 86    accepts desired abscissas from file.
  35.         -- ver 1.4 --
  36.         7 Jun 86    echoes a 2-point line.
  37.         3 Jun 86    -s switch allows specified slopes at ends.
  38.         2 Jun 86    Ver 1.3: Fixed several embarrassing bugs in -b option.
  39.                     -q switch quadruples number of intervals.
  40.         6 Apr 86    -b switch breaks interpolation at labels.
  41.         21 Sep 85    Added -xl and -yl switches for taking logs
  42.         23 Nov 85    Passing lines starting with ';' unchanged, otherwise
  43.                     ignoring them.
  44.     Bugs...
  45.         Should have -e command for evenly spaced abscissas.
  46.         -q switch should accept argument for desired fractional
  47.         increase in # points.
  48.         Should accept the output file as an optional second argument.
  49.  
  50.     Author...
  51.         Copyright (c) 1985, 1986  James R. Van Zandt (jrv@mitre-bedford)
  52.  
  53.         All rights reserved.
  54.         This program may be copied for personal, non-profit use only.
  55.  
  56.         Based on algorithms by A.  K.  Cline [Communications of the ACM
  57.         v 17 n 4 pp 218-223 (Apr 74)].
  58.  
  59. */
  60.  
  61. #include <stdio.h>
  62. #include <math.h>
  63.  
  64. #define VERSION "2.1"
  65. char *RIGHT=        "Copyright (c) 1985, 1986  James R. Van Zandt";
  66.  
  67. #define DESMET        /* this enables echoing of command line to output */
  68.  
  69. #define ENTRIES 200
  70. #define MAXLABELS 50
  71. #define BUFSIZE 200
  72.  
  73. double *x, *y, *z, *temp, *xp, *yp, *zp, *path;
  74. double slp1=0.,slpn=0.,sigma=1.;
  75. double abscissa=0.;
  76. double abscissa_step=1.;
  77. double slope1=0.,slopen=0.;    /* slopes supplied by user */
  78. double xmin=0.;
  79. double xmax=0.;
  80. double curv2(), mylog();
  81.  
  82. int abscissa_arguments=0;
  83. int automatic_abscissas=0;
  84. int breaking=0;        /* nonzero if breaking at labels */
  85. int curve=0;        /* nonzero if general curve permitted */
  86. int debugging=0;
  87. int given_abscissas=0;    /* nonzero if abscissas are given in a file */
  88. int labels=0;        /* number of labels in input data */
  89. int quadruple=0;    /* nonzero if user asking for 4 times the resolution */
  90. int slopes=0;        /* nonzero if slopes supplied by user */
  91. int threed=0;        /* nonzero if 3D curve (x, y, and z given for each point) */
  92. int xlog=0;        /* nonzero if taking log of x values */
  93. int ylog=0;        /* nonzero if taking log of y values */
  94. int zlog=0;        /* nonzero if taking log of z values */
  95.  
  96. int x_arguments=0;
  97. int n;                /* number of entries in x, y, yp, and temp */
  98. int index_array[MAXLABELS];    /* indexes into x and y */
  99. int *p_data=index_array;
  100. int total=100;
  101.  
  102. FILE *file, *afile;
  103.  
  104. char buf[BUFSIZE];
  105. char *label_array[MAXLABELS];
  106. char **p_text=label_array;
  107. char newline[]="\n", nothing[]="", *leadin;
  108.  
  109. main(argc,argv) int argc; char **argv;
  110. {    int nn,origin;
  111.     read_data(argc,argv);
  112.     if(breaking && labels)
  113.         {origin=0;
  114.         while(labels--)
  115.             {p_data[0] -= origin;
  116.             nn=p_data[0]+1;
  117.             if(quadruple) total=nn*4-3;
  118.             if(nn) curv0(nn,total);
  119.             origin += nn;
  120.             path += nn;
  121.             x += nn;
  122.             y += nn;
  123.             p_data++;
  124.             p_text++;
  125.             }
  126.         }
  127.     else
  128.         {if(quadruple) total=n*4-3;
  129.         curv0(n,total);
  130.         }
  131. }
  132.  
  133. curv0(n,requested)
  134. int n,requested;
  135. {    int i, j, each, stop, seg=0, regressing=0;
  136.     double s, ds, xx, yy, zz;
  137.     char *t;
  138.  
  139.     if(n<3)
  140.         {for (j=0; j<n; j++)
  141.             {xx=x[n-1]; if(xlog) xx=exp(xx);
  142.             yy=y[n-1]; if(ylog) yy=exp(yy);
  143.             printf("%g %g ",xx,yy);
  144.             if(threed)
  145.                 {zz=z[n-1]; if(zlog) zz=exp(zz);
  146.                 printf("%g ",zz);
  147.                 }
  148.             if(j==p_data[seg]) puts(p_text[seg++]);
  149.             putchar('\n');
  150.             }
  151.         return;
  152.         }
  153.     for(i=1; i<n; i++) if(path[i]<=path[i-1]) regressing++;
  154.     if (regressing)
  155.         {fprintf(stderr,"ERROR: independent variable not strictly increasing\n");
  156.         return;
  157.         }
  158.     if (curve && (xmin<0. || xmax>1.))
  159.         {fprintf(stderr,"ERROR: independent variable not in range 0 to 1\n");
  160.         return;
  161.         }
  162.     if(curve)
  163.         {curv1(path,x,xp,n); curv1(path,y,yp,n);
  164.         if(threed) curv1(path,z,zp,n);
  165.         }
  166.     else
  167.         {slp1=slope1; slpn=slopen; curv1(x,y,yp,n);
  168.         if(threed) 
  169.             {slp1=slpn=0.; curv1(x,z,zp,n);
  170.             }
  171.         }
  172.     s=path[0];
  173.     seg=0;
  174.     if(requested<n) requested=n;
  175.     if(given_abscissas)                /* abscissas specified in a file */
  176.         {leadin=nothing;
  177.         while(fgets(buf,BUFSIZE,afile))
  178.             {if(buf[0]==';') continue;  /* skip comment */
  179.             t=buf;
  180.             while(*t && isspace(*t)) t++;
  181.             if(*t == 0) continue;        /* skip blank lines */
  182.             buf[strlen(buf)-1]=0;        /* zero out the line feed */
  183.             sscanf(buf,"%F",&s);
  184.             if(xlog) s=mylog(s);
  185.             if(s<path[0] || path[n-1]<s) continue; /* ignore outlying points */
  186.             value(s, &xx, &yy, &zz, n);
  187.             if(threed) printf("%s%g %g %g", leadin, xx, yy, zz);
  188.             else printf("%s%g %g", leadin, xx, yy);
  189.             leadin=newline;
  190.             }
  191.         if(breaking && labels) printf(" %s\n",p_text[0]);
  192.         else printf(" \"\"\n");
  193.         rewind(afile);
  194.         }
  195.     else if(x_arguments==2)            /* specific range was requested */
  196.         {ds=(xmax-xmin)/requested;
  197.         if(curve)
  198.             {ds*=(path[n-1]-path[0]);
  199.             s=xmin*(path[n-1]-path[0])+path[0];
  200.             }
  201.         else s=xmin;
  202.         for(i=requested+1; i; i--)
  203.             {value(s,&xx,&yy,&zz,n);
  204.             if(threed) printf("%g %g %g ", xx, yy, zz);
  205.             else printf("%g %g ", xx, yy);
  206.             if(j==p_data[seg]) puts(p_text[seg++]);
  207.             putchar('\n');
  208.             s+=ds;
  209.             }
  210.         }
  211.     else                        /* spline for entire curve was requested */
  212.         {for (i=j=0; j<n-1; j++)
  213.             {stop=requested*(j+1)/(n-1);
  214.             ds=(path[j+1]-path[j])/(stop-i);
  215.             for ( ; i<stop; i++)
  216.                 {value(s,&xx,&yy,&zz,n);
  217.                 if(threed) printf("%g %g %g ", xx, yy, zz);
  218.                 else printf("%g %g ",xx,yy);
  219.                 if(j==p_data[seg]) puts(p_text[seg++]);
  220.                 putchar('\n');
  221.                 s+=ds;
  222.                 }
  223.             }
  224.         xx=x[n-1]; if(xlog) xx=exp(xx);
  225.         yy=y[n-1]; if(ylog) yy=exp(yy);
  226.         if(threed)
  227.             {zz=z[n-1]; if(zlog) zz=exp(zz);
  228.             printf("%g %g %g ", xx, yy, zz);
  229.             }
  230.         else printf("%g %g ",xx,yy);
  231.         if(j==p_data[seg]) puts(p_text[seg++]);
  232.         putchar('\n');
  233.         }
  234. }
  235.  
  236. value(s,q,r,v,n) double s,*q,*r,*v; int n;
  237. {    if(curve)
  238.         {*q=curv2(path,x,xp,s,n);
  239.         *r=curv2(path,y,yp,s,n);
  240.         if(threed) *v=curv2(path,z,zp,s,n);
  241.         }
  242.     else
  243.         {*q=s;
  244.         *r=curv2(x,y,yp,s,n);
  245.         if(threed) *v=curv2(x,z,zp,s,n);
  246.         }
  247.     if(xlog) *q=exp(*q);
  248.     if(ylog) *r=exp(*r);
  249.     if(zlog) *v=exp(*v);
  250. }
  251.  
  252. read_data(argc,argv) int argc; char **argv;
  253. {    int i, j, length, skipping;
  254.     double xx, yy, zz, ss, d, *pd, sum;
  255.     char *s,*t;
  256.  
  257. #ifdef DESMET
  258.     fprint_cmd(stdout, "; spline %s\n");
  259. #endif
  260.  
  261.     x=path=malloc(ENTRIES*sizeof(double));
  262.     y=malloc(ENTRIES*sizeof(double));
  263.     if(x==0 || y==0) {fprintf(stderr,"can\'t allocate buffer"); exit();}
  264.     if(argc>1 && streq(argv[1],"?")) help();
  265.     if(argc<=1 || *argv[1]=='-') file=stdin;
  266.     else
  267.         {if(argc>1)
  268.             {file=fopen(argv[1],"r");
  269.             if(file==0) {printf("file %s not found\n",argv[1]); exit();}
  270.             argc--; argv++;
  271.             }
  272.         else help();
  273.         }
  274.     argc--; argv++;
  275.     while(argc>0)
  276.         {i=get_parameter(argc,argv);
  277.         argv=argv+i; argc=argc-i;
  278.         }
  279.     if(threed)
  280.         {z=malloc(ENTRIES*sizeof(double));
  281.         if(z==0) {fprintf(stderr,"can\'t allocate buffer"); exit();}
  282.         }
  283.     if(sigma<=0.)
  284.         {fprintf(stderr,"ERROR: tension must be positive\n");
  285.         exit(1);
  286.         }
  287.     if(slopes && curve)
  288.         {fprintf(stderr,"ERROR: slopes can\'t be specified for general curve\n");
  289.         exit(1);
  290.         }
  291.     sigma *= -1.;
  292.     if(xlog && !curve)
  293.         {if(x_arguments>1) xmax=mylog(xmax);
  294.         if(x_arguments>=1) xmin=mylog(xmin);
  295.         }
  296.     if(automatic_abscissas && abscissa_arguments<2 && x_arguments>=1)
  297.         abscissa=xmin;
  298.     p_data[0]=-1;
  299.     skipping=2;
  300.     if(automatic_abscissas) skipping--;
  301.     if(threed) skipping++;
  302.     i=0;
  303.     while(i<ENTRIES)
  304.         {if(fgets(buf,BUFSIZE,file)==0)break;
  305.         t=buf;
  306.         while(*t && isspace(*t)) t++;
  307.         if(*t == 0) continue;        /* skip blank lines */
  308.         buf[strlen(buf)-1]=0;        /* zero out the line feed */
  309.         if(buf[0]==';') {printf("%s\n",buf); continue;}  /* skip comment */
  310.         if(automatic_abscissas)
  311.             {x[i]=abscissa;
  312.             abscissa+=abscissa_step;
  313.             if(threed) sscanf(buf,"%F %F",&y[i],&z[i]);
  314.             else sscanf(buf,"%F",&y[i]);
  315.             if(ylog) y[i]=mylog(y[i]);
  316.             if(zlog) z[i]=mylog(z[i]);
  317.             }
  318.         else
  319.             {if(threed) sscanf(buf,"%F %F %F",&x[i],&y[i],&z[i]);
  320.             else sscanf(buf,"%F %F",&x[i],&y[i]);
  321.             if(xlog) x[i]=mylog(x[i]);
  322.             if(ylog) y[i]=mylog(y[i]);
  323.             if(zlog) z[i]=mylog(z[i]);
  324.             }
  325.         s=buf;                      /* start looking for label */
  326.         for (j=skipping; j--; )
  327.             {while(*s==' ')s++;            /* skip a number */
  328.             while(*s && (*s!=' '))s++;
  329.             }
  330.         while(*s==' ')s++;
  331.         if((length=strlen(s))&&(labels<MAXLABELS))
  332.             {if(*s=='\"')           /* label in quotes */
  333.                 {t=s+1;
  334.                 while(*t && (*t!='\"')) t++;
  335.                 t++;
  336.                 }
  337.             else                    /* label without quotes */
  338.                 {t=s;
  339.                 while(*t && (*t!=' '))t++;
  340.                 }
  341.             *t=0;
  342.             length=t-s;
  343.             p_data[labels]=i;
  344.             p_text[labels]=malloc(length+1);
  345.             if(p_text[labels]) strcpy(p_text[labels++],s);
  346.             }
  347.         i++;
  348.         }
  349.     n=i;
  350.     if(breaking && (!labels || p_data[labels-1]!=n-1))
  351.         {p_data[labels]=i-1;
  352.         if(p_text[labels]=malloc(1)) *p_text[labels++]=0;
  353.         }
  354.     yp=malloc(n*sizeof(double));
  355.     temp=malloc(n*sizeof(double));
  356.     if(temp==0 || yp==0)
  357.         {fprintf(stderr,"can\'t allocate buffer"); exit();
  358.         }
  359.     if(threed) 
  360.         {zp=malloc(n*sizeof(double));
  361.         if(zp==0) {fprintf(stderr,"can\'t allocate buffer"); exit();}
  362.         }
  363.     if(curve)
  364.         {xp=malloc(n*sizeof(double));
  365.         path=malloc(n*sizeof(double));
  366.         if(xp==0|| path==0) {fprintf(stderr,"can\'t allocate buffer"); exit();}
  367.         path[0]=sum=0.;
  368.         if(threed)
  369.             for (i=1; i<n; i++)
  370.                 {xx=x[i]-x[i-1];
  371.                 yy=y[i]-y[i-1];
  372.                 zz=z[i]-z[i-1];
  373.                 path[i]=(sum+=sqrt(xx*xx + yy*yy + zz*zz));
  374.                 }
  375.         else
  376.             for (i=1; i<n; i++)
  377.                 {xx=x[i]-x[i-1];
  378.                 yy=y[i]-y[i-1];
  379.                 path[i]=(sum+=sqrt(xx*xx + yy*yy));
  380.                 }
  381. /*        for(i=0; i<n; i++)
  382.             printf("path[%d]=%g  x[%d]=%g \n",i,path[i],i,x[i]); */
  383.         }
  384. }
  385.  
  386.  
  387. /* get_parameter - process one command line option
  388.         (return # parameters used) */
  389. get_parameter(argc,argv) int argc; char **argv;
  390. {    int i;
  391.     if(streq(*argv,"-a"))
  392.         {i=get_double(argc,argv,2,&abscissa_step,&abscissa,&abscissa);
  393.         abscissa_arguments=i-1;
  394.         automatic_abscissas=1;
  395.         return i;
  396.         }
  397.     else if(streq(*argv,"-b")) {breaking=1; return 1;}
  398.     else if(streq(*argv,"-c")) {curve=1; return 1;}
  399.     else if(streq(*argv,"-d")) {debugging=1; return 1;}
  400.     else if(streq(*argv,"-i"))
  401.         {if(argc>1)
  402.             {afile=fopen(argv[1],"r");
  403.             if(afile==0)
  404.                 {fwrite(stderr,"can\'t open abscissa file %s",argv[1]);
  405.                 exit(1);
  406.                 }
  407.             given_abscissas=1;
  408.             }
  409.         return 2;
  410.         }
  411.     else if(streq(*argv,"-n"))
  412.         {if((argc>1) && numeric(argv[1])) total=atoi(argv[1]);
  413.         return 2;
  414.         }
  415.     else if(streq(*argv,"-q")) {quadruple=1; return 1;}
  416.     else if(streq(*argv,"-t"))
  417.         {return(get_double(argc,argv,1,&sigma,&abscissa,&abscissa));
  418.         }
  419.     else if(streq(*argv,"-s"))
  420.         {slopes++;
  421.         return(get_double(argc,argv,2,&slope1,&slopen,&slopen));
  422.         }
  423.     else if(streq(*argv,"-x"))
  424.         {i=get_double(argc,argv,2,&xmin,&xmax,&xmax);
  425.         x_arguments=i-1;
  426.         return i;
  427.         }
  428.     else if(streq(*argv,"-xl")) {xlog++; return 1;}
  429.     else if(streq(*argv,"-yl")) {ylog++; return 1;}
  430.     else if(streq(*argv,"-zl")) {zlog++; threed++; return 1;}
  431.     else if(streq(*argv,"-3")) {threed++; return 1;}
  432.     else gripe(argv);
  433. }
  434.  
  435. get_double(argc,argv,permitted,a,b,c)
  436. int argc,permitted; char **argv; double *a,*b,*c;
  437. {    int i=1;
  438.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *a=atof(argv[i++]);
  439.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *b=atof(argv[i++]);
  440.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *c=atof(argv[i++]);
  441.     return i;
  442. }
  443.  
  444. int streq(a,b) char *a,*b;
  445. {    while(*a)
  446.         {if(*a!=*b)return 0;
  447.         a++; b++;
  448.         }
  449.     return 1;
  450. }
  451.  
  452. gripe_arg(s) char *s;
  453. {    fprintf(stderr,"argument missing for switch %s",s);
  454.     help();
  455. }
  456.  
  457. gripe(argv) char **argv;
  458. {    fprintf(stderr,*argv); fprintf(stderr," isn\'t a legal argument \n\n");
  459.     help();
  460. }
  461.  
  462. help()
  463. {    fprintf(stderr,"spline   version %s",VERSION);
  464.     fprintf(stderr,"\nusage: spline [file][options]\n");
  465.     fprintf(stderr,"options are:\n");
  466.     fprintf(stderr,"  -a  [step [start]] automatic abscissas \n");
  467.     fprintf(stderr,"  -b                 break interpolation after each label \n");
  468.     fprintf(stderr,"  -c                 general curve \n");
  469.     fprintf(stderr,"  -i  file           interpolate at abscissas from the file\n");
  470.     fprintf(stderr,"  -n  num            interpolate over num intervals \n");
  471.     fprintf(stderr,"  -q                 increase intervals fourfold \n");
  472.     fprintf(stderr,"  -s  [num [num]]    specify slopes at beginning and end \n");
  473.     fprintf(stderr,"  -t  num            tension in interpolating curve \n");
  474.     fprintf(stderr,"  -x  [min [max]]    interpolate from min to max \n");
  475.     fprintf(stderr,"  -xl                take logs of x values before interpolating \n");
  476.     fprintf(stderr,"  -yl                take logs of y values before interpolating \n");
  477.     fprintf(stderr,"  -zl                take logs of z values before interpolating \n");
  478.     fprintf(stderr,"                     (implies -3)\n");
  479.     fprintf(stderr,"  -3                 3D curve: x, y, and z given for each points\n");
  480.     exit();
  481. }
  482.  
  483. numeric(s) char *s;
  484. {    char c;
  485.     while(c=*s++)
  486.         {if((c<='9' && c>='0') || c=='+' || c=='-' || c=='.') continue;
  487.         return 0;
  488.         }
  489.     return 1;
  490. }
  491.  
  492. curv1(x,y,yp,n) double *x,*y,*yp; int n;
  493. {    int i;
  494.     double c1,c2,c3,deln,delnm1,delnn,dels,delx1,delx2,delx12;
  495.     double diag1,diag2,diagin,dx1,dx2,exps;
  496.     double sigmap,sinhs,sinhin,slpp1,slppn,spdiag;
  497.  
  498.     delx1=x[1] - x[0];
  499.     dx1=(y[1] - y[0])/delx1;
  500.     if(slopes) {slpp1=slp1; slppn=slpn;}
  501.     else
  502.         {if(n!=2)
  503.             {delx2= x[2] - x[1];
  504.             delx12= x[2] - x[0];
  505.             c1= -(delx12 + delx1)/delx12/delx1;
  506.             c2= delx12/delx1/delx2;
  507.             c3= -delx1/delx12/delx2;
  508.             slpp1= c1*y[0] + c2*y[1] + c3*y[2];
  509.             deln= x[n-1] - x[n-2];
  510.             delnm1= x[n-2] - x[n-3];
  511.             delnn= x[n-1] - x[n-3];
  512.             c1= (delnn + deln)/delnn/deln;
  513.             c2= -delnn/deln/delnm1;
  514.             c3= deln/delnn/delnm1;
  515.             slppn= c3*y[n-3] + c2*y[n-2] + c1*y[n-1];
  516.             }
  517.         else yp[0]=yp[1]=0.;
  518.         }
  519.                         /* denormalize tension factor */
  520.     sigmap=fabs(sigma)*(n-1)/(x[n-1]-x[0]);
  521.                         /* set up right hand side and tridiagonal system for
  522.                            yp and perform forward elimination                */
  523.     dels=sigmap*delx1;
  524.     exps=exp(dels);
  525.     sinhs=.5*(exps-1./exps);
  526.     sinhin=1./(delx1*sinhs);
  527.     diag1=sinhin*(dels*.5*(exps+1./exps)-sinhs);
  528.     diagin=1./diag1;
  529.     yp[0]=diagin*(dx1-slpp1);
  530.     spdiag=sinhin*(sinhs-dels);
  531.     temp[0]=diagin*spdiag;
  532.     if(n!=2)
  533.         {for(i=1; i<=n-2; i++)
  534.             {delx2=x[i+1]-x[i];
  535.             dx2=(y[i+1]-y[i])/delx2;
  536.             dels=sigmap*delx2;
  537.             exps=exp(dels);
  538.             sinhs=.5*(exps-1./exps);
  539.             sinhin=1./(delx2*sinhs);
  540.             diag2=sinhin*(dels*(.5*(exps+1./exps))-sinhs);
  541.             diagin=1./(diag1+diag2-spdiag*temp[i-1]);
  542.             yp[i]=diagin*(dx2-dx1-spdiag*yp[i-1]);
  543.             spdiag=sinhin*(sinhs-dels);
  544.             temp[i]=diagin*spdiag;
  545.             dx1=dx2;
  546.             diag1=diag2;
  547.             }
  548.         }
  549.     diagin=1./(diag1-spdiag*temp[n-2]);
  550.     yp[n-1]=diagin*(slppn-dx2-spdiag*yp[n-2]);
  551.                     /* perform back substitution */
  552.     for (i=n-2; i>=0; i--) yp[i] -= temp[i]*yp[i+1];
  553. }
  554.  
  555.  
  556. double curv2(x,y,yp,t,n) double *x,*y,*yp,t; int n;
  557. {    int i,j;
  558.     static int i1=1;
  559.     double del1,del2,dels,exps,exps1,s,sigmap,sinhd1,sinhd2,sinhs;
  560.  
  561.     s=x[n-1]-x[0];
  562.     sigmap=fabs(sigma)*(n-1)/s;
  563.     i=i1;
  564.     if(i>n) i=1;        /* want: x[i-1] <= t < x[i], 0 < i <= n */
  565.     while(i<n && t>=x[i]) i++;
  566.     while(i>1 && x[i-1]>t) i--;
  567.     i1=i;
  568.     del1=t-x[i-1];
  569.     del2=x[i]-t;
  570.     dels=x[i] - x[i-1];
  571.     exps1=exp(sigmap*del1); sinhd1=.5*(exps1-1./exps1);
  572.     exps= exp(sigmap*del2); sinhd2=.5*(exps-1./exps);
  573.     exps=exps1*exps;        sinhs=.5*(exps-1./exps);
  574.     return ((yp[i]*sinhd1 + yp[i-1]*sinhd2)/sinhs +
  575.             ((y[i] - yp[i])*del1 + (y[i-1] - yp[i-1])*del2)/dels);
  576. }
  577.  
  578. double mylog(x) double x;
  579. {    if(x>0.) return log(x);
  580.     fprintf(stderr,"can%'t take log of nonpositive number");
  581.     exit(1);
  582.     return 0.;
  583. }
  584.  
  585. #ifdef DESMET
  586. /*    write out the command line */
  587. fprint_cmd(fp, format) 
  588. FILE *fp;        /* pointer to output file */
  589. char *format;     /* desired print format, possibly including program name */
  590. {    char buf[129];
  591.     int i=0;
  592.     _lmove(129, 130, _showcs()-0x10, &buf, _showds());
  593.     buf[128]=13;
  594.     while(buf[i]!=13) i++;
  595.     buf[i]=0;
  596.     fprintf(fp, format, buf);
  597. }
  598.  
  599. #endif
  600.