home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Programmierung / SOURCE.mdf / programm / msdos / c / spline29 / spline.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-11  |  23.7 KB  |  827 lines
  1. /*    spline - interpolate points in a tabulated function
  2.  
  3.     Usage...
  4.         usage: spline  [infile]  [options]
  5.         or:    spline   infile  [outfile]  [options]
  6.         options are:
  7.           -a  [step [start]] automatic abscissas 
  8.           -b                 break interpolation after each label 
  9.           -c                 general curve 
  10.           -i  file           interpolate at abscissas from the file.  Only
  11.                                1st number on each line of file is used -
  12.                                rest of line is ignored.  Blank lines and
  13.                                lines starting with ';' are ignored.
  14.           -n  num            interpolate over num intervals 
  15.           -q                 increase intervals fourfold 
  16.           -s  [num [num]]    specify slopes at beginning and end 
  17.                                (where "n" implies "natural" boundary)
  18.           -t  num            tension in interpolating curve 
  19.           -x  [min [max]]    interpolate from min to max 
  20.           -xt -yt -zt        take atanh of x, y, or z values before
  21.                              interpolating. Data values are restricted
  22.                              to (0, 1).  (-zt implies -3)
  23.           -xl -yl -zl        take log of x, y, or z values before
  24.                              interpolating.  Data values are positive. 
  25.                              (-zt implies -3)
  26.           -3                 3D curve: x, y, and z given for each points
  27.  
  28.     Notes...
  29.         This program fits a smooth curve through a given set of points,
  30.         using the splines under tension introduced by Schweikert [J.
  31.         Math. and Physics v 45 pp 312-317 (1966)].  They are similar
  32.         to cubic splines, but are less likely to introduce spurious
  33.         inflection points.
  34.  
  35.     History...
  36.         11 Nov 92    Ver 2.9: -q yields EXACTLY four times the number
  37.                     of intervals.  ';' anywhere in a data line starts
  38.                     a comment.  Any line not starting with a number is
  39.                     a comment.  Abscissas can be either decreasing or
  40.                     increasing.
  41.         10 Nov 92    Ver 2.8: Implemented implicit interpolation
  42.                     using zero().
  43.         31 Oct 92    Ver 2.7: Can specify only one slope.        
  44.         18 Aug 92    Ver 2.6: Correctly handles two point curves.
  45.         11 Aug 92    Ver 2.5: Echos command line under Turbo C.  atanh 
  46.                     transformations implemented.
  47.         14 Jan 88    Ver 2.4: Fixed integer multiplication overflow in curv0 
  48.                     which ruined end of large output files.
  49.         8 Jan 88    Ver 2.3: Fixed off-by-one bug in curv2 which sometimes 
  50.                     permitted an out-of-range index near end of curve.
  51.                     Output file name can be specified on command line.
  52.                     Switches can appear anywhere on command line.
  53.         6 Jan 88    600 input points permitted.
  54.         -- ver 2.2 --
  55.         5 Jan 88    numbers printed out with fixed width of 15 characters
  56.         15 Jun 87    echoes command line to output file.
  57.         -- ver 2.1 --
  58.         25 Nov 86    3D
  59.         -- ver 2.0 --
  60.         16 Oct 86    accepts desired abscissas from file.
  61.         -- ver 1.4 --
  62.         7 Jun 86    echoes a 2-point line.
  63.         3 Jun 86    -s switch allows specified slopes at ends.
  64.         2 Jun 86    Ver 1.3: Fixed several embarrassing bugs in -b option.
  65.                     -q switch quadruples number of intervals.
  66.         6 Apr 86    -b switch breaks interpolation at labels.
  67.         21 Sep 85    Added -xl and -yl switches for taking logs
  68.         23 Nov 85    Passing lines starting with ';' unchanged, otherwise
  69.                     ignoring them.
  70.     Bugs...
  71.         Should have -e command for evenly spaced abscissas.
  72.         -q switch should accept argument for desired fractional
  73.         increase in # points.
  74.         Should accept the output file as an optional second argument.
  75.  
  76.     Author...
  77.         Copyright (c) 1985, 1986  James R. Van Zandt (jrv@mitre-bedford)
  78.  
  79.         All rights reserved.
  80.         This program may be copied for personal, non-profit use only.
  81.  
  82.         Based on algorithms by A.  K.  Cline [Communications of the ACM
  83.         v 17 n 4 pp 218-223 (Apr 74)].
  84.  
  85. */
  86.  
  87. #include <stdio.h>
  88. #include <string.h>
  89. #include <math.h>
  90.  
  91. #define VERSION "2.9"
  92.  
  93. #ifndef __TURBOC__
  94. #define __DESMET__
  95. #endif
  96.  
  97. char *RIGHT=        "Copyright (c) 1985 - 1992  James R. Van Zandt";
  98.  
  99. #define TR(a) /*  if(debugging) printf a   /* debugging printouts */
  100.  
  101. #define ENTRIES 600
  102. #define MAXLABELS 50
  103. #define BUFSIZE 200
  104.  
  105. double *x, *y, *z, *temp, *xp, *yp, *zp, *path;
  106. double slp1=0.,slpn=0.,sigma=1.;
  107. double abscissa=0.;
  108. double abscissa_step=1.;
  109. double slope1=0.,slopen=0.;    /* slopes supplied by user */
  110. double xmin=0.;
  111. double xmax=0.;
  112. double curv2(), mylog(), tanh(), atanh();
  113.  
  114. int abscissa_arguments=0;
  115. int automatic_abscissas=0;
  116. int breaking=0;        /* nonzero if breaking at labels */
  117. int curve=0;        /* nonzero if general curve permitted */
  118. int debugging=0;
  119. int given_x=0;        /* nonzero if x values are given in a file */
  120. int given_y=0;        /* nonzero if y values are given in a file */
  121. int given_z=0;        /* nonzero if z values are given in a file */
  122. int given_abscissas=0;    /* nonzero if abscissas are given in a file 
  123.         (given_abscissas and given_x are equivalent if curve=0) */
  124. int labels=0;        /* number of labels in input data */
  125. int quadruple=0;    /* nonzero if user asking for 4 times the resolution */
  126. int specified_slope1=0;    /* nonzero if first slope specified */
  127. int specified_slopen=0;    /* nonzero if final slope specified */
  128. int threed=0;        /* nonzero if 3D curve (x, y, and z given for each point) */
  129. int xlog=0;        /* nonzero if taking log of x values */
  130. int ylog=0;        /* nonzero if taking log of y values */
  131. int zlog=0;        /* nonzero if taking log of z values */
  132. int xtanh=0;        /* nonzero if taking tanh of x values */
  133. int ytanh=0;        /* nonzero if taking tanh of y values */
  134. int ztanh=0;        /* nonzero if taking tanh of z values */
  135.  
  136. int x_arguments=0;
  137. int n;                /* number of entries in x, y, yp, and temp */
  138. int index_array[MAXLABELS];    /* indexes into x and y */
  139. int *p_data=index_array;
  140. int total=100;
  141.  
  142. FILE *ifile=stdin, *afile, *ofile=stdout;
  143.  
  144. char buf[BUFSIZE];
  145. char *label_array[MAXLABELS];
  146. char **p_text=label_array;
  147. char newline[]="\n", nothing[]="", *leadin;
  148.  
  149. main(argc,argv) int argc; char **argv;
  150. {    int nn,origin;
  151.     read_data(argc,argv);
  152.     if(breaking && labels)
  153.         {origin=0;
  154.         while(labels--)
  155.             {p_data[0] -= origin;
  156.             nn=p_data[0]+1;
  157.             if(quadruple) total=(nn-1)*4;
  158.             if(nn) curv0(nn,total);
  159.             origin += nn;
  160.             path += nn;
  161.             x += nn;
  162.             y += nn;
  163.             p_data++;
  164.             p_text++;
  165.             }
  166.         }
  167.     else
  168.         {if(quadruple) total=(n-1)*4;
  169.         curv0(n,total);
  170.         }
  171.     if(ofile!=stdout) fclose(ofile);
  172. }
  173.  
  174. curv0(n,requested)
  175. int n,requested;
  176. {    int i, j, each, stop, seg=0, regressing=0;
  177.     double s, ds, xx, yy, zz;
  178.     char *t;
  179.  
  180.     if(n<3)
  181.         {for (j=0; j<n; j++)
  182.             {xx=x[j]; if(xlog) xx=exp(xx); if(xtanh) xx = tanh(xx);
  183.             yy=y[j]; if(ylog) yy=exp(yy); if(ytanh) yy = tanh(yy);
  184.             fprintf(ofile, "%15.6g %15.6g ",xx,yy);
  185.             if(threed)
  186.                 {zz=z[n-1]; if(zlog) zz=exp(zz); if(ztanh) zz = tanh(zz);
  187.                 fprintf(ofile, "%15.6g ",zz);
  188.                 }
  189.             if(j==p_data[seg]) fputs(p_text[seg++], ofile);
  190.             fputc('\n', ofile);
  191.             }
  192.         return;
  193.         }
  194.     if(path[n-1] < path[0])
  195.         {reverse_array(x, n);
  196.         reverse_array(y, n);
  197.         if(threed) reverse_array(z, n);
  198.         }
  199.     for(i=1; i<n; i++) if(path[i]<=path[i-1]) regressing++;
  200.     if (regressing)
  201.         {fprintf(stderr,"ERROR: independent variable not strictly increasing\n");
  202.         return;
  203.         }
  204.     if (curve && (xmin<0. || xmax>1.))
  205.         {fprintf(stderr,"ERROR: independent variable not in range 0 to 1\n");
  206.         return;
  207.         }
  208.     if(curve)
  209.         {curv1(path,x,xp,n); curv1(path,y,yp,n);
  210.         if(threed) curv1(path,z,zp,n);
  211.         }
  212.     else
  213.         {slp1=slope1; slpn=slopen; curv1(x,y,yp,n);
  214.         if(threed) 
  215.             {slp1=slpn=0.; curv1(x,z,zp,n);
  216.             }
  217.         }
  218.     s=path[0];
  219.     seg=0;
  220.     if(requested<n) requested=n;
  221.     if(given_abscissas || given_x || given_y || given_z)
  222.         {                            /* abscissas specified in a file */
  223.         int mini, maxi, numbers;
  224.         double *target, *given, given_min, given_max, dif, tolerance, p=0.;
  225.         char *inversion=0, *init_zero();
  226.  
  227.         leadin=nothing;
  228.         if(given_x){target = &xx; given = x;}
  229.         if(given_y){target = &yy; given = y;}
  230.         if(given_z){target = &zz; given = z;}
  231.         if(given_x || given_y || given_z)
  232.             {given_min = given_max = given[0];
  233.             mini = maxi = 0;
  234.             for (i = 1; i < n; i++)
  235.                 {if(given[i] < given_min) given_min = given[mini=i];
  236.                 if(given[i] > given_max) given_max = given[maxi=i];
  237.                 }
  238.             }
  239.         while(fgets(buf,BUFSIZE,afile))
  240.             {if(t = strchr(buf, ';')) 
  241.                 *t = 0;    /* ';' starts a comment */
  242.             t=buf;
  243.             while(*t && isspace(*t)) t++;
  244.             if(*t == 0) continue;        /* skip blank lines */
  245.             numbers=sscanf(buf,"%lf",&s);
  246.             if(numbers==0) continue;    /* skip comment line */
  247.  
  248.             if(given_abscissas)
  249.                 {if(xlog) s=mylog(s);
  250.                 if(xtanh) s = atanh(s);
  251.                 if(s<path[0] || path[n-1]<s) continue; /* don't extrapolate */
  252.                 value(s, &xx, &yy, &zz, n);
  253.                 goto PRINTIT;
  254.                 }
  255.             if(given_x)
  256.                 {if(xlog) s = mylog(s);
  257.                 else if(xtanh) s = atanh(s);
  258.                 }
  259.             if(given_y)
  260.                 {if(ylog) s = mylog(s);
  261.                 else if(ytanh) s = atanh(s);
  262.                 }
  263.             if(given_z)
  264.                 {if(zlog) s = mylog(s);
  265.                 else if(ztanh) s = atanh(s);
  266.                 }
  267.             if((given_max - s)*(given_min - s) > 0.) continue;    /* hopeless */
  268.             inversion = init_zero(&p, inversion);
  269.             p = path[mini]; zero(given_min - s, inversion);
  270.             p = path[maxi]; zero(given_max - s, inversion);
  271.             if(fabs(s) < 1.e-10) tolerance = 1.e-20;
  272.             else tolerance = fabs(s)*1.e-6;
  273.             do
  274.                 {value(p, &xx, &yy, &zz, n);
  275.                 dif = *target - s;
  276.                 if(fabs(dif) <= tolerance) goto PRINTIT;
  277.                 } while(!zero(dif, inversion));
  278.             continue;    /* search failed */
  279. PRINTIT:
  280.  
  281.             if(threed) 
  282.                 fprintf(ofile, "%s%15.6g %15.6g %15.6g", leadin, xx, yy, zz);
  283.             else 
  284.                 fprintf(ofile, "%s%15.6g %15.6g", leadin, xx, yy);
  285.             leadin=newline;
  286.             }
  287.         if(breaking && labels) fprintf(ofile, " %s\n",p_text[0]);
  288.         else fprintf(ofile, " \"\"\n");
  289.         rewind(afile);
  290.         }
  291.     else if(x_arguments==2)            /* specific range was requested */
  292.         {ds=(xmax-xmin)/requested;
  293.         if(curve)
  294.             {ds*=(path[n-1]-path[0]);
  295.             s=xmin*(path[n-1]-path[0])+path[0];
  296.             }
  297.         else s=xmin;
  298.         for(i=requested+1; i; i--)
  299.             {value(s,&xx,&yy,&zz,n);
  300.             if(threed) fprintf(ofile, "%15.6g %15.6g %15.6g ", xx, yy, zz);
  301.             else fprintf(ofile, "%15.6g %15.6g ", xx, yy);
  302.             if(j==p_data[seg]) fputs(p_text[seg++], ofile);
  303.             fputc('\n', ofile);
  304.             s+=ds;
  305.             }
  306.         }
  307.     else                        /* spline for entire curve was requested */
  308.         {for (i=j=0; j<n-1; j++)
  309.             {stop=(int)((long)requested*(j+1)/(n-1));
  310.             ds=(path[j+1]-path[j])/(stop-i);
  311.             for ( ; i<stop; i++)
  312.                 {value(s,&xx,&yy,&zz,n);
  313.                 if(threed) fprintf(ofile, "%15.6g %15.6g %15.6g ", xx, yy, zz);
  314.                 else fprintf(ofile, "%15.6g %15.6g ",xx,yy);
  315.                 if(j==p_data[seg]) fputs(p_text[seg++], ofile);
  316.                 fputc('\n', ofile);
  317.                 s+=ds;
  318.                 }
  319.             }
  320.         xx=x[n-1]; if(xlog) xx=exp(xx); if(xtanh) xx = tanh(xx);
  321.         yy=y[n-1]; if(ylog) yy=exp(yy); if(ytanh) yy = tanh(yy);
  322.         if(threed)
  323.             {zz=z[n-1]; if(zlog) zz=exp(zz); if(ztanh) zz = tanh(zz);
  324.             fprintf(ofile, "%15.6g %15.6g %15.6g ", xx, yy, zz);
  325.             }
  326.         else fprintf(ofile, "%15.6g %15.6g ",xx,yy);
  327.         if(j==p_data[seg]) fputs(p_text[seg++], ofile);
  328.         fputc('\n', ofile);
  329.         }
  330. }
  331.  
  332. value(s,q,r,v,n) double s,*q,*r,*v; int n;
  333. {    TR(("\nvalue(%f,,,%d)\n", s, n));
  334.     if(curve)
  335.         {*q=curv2(path,x,xp,s,n);
  336.         *r=curv2(path,y,yp,s,n);
  337.         if(threed) *v=curv2(path,z,zp,s,n);
  338.         }
  339.     else
  340.         {*q=s;
  341.         *r=curv2(x,y,yp,s,n);
  342.         if(threed) *v=curv2(x,z,zp,s,n);
  343.         }
  344.     if(xlog) *q=exp(*q); if(xtanh) *q = tanh(*q);
  345.     if(ylog) *r=exp(*r); if(ytanh) *r = tanh(*r);
  346.     if(zlog) *v=exp(*v); if(ztanh) *v = tanh(*v);
  347. }
  348.  
  349. char *gmem(n) int n;
  350. {    char *p, *malloc();
  351.     p=malloc(n);
  352.     if(p==0)
  353.         {fprintf(stderr, "out of memory"); exit(1);
  354.         }
  355.     return p;
  356. }
  357.  
  358. reverse_array(x, n) double x[]; int n;
  359. {
  360.     int i, j;
  361.     double t;
  362.     j = n - 1;
  363.     for (i = 0; i < j; i++, j--)
  364.         {
  365.         t = x[i]; x[i] = x[j]; x[j] = t;
  366.         }
  367. }
  368.  
  369. read_data(argc,argv) int argc; char **argv;
  370. {    int i, j, length, skipping, ac, numbers;
  371.     double xx, yy, zz, ss, d, *pd, sum;
  372.     char *s, *t, **av;
  373.  
  374.     x=path=(double *)gmem(ENTRIES*sizeof(double));
  375.     y=(double *)gmem(ENTRIES*sizeof(double));
  376.     argc--; argv++;
  377. /*        handle switches on command line */
  378.     ac=argc; av=argv; argc=0;
  379.     while(ac>0)
  380.         {if(**av=='?') help();
  381.         if(**av=='-')
  382.             {i=get_parameter(ac, av);
  383.             av += i; ac -= i;
  384.             }
  385.         else {argv[argc++] = *av++; ac--;}
  386.         }
  387.     if(argc)
  388.         {ifile=fopen(*argv,"r");
  389.         if(ifile==0) {fprintf(stderr, "file %s not found\n",*argv); exit();}
  390.         argc--; argv++;
  391.         }
  392.     if(argc)
  393.         {unlink(*argv);
  394.         ofile=fopen(*argv,"w");
  395.         if(ofile==0) {fprintf(stderr, "can't open file %s\n",*argv); exit();}
  396.         argc--; argv++;
  397.         }
  398.  
  399.     fprint_cmd(ofile, "; spline %s\n");
  400.  
  401.     if(given_x && !curve){given_abscissas = 1; given_x = 0;}
  402.     if(given_abscissas && curve){given_x = 1; given_abscissas = 0;}
  403.     if(threed)
  404.         z=(double *)gmem(ENTRIES*sizeof(double));
  405.     if(sigma<=0.)
  406.         {fprintf(stderr,"ERROR: tension must be positive\n");
  407.         exit(1);
  408.         }
  409.     if((specified_slope1 || specified_slopen) && curve)
  410.         {fprintf(stderr,"ERROR: slopes can\'t be specified for general curve\n");
  411.         exit(1);
  412.         }
  413.     sigma *= -1.;
  414.     if(xlog && !curve)
  415.         {if(x_arguments>1) xmax=mylog(xmax);
  416.         if(x_arguments>=1) xmin=mylog(xmin);
  417.         }
  418.     if(xtanh && !curve)
  419.         {if(x_arguments>1) xmax=atanh(xmax);
  420.         if(x_arguments>=1) xmin=atanh(xmin);
  421.         }
  422.     if(automatic_abscissas && abscissa_arguments<2 && x_arguments>=1)
  423.         abscissa=xmin;
  424.     p_data[0]=-1;
  425.     skipping=2;
  426.     if(automatic_abscissas) skipping--;
  427.     if(threed) skipping++;
  428.     i=0;
  429.     while(i<ENTRIES)
  430.         {if(fgets(buf,BUFSIZE,ifile)==0)break;
  431.         if(buf[0]==';') {fprintf(ofile, "%s",buf); continue;}  /* echo comment */
  432.         buf[strlen(buf)-1]=0;                /* zero out the line feed */
  433.         if(t = strchr(buf, ';')) *t = 0;    /* ';' starts a comment */
  434.         t=buf;
  435.         while(*t && isspace(*t)) t++;
  436.         if(*t == 0) continue;                /* skip blank lines */
  437.         if(automatic_abscissas)
  438.             {x[i]=abscissa;
  439.             abscissa+=abscissa_step;
  440.             if(threed) numbers=sscanf(buf,"%lf %lf",&y[i],&z[i]);
  441.             else numbers=sscanf(buf,"%lf",&y[i]);
  442.             if(numbers==0) continue;
  443.             if(ylog) y[i]=mylog(y[i]);
  444.             if(ytanh) y[i]=atanh(y[i]);
  445.             if(zlog) z[i]=mylog(z[i]);
  446.             if(ztanh) z[i]=atanh(z[i]);
  447.             }
  448.         else
  449.             {if(threed) numbers=sscanf(buf,"%lf %lf %lf",&x[i],&y[i],&z[i]);
  450.             else numbers=sscanf(buf,"%lf %lf",&x[i],&y[i]);
  451.             if(numbers==0) continue;
  452.             if(xlog) x[i]=mylog(x[i]); if(xtanh) x[i] = atanh(x[i]);
  453.             if(ylog) y[i]=mylog(y[i]); if(ytanh) y[i] = atanh(y[i]);
  454.             if(zlog) z[i]=mylog(z[i]); if(ztanh) z[i] = atanh(z[i]);
  455.             }
  456.         s=buf;                      /* start looking for label */
  457.         for (j=skipping; j--; )
  458.             {while(*s==' ')s++;            /* skip a number */
  459.             while(*s && (*s!=' '))s++;
  460.             }
  461.         while(*s==' ')s++;
  462.         if((length=strlen(s))&&(labels<MAXLABELS))
  463.             {if(*s=='\"')           /* label in quotes */
  464.                 {t=s+1;
  465.                 while(*t && (*t!='\"')) t++;
  466.                 t++;
  467.                 }
  468.             else                    /* label without quotes */
  469.                 {t=s;
  470.                 while(*t && (*t!=' '))t++;
  471.                 }
  472.             *t=0;
  473.             length=t-s;
  474.             p_data[labels]=i;
  475.             p_text[labels]=gmem(length+1);
  476.             strcpy(p_text[labels++],s);
  477.             }
  478.         i++;
  479.         }
  480.     n=i;
  481.     if(n==ENTRIES) 
  482.         fprintf(stderr, "WARNING: only the first %d input points are used\n", 
  483.             ENTRIES);
  484.     if(breaking && (!labels || p_data[labels-1]!=n-1))
  485.         {p_data[labels]=i-1;
  486.         p_text[labels]=gmem(1);
  487.         *p_text[labels++]=0;
  488.         }
  489.     yp=(double *)gmem(n*sizeof(double));
  490.     temp=(double *)gmem(n*sizeof(double));
  491.     if(threed) 
  492.         zp=(double *)gmem(n*sizeof(double));
  493.     if(curve)
  494.         {xp=(double *)gmem(n*sizeof(double));
  495.         path=(double *)gmem(n*sizeof(double));
  496.         path[0]=sum=0.;
  497.         if(threed)
  498.             for (i=1; i<n; i++)
  499.                 {xx=x[i]-x[i-1];
  500.                 yy=y[i]-y[i-1];
  501.                 zz=z[i]-z[i-1];
  502.                 path[i]=(sum+=sqrt(xx*xx + yy*yy + zz*zz));
  503.                 }
  504.         else
  505.             for (i=1; i<n; i++)
  506.                 {xx=x[i]-x[i-1];
  507.                 yy=y[i]-y[i-1];
  508.                 path[i]=(sum+=sqrt(xx*xx + yy*yy));
  509.                 }
  510. /*        for(i=0; i<n; i++)
  511.             printf("path[%d]=%g  x[%d]=%g \n",i,path[i],i,x[i]); */
  512.         }
  513. }
  514.  
  515.  
  516. /* get_parameter - process one command line option
  517.         (return # parameters used) */
  518. get_parameter(argc,argv) int argc; char **argv;
  519. {    int i;
  520.     if(streq(*argv,"-a"))
  521.         {i=get_double(argc,argv,2,&abscissa_step,&abscissa,&abscissa);
  522.         abscissa_arguments=i-1;
  523.         automatic_abscissas=1;
  524.         return i;
  525.         }
  526.     else if(streq(*argv,"-b")) {breaking=1; return 1;}
  527.     else if(streq(*argv,"-c")) {curve=1; return 1;}
  528.     else if(streq(*argv,"-d")) {debugging=1; return 1;}
  529.     else if(streq(*argv,"-i"))
  530.         {if(argc>1)
  531.             {afile=fopen(argv[1],"r");
  532.             if(afile==0)
  533.                 {fprintf(stderr,"can\'t open abscissa file %s",argv[1]);
  534.                 exit(1);
  535.                 }
  536.             given_abscissas=1;
  537.             }
  538.         return 2;
  539.         }
  540.     else if(streq(*argv,"-ix"))
  541.         {if(argc>1)
  542.             {afile=fopen(argv[1],"r");
  543.             if(afile==0)
  544.                 {fprintf(stderr,"can\'t open abscissa file %s",argv[1]);
  545.                 exit(1);
  546.                 }
  547.             given_x=1;
  548.             }
  549.         return 2;
  550.         }
  551.     else if(streq(*argv,"-iy"))
  552.         {if(argc>1)
  553.             {afile=fopen(argv[1],"r");
  554.             if(afile==0)
  555.                 {fprintf(stderr,"can\'t open abscissa file %s",argv[1]);
  556.                 exit(1);
  557.                 }
  558.             given_y=1;
  559.             }
  560.         return 2;
  561.         }
  562.     else if(streq(*argv,"-iz"))
  563.         {if(argc>1)
  564.             {afile=fopen(argv[1],"r");
  565.             if(afile==0)
  566.                 {fprintf(stderr,"can\'t open abscissa file %s",argv[1]);
  567.                 exit(1);
  568.                 }
  569.             threed=given_z=1;
  570.             }
  571.         return 2;
  572.         }
  573.     else if(streq(*argv,"-n"))
  574.         {if((argc>1) && numeric(argv[1])) total=atoi(argv[1]);
  575.         return 2;
  576.         }
  577.     else if(streq(*argv,"-q")) {quadruple=1; return 1;}
  578.     else if(streq(*argv,"-t"))
  579.         {return(get_double(argc,argv,1,&sigma,&abscissa,&abscissa));
  580.         }
  581.     else if(streq(*argv,"-s"))
  582.         {if(streq(argv[1], "n")) i = 2;
  583.         else            
  584.             {i = get_double(argc, argv, 1, &slope1,&slopen,&slopen);
  585.             if(i == 1) return 1;
  586.             specified_slope1 = 1;
  587.             }
  588.         argv++;
  589.         argc--;
  590.         if(streq(argv[1], "n")) i++;
  591.         else            
  592.             {i += get_double(argc, argv, 1, &slopen,&slopen,&slopen) - 1;
  593.             if(i > 2) specified_slopen = 1;
  594.             }
  595.         return i;
  596.         }
  597.     else if(streq(*argv,"-x"))
  598.         {i=get_double(argc,argv,2,&xmin,&xmax,&xmax);
  599.         x_arguments=i-1;
  600.         return i;
  601.         }
  602.     else if(streq(*argv,"-xl")) {xlog=1; return 1;}
  603.     else if(streq(*argv,"-yl")) {ylog=1; return 1;}
  604.     else if(streq(*argv,"-zl")) {zlog=1; threed=1; return 1;}
  605.     else if(streq(*argv,"-xt")) {xtanh=1; return 1;}
  606.     else if(streq(*argv,"-yt")) {ytanh=1; return 1;}
  607.     else if(streq(*argv,"-zt")) {ztanh=1; threed=1; return 1;}
  608.     else if(streq(*argv,"-3")) {threed=1; return 1;}
  609.     else gripe(argv);
  610. }
  611.  
  612. get_double(argc,argv,permitted,a,b,c)
  613. int argc,permitted; char **argv; double *a,*b,*c;
  614. {    int i=1;
  615.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *a=atof(argv[i++]);
  616.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *b=atof(argv[i++]);
  617.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *c=atof(argv[i++]);
  618.     return i;
  619. }
  620.  
  621. int streq(a,b) char *a,*b;
  622. {    while(*a)
  623.         {if(*a!=*b)return 0;
  624.         a++; b++;
  625.         }
  626.     return 1;
  627. }
  628.  
  629. gripe_arg(s) char *s;
  630. {    fprintf(stderr,"argument missing for switch %s",s);
  631.     help();
  632. }
  633.  
  634. gripe(argv) char **argv;
  635. {    fprintf(stderr,*argv); fprintf(stderr," isn\'t a legal argument \n\n");
  636.     help();
  637. }
  638.  
  639. char *message[]=
  640. {"spline   version ", VERSION,
  641. "\nusage: spline  [infile]  [options]\n",
  642. "or:    spline   infile  [outfile]  [options]\n",
  643. "options are:\n",
  644. "  -a  [step [start]] automatic abscissas \n",
  645. "  -b                 break interpolation after each label \n",
  646. "  -c                 general curve \n",
  647. "  -i  file           explicit interpolation at abscissas from the file\n",
  648. "  -ix file\n",
  649. "  -iy file\n",
  650. "  -iz file           implicit interpolation at values from the file\n",
  651. "  -n  num            interpolate over num intervals \n",
  652. "  -q                 increase intervals fourfold \n",
  653. "  -s  [num [num]]    specify slopes at beginning and end \n",
  654. "                     (\"n\" for \"natural\")\n",
  655. "  -t  num            tension in interpolating curve \n",
  656. "  -x  [min [max]]    interpolate from min to max \n",
  657. "  -xt -yt -zt        take atanh of x, y, or z values before\n",
  658. "                     interpolating. Data values are restricted\n",
  659. "                     to (0, 1).  (-zt implies -3)\n",
  660. "  -xl -yl -zl        take log of x, y, or z values before\n",
  661. "                     interpolating.  Data values are positive. \n",
  662. "                     (-zt implies -3)\n",
  663. "  -3                 3D curve: x, y, and z given for each points\n",
  664. 0
  665. };
  666.  
  667. help()
  668. {    char **s;
  669.     for (s=message; *s; s++) printf(*s);
  670.     exit();
  671. }
  672.  
  673. numeric(s) char *s;
  674. {    char c;
  675.     while(c=*s++)
  676.         {if((c<='9' && c>='0') || c=='+' || c=='-' || c=='.') continue;
  677.         return 0;
  678.         }
  679.     return 1;
  680. }
  681.  
  682. curv1(x,y,yp,n) double *x,*y,*yp; int n;
  683. {    int i;
  684.     double c1,c2,c3,deln,delnm1,delnn,dels,delx1,delx2,delx12;
  685.     double diag1,diag2,diagin,dx1,dx2,exps;
  686.     double sigmap,sinhs,sinhin,slpp1,slppn,spdiag;
  687.  
  688.     delx1=x[1] - x[0];
  689.     dx1=(y[1] - y[0])/delx1;
  690.  
  691.     if(specified_slope1) slpp1 = slp1;
  692.     else if(n!=2)
  693.         {delx2= x[2] - x[1];
  694.         delx12= x[2] - x[0];
  695.         c1= -(delx12 + delx1)/delx12/delx1;
  696.         c2= delx12/delx1/delx2;
  697.         c3= -delx1/delx12/delx2;
  698.         slpp1= c1*y[0] + c2*y[1] + c3*y[2];
  699.         }
  700.     else yp[0] = yp[1] = 0.;
  701.  
  702.     if(specified_slopen) slppn = slpn;
  703.     else if(n!=2)
  704.         {deln= x[n-1] - x[n-2];
  705.         delnm1= x[n-2] - x[n-3];
  706.         delnn= x[n-1] - x[n-3];
  707.         c1= (delnn + deln)/delnn/deln;
  708.         c2= -delnn/deln/delnm1;
  709.         c3= deln/delnn/delnm1;
  710.         slppn= c3*y[n-3] + c2*y[n-2] + c1*y[n-1];
  711.         }
  712.     else yp[0] = yp[1] = 0.;
  713.                         /* denormalize tension factor */
  714.     sigmap=fabs(sigma)*(n-1)/(x[n-1]-x[0]);
  715.                         /* set up right hand side and tridiagonal system for
  716.                            yp and perform forward elimination                */
  717.     dels=sigmap*delx1;
  718.     exps=exp(dels);
  719.     sinhs=.5*(exps-1./exps);
  720.     sinhin=1./(delx1*sinhs);
  721.     diag1=sinhin*(dels*.5*(exps+1./exps)-sinhs);
  722.     diagin=1./diag1;
  723.     yp[0]=diagin*(dx1-slpp1);
  724.     spdiag=sinhin*(sinhs-dels);
  725.     temp[0]=diagin*spdiag;
  726.     if(n!=2)
  727.         {for(i=1; i<=n-2; i++)
  728.             {delx2=x[i+1]-x[i];
  729.             dx2=(y[i+1]-y[i])/delx2;
  730.             dels=sigmap*delx2;
  731.             exps=exp(dels);
  732.             sinhs=.5*(exps-1./exps);
  733.             sinhin=1./(delx2*sinhs);
  734.             diag2=sinhin*(dels*(.5*(exps+1./exps))-sinhs);
  735.             diagin=1./(diag1+diag2-spdiag*temp[i-1]);
  736.             yp[i]=diagin*(dx2-dx1-spdiag*yp[i-1]);
  737.             spdiag=sinhin*(sinhs-dels);
  738.             temp[i]=diagin*spdiag;
  739.             dx1=dx2;
  740.             diag1=diag2;
  741.             }
  742.         }
  743.     diagin=1./(diag1-spdiag*temp[n-2]);
  744.     yp[n-1]=diagin*(slppn-dx2-spdiag*yp[n-2]);
  745.                     /* perform back substitution */
  746.     for (i=n-2; i>=0; i--) yp[i] -= temp[i]*yp[i+1];
  747. }
  748.  
  749.  
  750. double curv2(x,y,yp,t,n) double *x,*y,*yp,t; int n;
  751. {    int i,j;
  752.     static int i1=1;
  753.     double del1,del2,dels,exps,exps1,s,sigmap,sinhd1,sinhd2,sinhs;
  754.  
  755.     TR(("curv2( , , ,%f,%d)\n", t, n));
  756.     s=x[n-1]-x[0];
  757.     sigmap=fabs(sigma)*(n-1)/s;
  758.     i=i1;
  759.     if(i>n) i=1;        /* want: x[i-1] <= t < x[i], 0 < i < n */
  760.     while(i<n-1 && t>=x[i]) i++;
  761.     while(i>1 && x[i-1]>t) i--;
  762.     TR(("  i=%d: x[i-1]=%f <= t=%f < x[i]=%f\n", i, x[i-1], t, x[i]));
  763.     i1=i;
  764.     del1=t-x[i-1];
  765.     del2=x[i]-t;
  766.     dels=x[i] - x[i-1];
  767.     exps1=exp(sigmap*del1); sinhd1=.5*(exps1-1./exps1);
  768.     exps= exp(sigmap*del2); sinhd2=.5*(exps-1./exps);
  769.     exps=exps1*exps;        sinhs=.5*(exps-1./exps);
  770.     return ((yp[i]*sinhd1 + yp[i-1]*sinhd2)/sinhs +
  771.             ((y[i] - yp[i])*del1 + (y[i-1] - yp[i-1])*del2)/dels);
  772. }
  773.  
  774. double mylog(x) double x;
  775. {    if(x>0.) return log(x);
  776.     fprintf(stderr,"can't take log of nonpositive number");
  777.     exit(1);
  778.     return 0.;
  779. }
  780.  
  781.     /* return hyperbolic tangent of x */
  782. double tanh(x) double x;
  783. {    double t;
  784.     t = exp(x);
  785.     return (t - 1./t)/(t + 1./t);
  786. }
  787.  
  788.     /* return inverse hyperbolic tangent of x */
  789. double atanh(x) double x;
  790. {    if(x <= 0. || x >= 1.) 
  791.         {fprintf(stderr, "can't transform %f", x);
  792.         exit(1);
  793.         }
  794.     return log(sqrt((1.+x)/(1.-x)));
  795. }
  796.  
  797. /*    write out the command line */
  798. fprint_cmd(fp, format) 
  799. FILE *fp;        /* pointer to output file */
  800. char *format;     /* desired print format, possibly including program name */
  801. {    
  802. #ifdef __DESMET__
  803.     char buf[129];
  804.     int i=0;
  805.  
  806.     _lmove(129, 130, _showcs()-0x10, &buf, _showds());
  807.     buf[128]=13;
  808.     while(buf[i] != 13) i++;
  809.     buf[i]=0;
  810.     fprintf(fp, format, buf);
  811. #endif
  812.  
  813. #ifdef __TURBOC__
  814. #include <dos.h>
  815.     char buf[129];
  816.     char far *bufptr;
  817.     int i=0;
  818.  
  819.     bufptr = buf;
  820.     movedata(_psp, 130, FP_SEG(bufptr), FP_OFF(bufptr), 129);
  821.     buf[128]=13;
  822.     while(buf[i] != 13) i++;
  823.     buf[i]=0;
  824.     fprintf(fp, format, buf);
  825. #endif 
  826. }
  827.