home *** CD-ROM | disk | FTP | other *** search
/ vis-ftp.cs.umass.edu / vis-ftp.cs.umass.edu.tar / vis-ftp.cs.umass.edu / pub / CMU / cmu_files / cal / coordtrans3.c next >
C/C++ Source or Header  |  1990-12-11  |  42KB  |  1,688 lines
  1. /* coordtrans, coordmap1, coordequations, coordsave, coordrestore,
  2.  * coordfree
  3.  *
  4.  * SYNOPSIS
  5.  * coordtrans (x,y,dx,dy,n,order)
  6.  * double *x, *y, *dx, *dy;
  7.  * int n, order;
  8.  *
  9.  * coordmap1 (x,y,dx,dy)
  10.  * double x, y, *dx, *dy;
  11.  *
  12.  * coordequations (device,file)
  13.  * char *device, *file;
  14.  *
  15.  * coordsave (file)
  16.  * coordrestore (file)
  17.  * char *file;
  18.  *
  19.  * coordfree ()
  20.  *
  21.  * coordsetup (n)
  22.  * int n;
  23.  *
  24.  * coordstore (n)
  25.  * int n;
  26.  *
  27.  * coordsel (n)
  28.  * int n;
  29.  *
  30.  * DESCRIPTION
  31.  * Coordtrans takes as input two matrices of corresponding points whose
  32.  * coordinates are defined in the arrays x, y, dx, and dy, where point
  33.  * (x[i],y[i]) corresponds to (dx[i],dy[i]).  The size of the arrays is given
  34.  * by n.  A best fit transformation (x,y) => (dx,dy) is determined and
  35.  * maintained in internal data structures for use by coordmap1,
  36.  * coordequations, and coordsave.  The parameter 'order' is the greatest
  37.  * power of any term in the computed transformation equation.
  38.  *
  39.  * Coordmap1 takes as input the coordinates (x,y) of a single point and,
  40.  * implicitly the transformation computed by coordtrans.  Its returns the
  41.  * coordinates of the new point in the locations specified by dx and dy.
  42.  *
  43.  * Coordequations writes the transformation equations determined by the
  44.  * coordtrans routine into the specified file in Scribe format.  The device
  45.  * parameter is put into a Scribe @device() command at the head of the file.
  46.  * If the device specified is "CRT", the Scribe output will produce a nicely
  47.  * formatted equation on a crt.  **HACK NOTE** if the device is "none",
  48.  * coordequations will output a non-Scribe file, with no @device command and
  49.  * with exponents specified as a**b instead of a@+(b).
  50.  *
  51.  * Coordsave writes into the specified file all of the internal data
  52.  * structures computed by coordtrans that are required for coordmap1.
  53.  *
  54.  * Coordrestore reads from the specified file all the internal data
  55.  * structures previously written out by coordsave.
  56.  *
  57.  * Coordsetup, coordstore, and coordsel provide the capability for more than
  58.  * one transform, for example a transform and its inverse.  Coordsetup
  59.  * allocates space for n transforms.  Coordstore stores the current transform
  60.  * as the nth stored transform.  Coordsel retrieves the nth stored transform
  61.  * and makes it the current transform.  All other operations (coordtrans,
  62.  * coordmap1, coordequations, coordsave, and coordrestore) work on the current
  63.  * transform.  The extra storage facilities provided by coordsetup, coordstore,
  64.  * and coordsel are is invisible to users who only need one transform.
  65.  *
  66.  * Coordfree releases all storage allocated by coordtrans and coordsetup.
  67.  *
  68.  * REFERENCE
  69.  * See chapters 5 and 7 of J. G. Hayes, ed., Numerical Approximation to
  70.  * Functions and Data, London, The Athlone Press, 1970, for a theoretical
  71.  * discusson of the algorithm used in coordtrans.
  72.  *
  73.  * HISTORY
  74.  * 15-Jul-85  Chuck Thorpe (cet) at Carnegie-Mellon University
  75.  *    Added the special device "none" to coordequations to save running
  76.  *    the output through Scribe.
  77.  *
  78.  * 22-Mar-84  Duane Williams (dtw) at Carnegie-Mellon University
  79.  *    Corrected some flaws in the way coordequations output the equations.
  80.  *
  81.  * 10-Mar-83  Chuck Thorpe (cet) at Carnegie-Mellon University
  82.  *    Added multiple in-memory transforms, and routines coordsetup (),
  83.  *    coordstore (), and coordsel ().
  84.  *
  85.  * 16-Feb-83  Duane Williams (dtw) at Carnegie-Mellon University
  86.  *    Added coordsave, coordrestore, and coordfree.
  87.  *
  88.  * 26-Jan-82  Duane Williams (dtw) at Carnegie-Mellon University
  89.  *    Created.
  90.  *
  91.  */
  92.  
  93. #include <stdio.h>
  94. #include <math.h>
  95. #include "matrix.h"
  96.  
  97. typedef enum {
  98.     FALSE, TRUE
  99. } boolean;
  100.  
  101. #define hidden static
  102. #define visible
  103. #define UNDEFINED -1
  104.  
  105. #define TCOEFX(l)  tcoefx[l][0]
  106. #define TCOEFY(l)  tcoefy[l][0]
  107. #define XDEG(l)      xdeg[l][0]
  108. #define YDEG(l)      ydeg[l][0]
  109. #define BEFORE(l)  before[l][0]
  110. #define A(l)        a[l][0]
  111. #define BP(l)           bp[l][0]
  112. #define CCX(l)          ccx[l][0]
  113. #define CCY(l)          ccy[l][0]
  114.  
  115. hidden dmat tcoefx_vector;
  116. hidden dmat tcoefy_vector;
  117. hidden double **tcoefx, **tcoefy;
  118. hidden int  Order;
  119. hidden dmat bp_vector;
  120. hidden dmat c_matrix;
  121. hidden imat xdeg_vector;
  122. hidden imat ydeg_vector;
  123. hidden imat before_vector;
  124. hidden imat num_matrix;
  125. hidden double **bp, **c;
  126. hidden int **xdeg, **ydeg;
  127. hidden int **before, **num;
  128. hidden boolean tcoefx_dirty     = FALSE;
  129. hidden boolean tcoefy_dirty     = FALSE;
  130. hidden boolean bp_dirty        = FALSE;
  131. hidden boolean c_dirty        = FALSE;
  132. hidden boolean xdeg_dirty    = FALSE;
  133. hidden boolean ydeg_dirty    = FALSE;
  134. hidden boolean before_dirty    = FALSE;
  135. hidden boolean num_dirty    = FALSE;
  136.  
  137. hidden int selected = -1;
  138. hidden int scount = -1;
  139. hidden boolean *stcoefx_dirty;
  140. hidden boolean *stcoefy_dirty;
  141. hidden boolean *sbp_dirty;
  142. hidden boolean *sc_dirty;
  143. hidden boolean *sxdeg_dirty;
  144. hidden boolean *sydeg_dirty;
  145. hidden boolean *sbefore_dirty;
  146. hidden boolean *snum_dirty;
  147. hidden dmat *stcoefx_vector;
  148. hidden dmat *stcoefy_vector;
  149. hidden dmat *sc_matrix;
  150. hidden dmat *sbp_vector;
  151. hidden imat *snum_matrix;
  152. hidden imat *sxdeg_vector;
  153. hidden imat *sydeg_vector;
  154. hidden imat *sbefore_vector;
  155.  
  156. extern char *malloc();
  157. extern double pol();
  158.  
  159. /*****************************************************************
  160.  * coordtrans
  161.  *****************************************************************
  162.  */
  163. visible int
  164. coordtrans (x,y,dx,dy,n,order)
  165. double *x, *y, *dx, *dy;
  166. int n, order;
  167. {
  168.     register int    i, j, k, l;
  169.     double  sum, p;
  170.     dmat a_vector,
  171.     s_matrix,
  172.     t_matrix;
  173.     double **a, **s, **t;
  174.     int     order1, pnum, pnum1, err;
  175.  
  176.  
  177. /* save order in global */
  178.     Order = order;
  179.  
  180.     order1 = order + 1;
  181.     pnum = order1 * (order1 + 1) / 2;
  182.     pnum1 = pnum - 1;
  183.  
  184. /*----------------------------------------------------------------
  185.  * Check that we have enough input (this is probably not right)
  186.  *----------------------------------------------------------------
  187.  */
  188.     if (n < pnum) {
  189.     fprintf (stderr,
  190.         "Warning: there may not be enough data for a valid result!\n");
  191.     }
  192.  
  193. /*----------------------------------------------------------------
  194.  * Storage allocation.
  195.  *----------------------------------------------------------------
  196.  */
  197.  /* storage allocated by coordalloc stays around and is shared with
  198.     coordmap1 */
  199.  
  200.     if (coordalloc(pnum1, order) == 0) {
  201.         return 0;
  202.     }
  203.     
  204.     
  205.  /* the following storage is local to coordtrans and is freed prior to
  206.     exit */
  207.  
  208.     a_vector = newdmat (0, pnum1, 0, 0, &err);
  209.     if (err) {
  210.     printf ("\tMemory exhausted.\n");
  211.     return 0;
  212.     }
  213.     a = a_vector.el;
  214.  
  215.     s_matrix = newdmat (0, pnum1, 0, pnum1, &err);
  216.     if (err) {
  217.     printf ("\tMemory exhausted.\n");
  218.     return 0;
  219.     }
  220.     s = s_matrix.el;
  221.  
  222.     t_matrix = newdmat (0, pnum1, 0, pnum1, &err);
  223.     if (err) {
  224.     printf ("\tMemory exhausted.\n");
  225.     return 0;
  226.     }
  227.     t = t_matrix.el;
  228.  
  229.  
  230. /*----------------------------------------------------------------
  231.  * Initialization (probably unnecessary, but ...)
  232.  *----------------------------------------------------------------
  233.  */
  234.     for (i = 0; i < pnum; i++) {
  235.     A (i) = BP (i) = 0.0;
  236.     XDEG (i) = YDEG (i) = BEFORE (i) = UNDEFINED;
  237.     for (j = 0; j < pnum; j++) {
  238.         s[i][j] = t[i][j] = c[i][j] = 0.0;
  239.     }
  240.     }
  241.  
  242.     for (i = 0; i <= order; i++) {
  243.     for (j = 0; j <= order; j++) {
  244.         num[i][j] = 0;
  245.     }
  246.     }
  247.  
  248. /*----------------------------------------------------------------
  249.  * Prepare num, xdeg, ydeg, and before.
  250.  *
  251.  * Num is a square matrix of size 'order' which is meaningful only
  252.  * above the main diagonal.  It is numbered as follows (for the
  253.  * case where order = 2):
  254.  *
  255.  *    0  2  5
  256.  *    1  4
  257.  *    3
  258.  *
  259.  * XDEG(i) is the row on which the number i is located in array
  260.  * num; similarly, YDEG(i) is the column for i.
  261.  *
  262.  * BEFORE(i) is the number that is "before" number i in the num
  263.  * array.  The relation BEFORE is defined as follows:
  264.  *
  265.  *    If i = 0, then BEFORE(i) is undefined;
  266.  *    Else If XDEG(i) = 0, then BEFORE(i) is the number
  267.  *        immediately to the left of i in array num;
  268.  *         Else BEFORE(i) is the number immediately above i.
  269.  *----------------------------------------------------------------
  270.  */
  271.     for (i = 0; i <= order; i++) {
  272.     for (j = 0; (k = i + j) <= order; j++) {
  273.         num[i][j] = l = (k * (k + 1) / 2) + j;
  274.         XDEG (l) = i;
  275.         YDEG (l) = j;
  276.         BEFORE (l) = (l == 0) ? UNDEFINED
  277.         : (i == 0) ? num[0][j - 1] : num[i - 1][j];
  278.     }
  279.     }
  280.  
  281. /*-----------------------------------------------------------------------
  282.  * a[l-1] = summation of ( pol(l-1,x,y) * pol(l-1,x,y) )
  283.  *-----------------------------------------------------------------------
  284.  */
  285.     for (l = 1; l < pnum; l++) {/* main for loop */
  286.     for (sum = i = 0; i < n; i++) {
  287.         p = pol (l - 1, x[i], y[i]);
  288.         sum += p * p;
  289.     }
  290.     A (l - 1) = sum;
  291.  
  292. /*-----------------------------------------------------------------------
  293.  * s[j][l-1] = s[l-1][j] = summation of ( x * pol(j,x,y) * pol(l-1,x,y) )
  294.  *-----------------------------------------------------------------------
  295.  */
  296.     for (j = 0; j < l; j++) {
  297.  
  298.         for (sum = i = 0; i < n; i++) {
  299.         sum += x[i] * pol (j, x[i], y[i]) * pol (l - 1, x[i], y[i]);
  300.         }
  301.         s[j][l - 1] = s[l - 1][j] = sum;
  302.  
  303. /*-----------------------------------------------------------------------
  304.  * t[j][l-1] = t[l-1][j] = summation of ( y * pol(j,x,y) * pol(l-1,x,y) )
  305.  *-----------------------------------------------------------------------
  306.  */
  307.         for (sum = i = 0; i < n; i++) {
  308.         sum += y[i] * pol (j, x[i], y[i]) * pol (l - 1, x[i], y[i]);
  309.         }
  310.         t[j][l - 1] = t[l - 1][j] = sum;
  311.     }
  312.  
  313.     for (j = 0; j < l; j++) {
  314.         if (A (j) == 0) {    /* prevent zero divide */
  315.         printf ("\tWarning: A[%d] = 0.  ", j);
  316.         printf ("c[%d][%d] will be undefined.\n", l, j);
  317.         continue;
  318.         }
  319.         if (XDEG (l) == 0) {
  320.         c[l][j] = t[BEFORE (l)][j] / A (j);
  321.         }
  322.         else {
  323.         c[l][j] = s[BEFORE (l)][j] / A (j);
  324.         }
  325.     }
  326.     }                /* end main for loop */
  327.  
  328. /*----------------------------------------------------------------
  329.  * a[pnum-1] = summation of ( pol(pnum-1,x,y) * pol(pnum-1,x,y) )
  330.  *----------------------------------------------------------------
  331.  */
  332.     for (sum = i = 0; i < n; i++) {
  333.     p = pol (pnum1, x[i], y[i]);
  334.     sum += p * p;
  335.     }
  336.     A (pnum1) = sum;
  337.  
  338. /*----------------------------------------------------------------
  339.  * Determine coefficients
  340.  *----------------------------------------------------------------
  341.  */
  342.     for (l = 0; l < pnum; l++) {
  343.     if (A (l) == 0) {
  344.         printf ("\tWarning: A[%d] = 0.  ", j);
  345.         printf ("coefs[%d] will be undefined.\n", l);
  346.         continue;
  347.     }
  348.     for (sum = i = 0; i < n; i++) {
  349.         sum += dx[i] * pol (l, x[i], y[i]);
  350.     }
  351.     TCOEFX(l) = sum / A (l);
  352.  
  353.     for (sum = i = 0; i < n; i++) {
  354.         sum += dy[i] * pol (l, x[i], y[i]);
  355.     }
  356.  
  357.     TCOEFY(l) = sum / A (l);
  358.     }
  359.  
  360. /*----------------------------------------------------------------
  361.  * Free temporary storage.
  362.  *----------------------------------------------------------------
  363.  */
  364.     freemat (a_vector);
  365.     freemat (s_matrix);
  366.     freemat (t_matrix); 
  367.  
  368.     return 1;
  369. }
  370.  
  371. /*----------------------------------------------------------------
  372.  * Compute polynomial
  373.  *----------------------------------------------------------------
  374.  */
  375. hidden double
  376. pol (l,x,y)
  377. int l;
  378. double x,y;
  379. {
  380.     int i,j;
  381.     double val;
  382.  
  383.     if (l == 0) {
  384.        return (1.0);
  385.     }
  386.     
  387.     BP(0) = 1.0;
  388.  
  389.     for (i=1; i <= l; i++) {
  390.        val = ((XDEG(i) != 0) ? x : y) * BP(BEFORE(i));
  391.     
  392.        for (j=0; j < i; j++) {
  393.           val -= c[i][j] * BP(j);
  394.        }
  395.     
  396.        BP(i) = val;
  397.     }
  398.     
  399.     return ( BP(l) );
  400. }
  401.  
  402. /*****************************************************************
  403.  * coordfree -- release all internal storage associated with 
  404.  *         coordtrans
  405.  *****************************************************************
  406.  */
  407. visible
  408. coordfree ()
  409. {
  410.     int i;
  411.  
  412.     if (tcoefx_dirty) {
  413.         freemat(tcoefx_vector);    
  414.     tcoefx_dirty = FALSE;
  415.     }
  416.     if (tcoefy_dirty) {
  417.         freemat(tcoefy_vector);    
  418.     tcoefy_dirty = FALSE;
  419.     }
  420.     if (c_dirty) {
  421.         freemat(c_matrix);       
  422.     c_dirty = FALSE;
  423.     }
  424.     if (bp_dirty) {
  425.         freemat(bp_vector);
  426.     bp_dirty = FALSE;
  427.     }
  428.     if (num_dirty) {
  429.         freemat(num_matrix);
  430.     num_dirty = FALSE;
  431.     }
  432.     if (xdeg_dirty) {
  433.         freemat(xdeg_vector);
  434.     xdeg_dirty = FALSE;
  435.     }
  436.     if (ydeg_dirty) {
  437.         freemat(ydeg_vector);
  438.     ydeg_dirty = FALSE;
  439.     }
  440.     if (before_dirty) {
  441.         freemat(before_vector);
  442.     before_dirty = FALSE;
  443.     }
  444.         
  445. for (i = 0; i < scount; i++) {
  446.     if (stcoefx_dirty[i]) {
  447.     freemat (stcoefx_vector[i]);
  448.     stcoefx_dirty[i] = FALSE;
  449.     }
  450.     if (stcoefy_dirty[i]) {
  451.     freemat (stcoefy_vector[i]);
  452.     stcoefy_dirty[i] = FALSE;
  453.     }
  454.     if (sc_dirty[i]) {
  455.     freemat (sc_matrix[i]);
  456.     sc_dirty[i] = FALSE;
  457.     }
  458.     if (sbp_dirty[i]) {
  459.     freemat (sbp_vector[i]);
  460.     sbp_dirty[i] = FALSE;
  461.     }
  462.     if (snum_dirty[i]) {
  463.     freemat (snum_matrix[i]);
  464.     snum_dirty[i] = FALSE;
  465.     }
  466.     if (sxdeg_dirty[i]) {
  467.     freemat (sxdeg_vector[i]);
  468.     sxdeg_dirty[i] = FALSE;
  469.     }
  470.     if (sydeg_dirty[i]) {
  471.     freemat (sydeg_vector[i]);
  472.     sydeg_dirty[i] = FALSE;
  473.     }
  474.     if (sbefore_dirty[i]) {
  475.     freemat (sbefore_vector[i]);
  476.     sbefore_dirty[i] = FALSE;
  477.     }
  478. }
  479. if (scount != -1) {
  480.     scount = -1;
  481.     free (stcoefx_vector);
  482.     free (stcoefy_vector);
  483.     free (sbp_vector);
  484.     free (c_matrix);
  485.     free (xdeg_vector);
  486.     free (ydeg_vector);
  487.     free (before_vector);
  488.     free (num_matrix);
  489.     free (stcoefx_dirty);
  490.     free (stcoefy_dirty);
  491.     free (sbp_dirty);
  492.     free (c_dirty);
  493.     free (xdeg_dirty);
  494.     free (ydeg_dirty);
  495.     free (before_dirty);
  496.     free (num_dirty);
  497.     } 
  498. }
  499.  
  500. /*-----------------------------------------------------------------
  501.  * coordalloc -- allocate storage common to coordtrans and coordmap1
  502.  *-----------------------------------------------------------------
  503.  */
  504. hidden int
  505. coordalloc (pnum1, order)
  506. int pnum1, order;
  507. {
  508.     int err;
  509.  
  510.     if (tcoefx_dirty) {
  511.         freemat (tcoefx_vector); 
  512.     tcoefx_dirty = FALSE;
  513.     }
  514.     tcoefx_vector = newdmat (0, pnum1, 0, 0, &err);
  515.     if (err) {
  516.         printf("\tMemory exhausted.\n");
  517.     return 0;
  518.     }
  519.     tcoefx = tcoefx_vector.el;
  520.     tcoefx_dirty = TRUE;
  521.     
  522.     if (tcoefy_dirty) {
  523.         freemat (tcoefy_vector); 
  524.     tcoefy_dirty = FALSE;
  525.     }
  526.     tcoefy_vector = newdmat (0, pnum1, 0, 0, &err);
  527.     if (err) {
  528.         printf("\tMemory exhausted.\n");
  529.     return 0;
  530.     }
  531.     tcoefy = tcoefy_vector.el;
  532.     tcoefy_dirty = TRUE;
  533.     
  534.     if (c_dirty) {
  535.     freemat (c_matrix);
  536.     c_dirty = FALSE;
  537.     }
  538.     c_matrix = newdmat (0, pnum1, 0, pnum1, &err);
  539.     if (err) {
  540.     printf ("\tMemory exhausted.\n");
  541.     return 0;
  542.     }
  543.     c = c_matrix.el;
  544.     c_dirty = TRUE;
  545.  
  546.     if (bp_dirty) {
  547.     freemat (bp_vector);
  548.     bp_dirty = FALSE;
  549.     }
  550.     bp_vector = newdmat (0, pnum1, 0, 0, &err);
  551.     if (err) {
  552.     printf ("\tMemory exhausted.\n");
  553.     return 0;
  554.     }
  555.     bp = bp_vector.el;
  556.     bp_dirty = TRUE;
  557.  
  558.     if (num_dirty) {
  559.     freemat (num_matrix); 
  560.     num_dirty = FALSE;
  561.     }
  562.     num_matrix = newimat (0, order, 0, order, &err);
  563.     if (err) {
  564.     printf ("\tMemory exhausted.\n");
  565.     return 0;
  566.     }
  567.     num = num_matrix.el;
  568.     num_dirty = TRUE;
  569.  
  570.     if (xdeg_dirty) {
  571.     freemat (xdeg_vector); 
  572.     xdeg_dirty = FALSE;
  573.     }
  574.     xdeg_vector = newimat (0, pnum1, 0, 0, &err);
  575.     if (err) {
  576.     printf ("\tMemory exhausted.\n");
  577.     return 0;
  578.     }
  579.     xdeg = xdeg_vector.el;
  580.     xdeg_dirty = TRUE;
  581.  
  582.     if (ydeg_dirty) {
  583.     freemat (ydeg_vector); 
  584.     ydeg_dirty = FALSE;
  585.     }
  586.     ydeg_vector = newimat (0, pnum1, 0, 0, &err);
  587.     if (err) {
  588.     printf ("\tMemory exhausted.\n");
  589.     return 0;
  590.     }
  591.     ydeg = ydeg_vector.el;
  592.     ydeg_dirty = TRUE;
  593.  
  594.     if (before_dirty) {
  595.     freemat (before_vector); 
  596.     before_dirty = FALSE;
  597.     }
  598.     before_vector = newimat (0, pnum1, 0, 0, &err);
  599.     if (err) {
  600.     printf ("\tMemory exhausted.\n");
  601.     return 0;
  602.     }
  603.     before = before_vector.el;
  604.     before_dirty = TRUE;
  605.     
  606.     return 1;
  607. }
  608.  
  609. /*-----------------------------------------------------------------
  610.  * coordsave, coordrestore  -- write and read data structures needed by
  611.  * coordmap1 to and from a file.
  612.  *
  613.  * The following data structures are saved and restored:
  614.  *    dmat tcoefx_vector
  615.  *    dmat tcoefy_vector
  616.  *    dmat c_matrix
  617.  *    dmat bp_vector
  618.  *    imat num_matrix
  619.  *    imat xdeg_vector
  620.  *    imat ydeg_vector
  621.  *    imat before_vector
  622.  *
  623.  * During the restoration the Order variable is recomputed.
  624.  * 
  625.  * Both routines return 0 on error and 1 if successful.
  626.  *-----------------------------------------------------------------
  627.  */
  628. #define WRITE(ptr,bytes) fwrite((char *)(ptr), sizeof(char), bytes, fd)
  629. #define READ(ptr,bytes)     fread ((char *)(ptr), sizeof(char), bytes, fd)
  630. #define HEADER "!<coord>\n"
  631. #define HEADERLENGTH 9
  632. #define ASCIIHEADER "Ascii_calibration_data_--_compliments_of_jdc"
  633. #define ASCIIHEADLEN 44
  634.  
  635. visible
  636. coordsave (file)
  637. char *file;
  638. {
  639.     FILE  *fd;    /* file descriptor for access to file */
  640.     char  *header = HEADER;
  641.     
  642.     fd = fopen (file, "w");    /* create for writing */
  643.     if (fd == NULL) {
  644.         printf("coordsave: can't create file '%s'\n",file);
  645.     return 0;
  646.     }
  647.     
  648.     if (WRITE (header, HEADERLENGTH) != HEADERLENGTH) {
  649.         printf("coordsave: encountered an error writing the file header\n");
  650.         return 0;
  651.     }
  652.     
  653.     if (writedmat (fd, tcoefx_vector) == 0 ||
  654.     writedmat (fd, tcoefy_vector) == 0 ||
  655.     writedmat (fd, c_matrix)      == 0 ||
  656.     writedmat (fd, bp_vector)     == 0 ||
  657.     writeimat (fd, num_matrix)    == 0 ||
  658.     writeimat (fd, xdeg_vector)   == 0 ||
  659.     writeimat (fd, ydeg_vector)   == 0 ||
  660.     writeimat (fd, before_vector) == 0) {
  661.         printf("coordsave: encountered an error writing the file\n");
  662.     return 0;
  663.     }
  664.     
  665.     if (fclose (fd) != EOF) {
  666.         return 0;
  667.     }
  668.     
  669.     return 1;
  670. }
  671.  
  672. visible
  673. coordrestore (file)
  674. char *file;
  675. {
  676.     FILE *fd;     /* file pointer for access to file */
  677.     char  header [ HEADERLENGTH ];
  678.     
  679.     selected = -1;
  680.  
  681.     fd = fopen (file, "r");    /* open for reading */
  682.     if (fd == NULL) {
  683.         printf("coordrestore: can't open file '%s'\n",file);
  684.     return 0;
  685.     }
  686.     
  687.     if (READ(header,HEADERLENGTH) != HEADERLENGTH) {
  688.     printf("coordrestore: error in reading the file header\n");
  689.         return 0;
  690.     }
  691.     
  692.     if (strncmp (header, HEADER, HEADERLENGTH) != 0) {
  693.         printf("coordrestore: invalid file header ... restoration aborted\n");
  694.     return 0;
  695.     }
  696.             
  697.     if (readdmat (fd, &tcoefx_vector) == 0 ||
  698.     readdmat (fd, &tcoefy_vector) == 0 ||
  699.     readdmat (fd, &c_matrix)      == 0 ||
  700.     readdmat (fd, &bp_vector)     == 0 ||
  701.     readimat (fd, &num_matrix)    == 0 ||
  702.     readimat (fd, &xdeg_vector)   == 0 ||
  703.     readimat (fd, &ydeg_vector)   == 0 ||
  704.     readimat (fd, &before_vector) == 0) {
  705.     printf("coordrestore: encountered an error reading the file\n");
  706.     return 0;
  707.     }
  708.     
  709.     if (fclose (fd) == EOF) {
  710.         return 0;
  711.     }
  712.     
  713.     tcoefx = tcoefx_vector.el;
  714.     tcoefy = tcoefy_vector.el;
  715.     c      = c_matrix.el;
  716.     bp     = bp_vector.el;
  717.     num       = num_matrix.el;
  718.     xdeg   = xdeg_vector.el;
  719.     ydeg   = ydeg_vector.el;
  720.     before = before_vector.el;
  721.     tcoefx_dirty = TRUE;
  722.     tcoefy_dirty = TRUE;
  723.     c_dirty     = TRUE;
  724.     bp_dirty     = TRUE;
  725.     num_dirty     = TRUE;
  726.     xdeg_dirty     = TRUE;
  727.     ydeg_dirty   = TRUE;
  728.     before_dirty = TRUE;
  729.     
  730.     Order = num_matrix.ub1;
  731.  
  732.     return 1;
  733. }
  734.  
  735. /*-----------------------------------------------------------------
  736.  * coordasciisave, coordasciirestore  -- write and read data structures 
  737.  * needed by coordmap1 to and from an ascii file.
  738.  *
  739.  * The following data structures are saved and restored:
  740.  *    dmat tcoefx_vector
  741.  *    dmat tcoefy_vector
  742.  *    dmat c_matrix
  743.  *    dmat bp_vector
  744.  *    imat num_matrix
  745.  *    imat xdeg_vector
  746.  *    imat ydeg_vector
  747.  *    imat before_vector
  748.  *
  749.  * During the restoration the Order variable is recomputed.
  750.  * 
  751.  * Both routines return 0 on error and 1 if successful.
  752.  *-----------------------------------------------------------------
  753.  */
  754.  
  755. visible
  756. coordasciisave (file)
  757. char *file;
  758. {
  759.     FILE  *fd;    /* file descriptor for access to file */
  760.     char  *header = ASCIIHEADER;
  761.     
  762.     fd = fopen (file, "w");    /* create for writing */
  763.     if (fd == NULL) {
  764.         printf("coordsave: can't create file '%s'\n",file);
  765.     return 0;
  766.     }
  767.     
  768.     if (fprintf(fd, "%s\n", header) != ASCIIHEADLEN + 1) {
  769.         printf("coordasciisave: error in writing the file header\n");
  770.         return 0;
  771.     }
  772.     
  773.     if (fprintfdmat (fd, tcoefx_vector) == 0 ||
  774.     fprintfdmat (fd, tcoefy_vector) == 0 ||
  775.     fprintfdmat (fd, c_matrix)      == 0 ||
  776.     fprintfdmat (fd, bp_vector)     == 0 ||
  777.     fprintfimat (fd, num_matrix)    == 0 ||
  778.     fprintfimat (fd, xdeg_vector)   == 0 ||
  779.     fprintfimat (fd, ydeg_vector)   == 0 ||
  780.     fprintfimat (fd, before_vector) == 0) {
  781.         printf("coordsave: encountered an error writing the file\n");
  782.     return 0;
  783.     }
  784.     
  785.     if (fclose (fd) != EOF) {
  786.         return 0;
  787.     }
  788.  
  789.     return 1;
  790. }
  791.  
  792. visible
  793. coordasciirestore (file)
  794. char *file;
  795. {
  796.     FILE *fd;     /* file pointer for access to file */
  797.     char  asciiheader [ ASCIIHEADLEN ];
  798.     int length;
  799.     selected = -1;
  800.  
  801.     fd = fopen (file, "r");    /* open for reading */
  802.     if (fd == NULL) {
  803.         printf("coordasciirestore: can't open file '%s'\n",file);
  804.     return 0;
  805.     }
  806.  
  807.     length = fscanf(fd, "%s\n", asciiheader);
  808. /*    if(length != ASCIIHEADLEN + 1) {
  809.     printf("coordasciirestore: error in reading the ascii file header\n");
  810.         return 0;
  811.     }
  812. */    
  813.     if (strncmp (asciiheader, ASCIIHEADER, ASCIIHEADLEN) != 0) {
  814.         printf("coordasciirestore: invalid ascii file header\n");
  815.     return 0;
  816.     }
  817.             
  818.     if (fscanfdmat (fd, &tcoefx_vector) == 0 ||
  819.     fscanfdmat (fd, &tcoefy_vector) == 0 ||
  820.     fscanfdmat (fd, &c_matrix)      == 0 ||
  821.     fscanfdmat (fd, &bp_vector)     == 0 ||
  822.     fscanfimat (fd, &num_matrix)    == 0 ||
  823.     fscanfimat (fd, &xdeg_vector)   == 0 ||
  824.     fscanfimat (fd, &ydeg_vector)   == 0 ||
  825.     fscanfimat (fd, &before_vector) == 0) {
  826.     printf("coordasciirestore: encountered an error reading the file\n");
  827.     return 0;
  828.     }
  829.     
  830.     if (fclose (fd) == EOF) {
  831.         return 0;
  832.     }
  833.     
  834.     tcoefx = tcoefx_vector.el;
  835.     tcoefy = tcoefy_vector.el;
  836.     c      = c_matrix.el;
  837.     bp     = bp_vector.el;
  838.     num       = num_matrix.el;
  839.     xdeg   = xdeg_vector.el;
  840.     ydeg   = ydeg_vector.el;
  841.     before = before_vector.el;
  842.     tcoefx_dirty = TRUE;
  843.     tcoefy_dirty = TRUE;
  844.     c_dirty     = TRUE;
  845.     bp_dirty     = TRUE;
  846.     num_dirty     = TRUE;
  847.     xdeg_dirty     = TRUE;
  848.     ydeg_dirty   = TRUE;
  849.     before_dirty = TRUE;
  850.     
  851.     Order = num_matrix.ub1;
  852.  
  853.     return 1;
  854. }
  855.  
  856. hidden writedmat (fd, mat)
  857. FILE *fd;
  858. dmat  mat;
  859. {
  860.     int   size, rows;
  861.     int   nbytes = 0;
  862.     char *data;
  863.     
  864.     /* write out the matrix bounds - 4 ints */
  865.  
  866.     nbytes += WRITE (&(mat.lb1), sizeof(int));
  867.     nbytes += WRITE (&(mat.ub1), sizeof(int));
  868.     nbytes += WRITE (&(mat.lb2), sizeof(int));
  869.     nbytes += WRITE (&(mat.ub2), sizeof(int));
  870.     
  871.     /* compute the matrix size in bytes and the virtual address of the
  872.        actual data */
  873.  
  874.     rows = mat.ub1 - mat.lb1 + 1;
  875.     size = rows * (mat.ub2 - mat.lb2 + 1) * sizeof(double);
  876.     data = mat.mat_sto + rows * sizeof(double **); /* (char *) */
  877.     
  878.     /* write out the matrix data -- the matrix rows are documented to be
  879.        stored contiguously; so we can write them all at once! */
  880.     
  881.     nbytes += WRITE (data, size);
  882.     
  883.     /* check number of bytes written */
  884.     
  885.     size = 4 * sizeof(int) + size;
  886.     
  887.     if (nbytes != size) {
  888.         printf("writedmat:  should have written %d bytes -- wrote %d\n",
  889.         size, nbytes);
  890.     return 0;
  891.     }
  892.     
  893.     return 1;
  894. }
  895.  
  896. hidden writeimat (fd, mat)
  897. FILE *fd;
  898. imat  mat;
  899. {
  900.     int   size, rows;
  901.     int      nbytes = 0;
  902.     char *data;
  903.     
  904.     /* write out the matrix bounds - 4 ints */
  905.  
  906.     nbytes += WRITE (&(mat.lb1), sizeof(int));
  907.     nbytes += WRITE (&(mat.ub1), sizeof(int));
  908.     nbytes += WRITE (&(mat.lb2), sizeof(int));
  909.     nbytes += WRITE (&(mat.ub2), sizeof(int));
  910.     
  911.     /* compute the matrix size in bytes and the virtual address of the
  912.        actual data */
  913.  
  914.     rows = mat.ub1 - mat.lb1 + 1;
  915.     size = rows * (mat.ub2 - mat.lb2 + 1) * sizeof(int);
  916.     data = mat.mat_sto + rows * sizeof(int **); /* (char *) */
  917.     
  918.     /* write out the matrix data -- the matrix rows are documented to be
  919.        stored contiguously; so we can write them all at once! */
  920.     
  921.     nbytes += WRITE (data, size);
  922.  
  923.     /* check number of bytes written */
  924.     
  925.     size = 4 * sizeof(int) + size;
  926.     
  927.     if (nbytes != size) {
  928.         printf("writedmat:  should have written %d bytes -- wrote %d\n",
  929.         size, nbytes);
  930.     return 0;
  931.     }
  932.     
  933.     return 1;
  934. }
  935.  
  936. hidden readdmat (fd, mat)
  937. FILE *fd;
  938. dmat *mat;
  939. {
  940.     dmat mymat;
  941.     int  lb1, ub1, lb2, ub2;
  942.     int  nbytes = 0;
  943.     int  err, size, rows;
  944.     char *data;
  945.  
  946.     /* read the matrix bounds - 4 ints */
  947.     
  948.     nbytes += READ (&lb1, sizeof(int));
  949.     nbytes += READ (&ub1, sizeof(int));
  950.     nbytes += READ (&lb2, sizeof(int));
  951.     nbytes += READ (&ub2, sizeof(int));
  952.     
  953.     /* allocate the matrix */
  954.     
  955.     mymat = newdmat (lb1, ub1, lb2, ub2, &err);
  956.     if (err) {
  957.         printf("readdmat: can't allocate storage for matrix\n");
  958.     return 0;
  959.     }
  960.     
  961.     /* find out where to put data and how much there is */
  962.     
  963.     rows = ub1 - lb1 + 1;
  964.     size = rows * (ub2 - lb2 + 1) * sizeof(double);
  965.     data = mymat.mat_sto + rows * sizeof(double **); /* (char *) */
  966.     
  967.     /* read in the matrix data -- matrix rows are stored contiguously */
  968.  
  969.     nbytes += READ (data, size);
  970.     
  971.     /* check number of bytes read */
  972.  
  973.     size = 4 * sizeof(int) + size;
  974.     
  975.     if (nbytes != size) {
  976.         printf("readdmat:  should have read %d bytes -- read %d\n",
  977.         size, nbytes);
  978.     return 0;
  979.     }
  980.  
  981.     /* copy matrix fields to the outside world */
  982.     
  983.     mat -> lb1 = lb1;
  984.     mat -> ub1 = ub1;
  985.     mat -> lb2 = lb2;
  986.     mat -> ub2 = ub2;
  987.     mat -> mat_sto = mymat.mat_sto;
  988.     mat -> el  = mymat.el;
  989.  
  990.     return 1;
  991. }
  992.  
  993. hidden readimat (fd, mat)
  994. FILE *fd;
  995. imat *mat;
  996. {
  997.     imat mymat;
  998.     int  lb1, ub1, lb2, ub2;
  999.     int  nbytes = 0;
  1000.     int  err, size, rows;
  1001.     char *data;
  1002.  
  1003.     /* read the matrix bounds - 4 ints */
  1004.     
  1005.     nbytes += READ (&lb1, sizeof(int));
  1006.     nbytes += READ (&ub1, sizeof(int));
  1007.     nbytes += READ (&lb2, sizeof(int));
  1008.     nbytes += READ (&ub2, sizeof(int));
  1009.     
  1010.     /* allocate the matrix */
  1011.     
  1012.     mymat = newimat (lb1, ub1, lb2, ub2, &err);
  1013.     if (err) {
  1014.         printf("readimat: can't allocate storage for matrix\n");
  1015.     return 0;
  1016.     }
  1017.     
  1018.     /* find out where to put data and how much there is */
  1019.     
  1020.     rows = ub1 - lb1 + 1;
  1021.     size = rows * (ub2 - lb2 + 1) * sizeof(int);
  1022.     data = mymat.mat_sto + rows * sizeof(int **); /* (char *) */
  1023.     
  1024.     /* read in the matrix data -- matrix rows are stored contiguously */
  1025.  
  1026.     nbytes += READ (data, size);
  1027.     
  1028.     /* check number of bytes read */
  1029.  
  1030.     size = 4 * sizeof(int) + size;
  1031.     
  1032.     if (nbytes != size) {
  1033.         printf("readimat:  should have read %d bytes -- read %d\n",
  1034.         size, nbytes);
  1035.     return 0;
  1036.     }
  1037.  
  1038.     /* copy matrix fields to the outside world */
  1039.     
  1040.     mat -> lb1 = lb1;
  1041.     mat -> ub1 = ub1;
  1042.     mat -> lb2 = lb2;
  1043.     mat -> ub2 = ub2;
  1044.     mat -> mat_sto = mymat.mat_sto;
  1045.     mat -> el  = mymat.el;
  1046.  
  1047.     return 1;
  1048. }
  1049.  
  1050. hidden fprintfdmat (fd, mat)
  1051. FILE *fd;
  1052. dmat  mat;
  1053. {
  1054.     int   size, rows;
  1055.     int   nbytes = 0;
  1056.     double *data;
  1057.     register int r, c;
  1058.     
  1059.     /* write out the matrix bounds - 4 ints */
  1060.  
  1061.     nbytes += fprintf (fd, "%d ", mat.lb1);
  1062.     nbytes += fprintf (fd, "%d ", mat.ub1);
  1063.     nbytes += fprintf (fd, "%d ", mat.lb2);
  1064.     nbytes += fprintf (fd, "%d \n", mat.ub2);
  1065.     
  1066.     /* write out the matrix data */
  1067.     for (r = mat.lb1; r <= mat.ub1; r++) {
  1068.         data = &mat.el[r][mat.lb2];
  1069.         for (c = mat.lb2; c <= mat.ub2; c++) {
  1070.             fprintf(fd, "%.12f ", *data++);
  1071.         }
  1072.     }
  1073.     fprintf(fd, "\n");
  1074.     return 1;
  1075. }
  1076.  
  1077. hidden fprintfimat (fd, mat)
  1078. FILE *fd;
  1079. imat  mat;
  1080. {
  1081.     int   size, rows;
  1082.     int   nbytes = 0;
  1083.     register int *data;
  1084.     register int r, c;
  1085.     
  1086.     /* write out the matrix bounds - 4 ints */
  1087.  
  1088.     nbytes += fprintf (fd, "%d ", mat.lb1);
  1089.     nbytes += fprintf (fd, "%d ", mat.ub1);
  1090.     nbytes += fprintf (fd, "%d ", mat.lb2);
  1091.     nbytes += fprintf (fd, "%d \n", mat.ub2);
  1092.     
  1093.     /* write out the matrix data */
  1094.     for (r = mat.lb1; r <= mat.ub1; r++) {
  1095.         data = &mat.el[r][mat.lb2];
  1096.         for (c = mat.lb2; c <= mat.ub2; c++) {
  1097.             fprintf(fd, "%d ", *data++);
  1098.         }
  1099.     }
  1100.     fprintf(fd, "\n");
  1101.     return 1;
  1102. }
  1103.  
  1104. hidden fscanfdmat (fd, mat)
  1105. FILE *fd;
  1106. dmat *mat;
  1107. {
  1108.     dmat mymat;
  1109.     int  lb1, ub1, lb2, ub2, r, c;
  1110.     int  nbytes = 0;
  1111.     int  err, size, rows;
  1112.     double  *data;
  1113.  
  1114.     /* read the matrix bounds - 4 ints */
  1115.     
  1116.     fscanf(fd, "%d ", &lb1);
  1117.     fscanf(fd, "%d ", &ub1);
  1118.     fscanf(fd, "%d ", &lb2);
  1119.     fscanf(fd, "%d \n", &ub2);
  1120.     /* allocate the matrix */
  1121.     
  1122.     mymat = newdmat (lb1, ub1, lb2, ub2, &err);
  1123.     if (err) {
  1124.         printf("fscanfdmat: can't allocate storage for matrix\n");
  1125.     return 0;
  1126.     }
  1127.     
  1128.     /* find out where to put data and how much there is */
  1129.  
  1130.     for(r = lb1; r <= ub1; r++) {
  1131.         data = &mymat.el[r][lb2];
  1132.         for (c = lb2; c <= ub2; c++) {
  1133.             fscanf(fd, "%F ", data++);
  1134.     }
  1135.     }
  1136.     fscanf(fd, "\n");
  1137.  
  1138.     /* copy matrix fields to the outside world */
  1139.     
  1140.     mat -> lb1 = lb1;
  1141.     mat -> ub1 = ub1;
  1142.     mat -> lb2 = lb2;
  1143.     mat -> ub2 = ub2;
  1144.     mat -> mat_sto = mymat.mat_sto;
  1145.     mat -> el  = mymat.el;
  1146.  
  1147.     return 1;
  1148. }
  1149.  
  1150. hidden fscanfimat (fd, mat)
  1151. FILE *fd;
  1152. imat *mat;
  1153. {
  1154.     imat mymat;
  1155.     int  lb1, ub1, lb2, ub2, r, c;
  1156.     int  nbytes = 0;
  1157.     int  err, size, rows;
  1158.     int  *data;
  1159.  
  1160.     /* read the matrix bounds - 4 ints */
  1161.     
  1162.     fscanf(fd, "%d ", &lb1);
  1163.     fscanf(fd, "%d ", &ub1);
  1164.     fscanf(fd, "%d ", &lb2);
  1165.     fscanf(fd, "%d \n", &ub2);
  1166.     /* allocate the matrix */
  1167.     
  1168.     mymat = newimat (lb1, ub1, lb2, ub2, &err);
  1169.     if (err) {
  1170.         printf("fscanfimat: can't allocate storage for matrix\n");
  1171.     return 0;
  1172.     }
  1173.     
  1174.     /* find out where to put data and how much there is */
  1175.  
  1176.     for(r = lb1; r <= ub1; r++) {
  1177.         data = &mymat.el[r][lb2];
  1178.         for (c = lb2; c <= ub2; c++) {
  1179.             fscanf(fd, "%d ", data++);
  1180.     }
  1181.     }
  1182.     fscanf(fd, "\n");
  1183.  
  1184.     /* copy matrix fields to the outside world */
  1185.     
  1186.     mat -> lb1 = lb1;
  1187.     mat -> ub1 = ub1;
  1188.     mat -> lb2 = lb2;
  1189.     mat -> ub2 = ub2;
  1190.     mat -> mat_sto = mymat.mat_sto;
  1191.     mat -> el  = mymat.el;
  1192.  
  1193.     return 1;
  1194. }
  1195.  
  1196.  
  1197. /*****************************************************************
  1198.  * coordmap1 -- transform one x,y point according to the transform
  1199.  *        computed by coordtrans.
  1200.  *****************************************************************
  1201.  */
  1202. visible
  1203. coordmap1 (x,y,dx,dy)
  1204. double x, y, *dx, *dy;
  1205. {
  1206.     int i, pnum;
  1207.     double sum_x, sum_y, p;
  1208.  
  1209.     if (tcoefx_dirty == FALSE) {
  1210.         printf("coordmap1: transformation not yet defined\n");
  1211.         return;
  1212.     }
  1213.     
  1214.     pnum = (Order+1) * (Order+2) / 2;
  1215.  
  1216.     sum_x = sum_y = 0.;
  1217.  
  1218.     for (i = 0; i < pnum; i++) {
  1219.        p = pol(i,x,y);
  1220.        sum_x += TCOEFX(i) * p;
  1221.        sum_y += TCOEFY(i) * p;
  1222.     }
  1223.     
  1224.     *dx = sum_x;
  1225.     *dy = sum_y;
  1226.     
  1227.     return;
  1228. }
  1229.  
  1230. /*****************************************************************
  1231.  * transformation cooefficients
  1232.  *    used by both coordequations and coordremap (not included)
  1233.  *****************************************************************
  1234.  */
  1235. hidden dmat ccx_vector, ccy_vector;
  1236. hidden double **ccx;
  1237. hidden double **ccy;
  1238. hidden boolean ccx_dirty = FALSE;
  1239. hidden boolean ccy_dirty = FALSE;
  1240.  
  1241. /*****************************************************************
  1242.  * coordequations -- determine and print fitting polynomial
  1243.  *****************************************************************
  1244.  */
  1245. visible
  1246. coordequations (device,file)
  1247. char *device;    /* scribe device (CRT for tty output) */
  1248. char *file;
  1249. {
  1250.     char    buf[128];
  1251.     int     fd;
  1252.     int     i, j, l, pnum, pnum1, err;
  1253.     int        previous_term;
  1254.     double  sum, sum_for_x, sum_for_y;
  1255.  
  1256.     dmat pc_matrix;
  1257.     double **pc;
  1258.  
  1259.     pnum = (Order + 1) * (Order + 2) / 2;
  1260.     pnum1 = pnum - 1;
  1261.  
  1262. /*----------------------------------------------------------------
  1263.  * Allocate temporary storage
  1264.  *----------------------------------------------------------------
  1265.  */
  1266.     pc_matrix = newdmat (0, pnum1, 0, pnum1, &err);
  1267.     if (err) {
  1268.     printf ("\tMemory exhausted.\n");
  1269.     return (1);
  1270.     }
  1271.     pc = pc_matrix.el;
  1272.  
  1273.     if (ccx_dirty) {
  1274.     freemat (ccx_vector); 
  1275.     ccx_dirty = FALSE;
  1276.     }
  1277.  
  1278.     ccx_vector = newdmat (0, pnum1, 0, 0, &err);
  1279.     if (err) {
  1280.     printf ("\tMemory exhausted.\n");
  1281.     return (1);
  1282.     }
  1283.     ccx = ccx_vector.el;
  1284.     ccx_dirty = TRUE;
  1285.  
  1286.     if (ccy_dirty) {
  1287.     freemat (ccy_vector); 
  1288.     ccy_dirty = FALSE;
  1289.     }
  1290.  
  1291.     ccy_vector = newdmat (0, pnum1, 0, 0, &err);
  1292.     if (err) {
  1293.     printf ("\tMemory exhausted.\n");
  1294.     return (1);
  1295.     }
  1296.     ccy = ccy_vector.el;
  1297.     ccy_dirty = TRUE;
  1298.  
  1299.     pc[0][0] = 1.0;
  1300.  
  1301.     for (i = 1; i < pnum; i++) {
  1302.     pc[0][i] = 0;
  1303.     }
  1304.  
  1305.     for (l = 1; l < pnum; l++) {
  1306.     if (XDEG (l) == 0) {
  1307.         for (i = 0; i < pnum; i++) {
  1308.         for (sum = j = 0; j < l; j++) {
  1309.             sum += c[l][j] * pc[j][i];
  1310.         }
  1311.  
  1312.         if (YDEG (i) == 0) {
  1313.             pc[l][i] = -sum;
  1314.         }
  1315.         else
  1316.             if (YDEG (i) > 0) {
  1317.             pc[l][i] = pc[BEFORE (l)][num[XDEG (i)][YDEG (i) - 1]] - sum;
  1318.             }
  1319.         }            /* end for i */
  1320.     }
  1321.     else
  1322.         if (XDEG (l) > 0) {
  1323.         for (i = 0; i < pnum; i++) {
  1324.             for (sum = j = 0; j < l; j++) {
  1325.             sum += c[l][j] * pc[j][i];
  1326.             }
  1327.  
  1328.             if (XDEG (i) == 0) {
  1329.             pc[l][i] = -sum;
  1330.             }
  1331.             else
  1332.             if (XDEG (i) > 0) {
  1333.                 pc[l][i] = pc[BEFORE (l)][num[XDEG (i) - 1][YDEG (i)]] - sum;
  1334.             }
  1335.         }        /* end for i */
  1336.         }
  1337.     }                /* end for l */
  1338.  
  1339.     for (i = 0; i < pnum; i++) {
  1340.     for (sum_for_x = sum_for_y = j = 0; j < pnum; j++) {
  1341.         sum_for_x += TCOEFX (j) * pc[j][i];
  1342.         sum_for_y += TCOEFY (j) * pc[j][i];
  1343.     }
  1344.     CCX (i) = sum_for_x;
  1345.     CCY (i) = sum_for_y;
  1346.     }
  1347.  
  1348. /*----------------------------------------------------------------
  1349.  * create new file for output
  1350.  *----------------------------------------------------------------
  1351.  */
  1352.     fd = creat (file, 0644);
  1353.     if (fd == -1) {
  1354.     printf ("Can't open file for output.\n");
  1355.     freemat (pc_matrix); 
  1356.     return;
  1357.     }
  1358.  
  1359. /*----------------------------------------------------------------
  1360.  * write out the equations (what a mess!)
  1361.  *----------------------------------------------------------------
  1362.  */
  1363.     if (strcmp (device, "none") != 0) {
  1364.     sprintf (buf, "@device(%s)\n\n", device);
  1365.     write (fd, buf, strlen (buf));
  1366.     }
  1367.  
  1368.     sprintf (buf, "dx = ");
  1369.     write (fd, buf, strlen (buf));
  1370.  
  1371.     previous_term = 0;
  1372.  
  1373.     if (CCX (0) != 0.0) {
  1374.     sprintf (buf, "%lf", CCX (0));
  1375.     write (fd, buf, strlen (buf));
  1376.     previous_term = 1;
  1377.     }
  1378.  
  1379.     for (i = 1; i < pnum; i++) {
  1380.     if (CCX (i) != 0.0) {
  1381.         if (CCX (i) > 0.0) {
  1382.         if (previous_term) {
  1383.             sprintf (buf, " + ");
  1384.             write (fd, buf, strlen (buf));
  1385.         }
  1386.         }
  1387.         else {
  1388.         sprintf (buf, " - ");
  1389.         write (fd, buf, strlen (buf));
  1390.         }
  1391.  
  1392.         sprintf (buf, "%lf", (CCX (i) >= 0.0) ? CCX (i) : -CCX (i));
  1393.         write (fd, buf, strlen (buf));
  1394.     }
  1395.     else {
  1396.         continue;
  1397.     }
  1398.  
  1399.     if (XDEG (i) != 0) {
  1400.         sprintf (buf, " x");
  1401.         write (fd, buf, strlen (buf));
  1402.         previous_term = 1;
  1403.         
  1404.         if (XDEG (i) > 1) {
  1405.         if (strcmp (device, "none") == 0) {
  1406.             sprintf (buf, "**%d", XDEG(i));
  1407.             }
  1408.         else {
  1409.             sprintf (buf, "@+(%d)", XDEG (i));
  1410.             }
  1411.         write (fd, buf, strlen (buf));
  1412.         }
  1413.         if (YDEG (i) != 0) {
  1414.         sprintf (buf, " y");
  1415.         write (fd, buf, strlen (buf));
  1416.         
  1417.         if (YDEG (i) > 1) {
  1418.             if (strcmp (device, "none") == 0) {
  1419.             sprintf (buf, "**%d", YDEG(i));
  1420.             }
  1421.             else {
  1422.              sprintf (buf, "@+(%d)", YDEG (i));
  1423.              }
  1424.             write (fd, buf, strlen (buf));
  1425.         }
  1426.         }
  1427.     }
  1428.     else {
  1429.         if (YDEG (i) != 0) {
  1430.         sprintf (buf, " y");
  1431.         write (fd, buf, strlen (buf));
  1432.         previous_term = 1;
  1433.  
  1434.         if (YDEG (i) > 1) {
  1435.             if (strcmp (device, "none") == 0) {
  1436.             sprintf (buf, "**%d", YDEG(i));
  1437.             }
  1438.             else {
  1439.              sprintf (buf, "@+(%d)", YDEG (i));
  1440.              }
  1441.             write (fd, buf, strlen (buf));
  1442.         }
  1443.         }
  1444.     }
  1445.     }
  1446.     sprintf (buf, "\n\n");
  1447.     write (fd, buf, strlen (buf));
  1448.  
  1449.     sprintf (buf, "dy = ");
  1450.     write (fd, buf, strlen (buf));
  1451.  
  1452.     previous_term = 0;
  1453.  
  1454.     if (CCY (0) != 0.0) {
  1455.     sprintf (buf, "%lf", CCY (0));
  1456.     write (fd, buf, strlen (buf));
  1457.     previous_term = 1;
  1458.     }
  1459.  
  1460.     for (i = 1; i < pnum; i++) {
  1461.     if (CCY (i) != 0.0) {
  1462.         if (CCY (i) > 0.0) {
  1463.         if (previous_term) {
  1464.             sprintf (buf, " + ");
  1465.             write (fd, buf, strlen (buf));
  1466.         }
  1467.         }
  1468.         else {
  1469.         sprintf (buf, " - ");
  1470.         write (fd, buf, strlen (buf));
  1471.         }
  1472.  
  1473.         sprintf (buf, "%lf", (CCY (i) >= 0.0) ? CCY (i) : -CCY (i));
  1474.         write (fd, buf, strlen (buf));
  1475.     }
  1476.     else {
  1477.         continue;
  1478.     }
  1479.  
  1480.     if (XDEG (i) != 0) {
  1481.         sprintf (buf, " x");
  1482.         write (fd, buf, strlen (buf));
  1483.         previous_term = 1;
  1484.         
  1485.         if (XDEG (i) > 1) {
  1486.         if (strcmp (device, "none") == 0) {
  1487.             sprintf (buf, "**%d", XDEG(i));
  1488.             }
  1489.         else {
  1490.             sprintf (buf, "@+(%d)", XDEG (i));
  1491.             }
  1492.         write (fd, buf, strlen (buf));
  1493.         }
  1494.         if (YDEG (i) != 0) {
  1495.         sprintf (buf, " y");
  1496.         write (fd, buf, strlen (buf));
  1497.  
  1498.         if (YDEG (i) > 1) {
  1499.             if (strcmp (device, "none") == 0) {
  1500.             sprintf (buf, "**%d", YDEG(i));
  1501.             }
  1502.             else {
  1503.             sprintf (buf, "@+(%d)", YDEG (i));
  1504.             }
  1505.             write (fd, buf, strlen (buf));
  1506.         }
  1507.         }
  1508.     }
  1509.     else
  1510.         if (YDEG (i) != 0) {
  1511.         sprintf (buf, " y");
  1512.         write (fd, buf, strlen (buf));
  1513.         previous_term = 1;
  1514.         
  1515.         if (YDEG (i) > 1) {
  1516.             if (strcmp (device, "none") == 0) {
  1517.             sprintf (buf, "**%d", YDEG(i));
  1518.             }
  1519.             else {
  1520.             sprintf (buf, "@+(%d)", YDEG (i));
  1521.             }
  1522.             write (fd, buf, strlen (buf));
  1523.         }
  1524.         }
  1525.     }
  1526.     sprintf (buf, "\n");
  1527.     write (fd, buf, strlen (buf));
  1528.  
  1529.     close (fd);
  1530.     freemat (pc_matrix); 
  1531. }
  1532. /******************************************************************************
  1533.  *  coordsetup -- allocate space for multiple transforms
  1534.  ******************************************************************************
  1535.  */
  1536. visible
  1537. coordsetup (n)
  1538. int n;
  1539. {
  1540. int i;
  1541. if (scount != -1) {
  1542.     printf ("coordsetup: space has already been allocated\n");
  1543.     return 0;
  1544.     }
  1545. if (n < 0) {
  1546.     printf ("coordsetup: can't allocate %d transforms\n", n);
  1547.     return 0;
  1548.     }
  1549. scount = n;
  1550. stcoefx_vector = (dmat *) malloc (n * sizeof (dmat));
  1551. stcoefy_vector = (dmat *) malloc (n * sizeof (dmat));
  1552. sc_matrix = (dmat *) malloc (n * sizeof (dmat));
  1553. sbp_vector = (dmat *) malloc (n * sizeof (dmat));
  1554. snum_matrix = (imat *) malloc (n * sizeof (imat));
  1555. sydeg_vector = (imat *) malloc (n * sizeof (imat));
  1556. sxdeg_vector = (imat *) malloc (n * sizeof (imat));
  1557. sbefore_vector = (imat *) malloc (n * sizeof (imat));
  1558.  
  1559. stcoefx_dirty = (boolean *) malloc (n * sizeof (boolean));
  1560. stcoefy_dirty = (boolean *) malloc (n * sizeof (boolean));
  1561. sc_dirty = (boolean *) malloc (n * sizeof (boolean));
  1562. sbp_dirty = (boolean *) malloc (n * sizeof (boolean));
  1563. snum_dirty = (boolean *) malloc (n * sizeof (boolean));
  1564. sydeg_dirty = (boolean *) malloc (n * sizeof (boolean));
  1565. sxdeg_dirty = (boolean *) malloc (n * sizeof (boolean));
  1566. sbefore_dirty = (boolean *) malloc (n * sizeof (boolean));
  1567. for (i = 0; i < n; i++)
  1568.     stcoefx_dirty[i] = stcoefy_dirty[i] = sc_dirty[i] = sbp_dirty[i] = 
  1569.     sydeg_dirty[i] = sxdeg_dirty[i] = sbefore_dirty[i] = snum_dirty[i] = FALSE;
  1570. return 1;
  1571. }
  1572. /******************************************************************************
  1573.  *  coordstore -- store the current transform
  1574.  ******************************************************************************
  1575.  */
  1576. visible
  1577. coordstore (n)
  1578. int n;
  1579. {
  1580. if (scount == -1) {
  1581.     printf ("coordstore: no space allocated\n");
  1582.     return 0;
  1583.     }
  1584.  
  1585. if ((n < 0) || (n > scount)) {
  1586.     printf ("coordstore: don't have %d records allocated\n", n);
  1587.     return 0;
  1588.     }
  1589.  
  1590. if (stcoefx_dirty[n]) freemat (stcoefx_vector[n]);
  1591. stcoefx_vector[n]  = tcoefx_vector;
  1592. tcoefx_dirty = FALSE;
  1593. stcoefx_dirty[n] = TRUE;
  1594.  
  1595. if (stcoefy_dirty[n]) freemat (stcoefy_vector[n]); 
  1596. stcoefy_vector[n]  = tcoefy_vector;
  1597. tcoefy_dirty = FALSE;
  1598. stcoefy_dirty[n] = TRUE;
  1599.  
  1600. if (sbp_dirty[n]) freemat (sbp_vector[n]);
  1601. sbp_vector[n] = bp_vector;
  1602. bp_dirty = FALSE;
  1603. sbp_dirty[n] = TRUE;
  1604.  
  1605. if (sc_dirty[n]) freemat (sc_matrix[n]);
  1606. sc_matrix[n] = c_matrix;
  1607. c_dirty = FALSE;
  1608. sc_dirty[n] = TRUE;
  1609.  
  1610. if (sxdeg_dirty[n]) freemat (sxdeg_vector[n]);
  1611. sxdeg_vector[n] = xdeg_vector;
  1612. xdeg_dirty = FALSE;
  1613. sxdeg_dirty[n] = TRUE;
  1614.  
  1615. if (sydeg_dirty[n]) freemat (sydeg_vector[n]);
  1616. sydeg_vector[n] = ydeg_vector;
  1617. ydeg_dirty = FALSE;
  1618. sydeg_dirty[n] = TRUE;
  1619.  
  1620. if (sbefore_dirty[n]) freemat (sbefore_vector[n]);
  1621. sbefore_vector[n] = before_vector;
  1622. before_dirty = FALSE;
  1623. sbefore_dirty[n] = TRUE;
  1624.  
  1625. if (snum_dirty[n]) freemat (snum_matrix[n]);
  1626. snum_matrix[n] = num_matrix;
  1627. num_dirty = FALSE;
  1628. snum_dirty[n] = TRUE;
  1629.  
  1630. return 1;
  1631. }
  1632.  
  1633. /******************************************************************************
  1634.  *  coordsel -- retrieve transform from storage
  1635.  ******************************************************************************
  1636.  */
  1637. visible
  1638. coordsel (n)
  1639. int n;
  1640. {
  1641. if (scount == -1) {
  1642.     printf ("coordsel: no transforms stored\n");
  1643.     return 0;
  1644.     }
  1645.  
  1646. if ((n < 0) || (n > scount)) {
  1647.     printf ("coordsel: don't have %d transforms\n", n);
  1648.     return 0;
  1649.     }
  1650.  
  1651. if (stcoefx_dirty[n] == FALSE) {
  1652.     printf ("coordsel: transform %d not stored\n", n);
  1653.     return 0;
  1654.     }
  1655.  
  1656. if (n == scount) return 1;
  1657.  
  1658. tcoefx_vector = stcoefx_vector[n];
  1659. tcoefy_vector = stcoefy_vector[n];
  1660. c_matrix = sc_matrix[n];
  1661. bp_vector = sbp_vector[n];
  1662. num_matrix = snum_matrix[n];
  1663. xdeg_vector = sxdeg_vector[n];
  1664. ydeg_vector = sydeg_vector[n];
  1665. before_vector = sbefore_vector[n];
  1666.  
  1667. tcoefx = tcoefx_vector.el;
  1668. tcoefy = tcoefy_vector.el;
  1669. c      = c_matrix.el;
  1670. bp     = bp_vector.el;
  1671. num       = num_matrix.el;
  1672. xdeg   = xdeg_vector.el;
  1673. ydeg   = ydeg_vector.el;
  1674. before = before_vector.el;
  1675. tcoefx_dirty = TRUE;
  1676. tcoefy_dirty = TRUE;
  1677. c_dirty     = TRUE;
  1678. bp_dirty     = TRUE;
  1679. num_dirty     = TRUE;
  1680. xdeg_dirty     = TRUE;
  1681. ydeg_dirty   = TRUE;
  1682. before_dirty = TRUE;
  1683.  
  1684. Order = num_matrix.ub1;
  1685.  
  1686. return 1;
  1687. }
  1688.