home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / pshistogram.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-03-13  |  15.6 KB  |  611 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)pshistogram.c    2.22  3/13/95
  3.  *
  4.  *    Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
  5.  *    See README file for copying and redistribution conditions.
  6.  *--------------------------------------------------------------------*/
  7. /*
  8.  * pshistogram.c -- a program for plotting histograms
  9.  *
  10.  *
  11.  * Author:    Walter H. F. Smith
  12.  * Date:    15 February, 1988
  13.  *
  14.  * Updated to v2.0 5-May-1991 Paul Wessel
  15.  * Updated to v3.0 1-Jan-1995 Paul Wessel
  16.  *
  17.  */
  18.  
  19. #include "gmt.h"
  20.  
  21. #define    COUNTS 0
  22. #define FREQ_PCT 1
  23. #define LOG_COUNTS 2
  24. #define LOG_FREQ_PCT 3
  25.  
  26. double    yy0, yy1, xmin, xmax;
  27.  
  28. float    *x;
  29. int    *boxh;
  30.  
  31. int    n = 0;
  32. int    n_boxes = 0;
  33. int    n_counted = 0;
  34. double    box_width = 0.0;
  35. double    west = 0.0, east = 0.0, south = 0.0, north = 0.0;
  36. float    stats[6];
  37.  
  38. FILE    *fp_in = NULL;
  39.  
  40. char    buffer[512];
  41.  
  42. int    center_box = FALSE, cumulative = FALSE, draw_outline = FALSE;
  43. int    second = FALSE, hist_type = COUNTS;
  44.  
  45. struct PEN pen;
  46. struct FILL fill;
  47.     
  48.  
  49. main (argc, argv)
  50. int    argc;
  51. char    **argv;
  52. {
  53.     int i;
  54.     int inquire = FALSE, error = FALSE, automatic = FALSE, stairs = FALSE;
  55.     double tmp;
  56.     
  57.     argc = gmt_begin (argc, argv);
  58.  
  59.     gmt_init_pen (&pen, 1);
  60.     gmt_init_fill (&fill, 0, 0, 0);
  61.     
  62.     for (i = 1; i < argc; i++) {
  63.         if (argv[i][0] == '-') {
  64.             switch (argv[i][1]) {
  65.         
  66.                 /* Common parameters */
  67.             
  68.                 case 'B':
  69.                 case 'H':
  70.                 case 'J':
  71.                 case 'K':
  72.                 case 'O':
  73.                 case 'P':
  74.                 case 'R':
  75.                 case 'U':
  76.                 case 'V':
  77.                 case 'X':
  78.                 case 'x':
  79.                 case 'Y':
  80.                 case 'y':
  81.                 case 'c':
  82.                 case '\0':
  83.                     error += get_common_args (argv[i], &west, &east, &south, &north);
  84.                     break;
  85.  
  86.                 /* Supplemental parameters */
  87.             
  88.                 case '2':
  89.                     second = TRUE;
  90.                     break;
  91.                 case 'C':
  92.                     center_box = TRUE;
  93.                     break;
  94.                 case 'E':
  95.                     sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation);
  96.                     break;
  97.                 case 'G':
  98.                     if (gmt_getfill (&argv[i][2], &fill)) {
  99.                         gmt_fill_syntax ('G');
  100.                         error++;
  101.                     }
  102.                     break;
  103.                 case 'L':        /* Set line attributes */
  104.                     if (gmt_getpen (&argv[i][2], &pen)) {
  105.                         gmt_pen_syntax ('L');
  106.                         error++;
  107.                     }
  108.                     draw_outline = TRUE;
  109.                     break;
  110.                 case 'I':
  111.                     inquire = TRUE;
  112.                     break;
  113.                 case 'Q':
  114.                     cumulative = TRUE;
  115.                     break;
  116.                 case 'S':
  117.                     stairs = TRUE;
  118.                     break;
  119.                 case 'W':
  120.                     box_width = atof (&argv[i][2]);
  121.                     break;
  122.                 case 'Z':
  123.                     hist_type = atoi (&argv[i][2]);
  124.                     break;
  125.                 default:
  126.                     error = TRUE;
  127.                     gmt_default_error (argv[i][1]);
  128.                     break;
  129.             }
  130.         }
  131.         else {
  132.             if ((fp_in = fopen(argv[i], "r")) == NULL) {
  133.                 fprintf (stderr, "pshistogram: Cannot open file %s\n", argv[i]);
  134.                 exit(-1);
  135.             }
  136.         }
  137.     }
  138.  
  139.     if (argc == 1 || gmt_quick) {
  140.         fprintf (stderr,"pshistogram %s - Calculate and plot histograms\n\n", GMT_VERSION);
  141.         fprintf (stderr, "usage: pshistogram [file] -Jx<scale> -W [-2] [-B<tickinfo>] [-C] [-Eaz/el] [-G<fill>]\n");
  142.         fprintf (stderr, "    [-H] [-I] [-K] [-L<pen>] [-O] [-P] [-Q] [-Rw/e/s/n] [-S] [-U[text]] [-V]\n");
  143.         fprintf (stderr, "    [-X<x_shift>] [-Y<y_shift>] [-Z[type]] [-c<ncopies>]\n\n");
  144.         
  145.         if (gmt_quick) exit(-1);
  146.         fprintf (stderr, "    -Jx for linear projection.  Scale in %s/units\n", gmt_unit_names[gmtdefs.measure_unit]);
  147.         fprintf (stderr, "        Use / to specify separate x/y scaling.\n");
  148.         fprintf (stderr, "        If -JX is used then give axes lengths rather than scales\n");
  149.  
  150.         fprintf (stderr, "    -W sets the bin width\n");
  151.         fprintf (stderr, "\n\tOPTIONS:\n");
  152.         fprintf (stderr, "    -2 read second rather than first column\n");
  153.         explain_option ('b');
  154.         fprintf (stderr, "    -C will center the bins\n");
  155.         fprintf (stderr, "    -E set azimuth and elevation of viewpoint for 3-D perspective [180/90]\n");
  156.         fprintf (stderr, "    -G specify Grayshade for columns.\n");
  157.         fprintf (stderr, "       Grayshade range is black (0) through white (255).  For color fill, use\n");
  158.         fprintf (stderr, "       -Gr/g/b, where r, g, and b each range from 0 to 255\n");
  159.         explain_option ('H');
  160.         fprintf (stderr, "    -I will inquire about min/max x and y.  No plotting is done\n");
  161.         explain_option ('K');
  162.         fprintf (stderr, "    -L <pen> is penwidth. Specify r/g/b for colored line.\n");
  163.         fprintf (stderr, "       Append O for dOtted line, A for dAshed [Default = 1, black, solid]\n");
  164.         explain_option ('O');
  165.         explain_option ('P');
  166.         fprintf (stderr, "    -Q plot a cumulative histogram\n");
  167.         explain_option ('r');
  168.         fprintf (stderr, "       If neither -R nor -I are set, w/e/s/n will be based on input data\n");
  169.         fprintf (stderr, "    -S Draws a stairs-step diagram instead of histogram.\n");
  170.         explain_option ('U');
  171.         explain_option ('V');
  172.         explain_option ('X');
  173.         fprintf (stderr, "    -Z to choose type of vertical axis.  Select from\n");
  174.         fprintf (stderr, "       0 - Counts [Default]\n");
  175.         fprintf (stderr, "       1 - Frequency percent\n");
  176.         fprintf (stderr, "       2 - Log10 (counts)\n");
  177.         fprintf (stderr, "       3 - Log10 (frequency percent)\n");
  178.         explain_option ('c');
  179.         explain_option ('.');
  180.  
  181.         exit(-1);
  182.     }
  183.     
  184.     if (project_info.projection < 0 || project_info.projection > 9) {
  185.         fprintf (stderr, "%s: GMT SYNTAX ERROR -J option:  Only linear projection supported.\n", gmt_program);
  186.         error++;
  187.     }
  188.     if (z_project.view_azimuth > 360.0 || z_project.view_elevation <= 0.0 || z_project.view_elevation > 90.0) {
  189.         fprintf (stderr, "%s: GMT SYNTAX ERROR -E option:  Enter azimuth in 0-360 range, elevation in 0-90 range\n", gmt_program);
  190.         error++;
  191.     }
  192.     if (hist_type < COUNTS || hist_type > LOG_FREQ_PCT) {
  193.         fprintf (stderr, "%s: GMT SYNTAX ERROR -Z option:  histogram type must be 0, 1, 2, or 3\n", gmt_program);
  194.         error++;
  195.     }
  196.     if (box_width <= 0.0) {
  197.         fprintf (stderr, "%s: GMT SYNTAX ERROR -W option:  bin width must be nonzero\n", gmt_program);
  198.         error++;
  199.     }
  200.     
  201.     if (error) exit (-1);
  202.  
  203.     if (inquire) gmtdefs.verbose = TRUE;
  204.     if (!inquire && (west == east || south == north)) automatic = TRUE;
  205.     
  206.     if (fp_in == NULL) fp_in = stdin;
  207.  
  208.     if ( read_data() ) {
  209.         fprintf (stderr, "pshistogram: Fatal error, read only 0 points.\n");
  210.         exit (-1);
  211.     }
  212.     
  213.     if (gmtdefs.verbose) {
  214.         fprintf (stderr, "pshistogram: Extreme values of the data :\t%lg\t%lg\n", x[0], x[n-1]);
  215.         fprintf (stderr, "pshistogram: Locations: L2, L1, LMS; Scales: L2, L1, LMS\t%lg\t%lg\t%lg\t%lg\t%lg\t%lg\n",
  216.             stats[0], stats[1], stats[2], stats[3], stats[4], stats[5]);
  217.     }
  218.     
  219.     if (west == east) {    /* Set automatic x range [ and tickmarks] */
  220.         if (frame_info.anot_int[0] == 0.0) {
  221.             tmp = pow (10.0, floor (d_log10 (xmax-xmin)));
  222.             if (((xmax-xmin) / tmp) < 3.0) tmp *= 0.5;
  223.         }
  224.         else
  225.             tmp = frame_info.anot_int[0];
  226.         west = floor (xmin / tmp) * tmp;
  227.         east = ceil (xmax / tmp) * tmp;
  228.         if (frame_info.anot_int[0] == 0.0) {
  229.             frame_info.anot_int[0] = frame_info.frame_int[0] = tmp;
  230.             frame_info.plot = TRUE;
  231.         }
  232.     }
  233.     
  234.     if ( fill_boxes () ) {
  235.         fprintf (stderr, "pshistogram: Fatal error during box fill.\n");
  236.         exit (-1);
  237.     }
  238.     
  239.     if (gmtdefs.verbose) fprintf (stderr, "\npshistogram: min/max values are :\t%lg\t%lg\t%lg\t%lg\n", xmin, xmax, yy0, yy1);
  240.     
  241.     if (inquire) {    /* Only info requested, quit before plotting */
  242.         free((char *) x);
  243.         free((char *) boxh);
  244.         exit (0);
  245.     }
  246.     
  247.     if (automatic) {    /* Set up s/n based on 'clever' rounding up of the minmax values */
  248.         project_info.region = TRUE;
  249.         south = 0.0;
  250.         if (frame_info.anot_int[1] == 0.0) {
  251.             tmp = pow (10.0, floor (d_log10 (yy1)));
  252.             if ((yy1 / tmp) < 3.0) tmp *= 0.5;
  253.         }
  254.         else
  255.             tmp = frame_info.anot_int[1];
  256.         north = ceil (yy1 / tmp) * tmp;
  257.         if (frame_info.anot_int[1] == 0.0) {    /* Tickmarks not set */
  258.             frame_info.anot_int[1] = frame_info.frame_int[1] = tmp;
  259.             frame_info.plot = TRUE;
  260.         }
  261.         if (project_info.pars[0] == 0.0) {    /* Must give default xscale */
  262.             project_info.pars[0] = gmtdefs.x_axis_length / (east - west);
  263.             project_info.projection = LINEAR;
  264.         }
  265.         if (project_info.pars[1] == 0.0) {    /* Must give default yscale */
  266.             project_info.pars[1] = gmtdefs.y_axis_length / north;
  267.             project_info.projection = LINEAR;
  268.         }
  269.     }
  270.     
  271.     if (automatic && gmtdefs.verbose) fprintf (stderr, "pshistogram: Use w/e/s/n = %lg/%lg/%lg/%lg and x-tic/y-tick = %lg/%lg\n",
  272.         west, east, south, north, frame_info.anot_int[0], frame_info.anot_int[1]);
  273.         
  274.     map_setup (west, east, south, north);
  275.  
  276.     ps_plotinit (CNULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin,
  277.         gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies,
  278.         gmtdefs.dpi, gmtdefs.measure_unit, gmtdefs.paper_width, gmtdefs.page_rgb, gmt_epsinfo (argv[0]));
  279.     echo_command (argc, argv);
  280.     if (gmtdefs.unix_time) timestamp (argc, argv);
  281.     
  282.     if (project_info.three_D) ps_transrotate (-z_project.xmin, -z_project.ymin, 0.0);
  283.  
  284.     map_clip_on (-1, -1, -1, 3);
  285.     if ( plot_boxes (stairs) ) {
  286.         fprintf (stderr, "pshistogram: Fatal error during box plotting.\n");
  287.         exit (-1);
  288.     }
  289.     
  290.     map_clip_off();
  291.     
  292.     if (frame_info.plot) {
  293.         ps_setpaint (gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  294.         map_basemap ();
  295.     }
  296.     ps_setpaint (0, 0, 0);
  297.  
  298.     if (project_info.three_D) ps_rotatetrans (z_project.xmin, z_project.ymin, 0.0);
  299.     ps_plotend (gmtdefs.last_page);
  300.     
  301.     free((char *) x);
  302.     free((char *) boxh);
  303.     
  304.     gmt_end (argc, argv);
  305. }
  306.  
  307. int    read_data () {
  308.     int i, n_alloc = GMT_CHUNK;
  309.     float tmp;
  310.     
  311.     x = (float *) memory (CNULL, n_alloc , sizeof (float), "pshistogram");
  312.     
  313.     n = 0;
  314.     xmin = MAXDOUBLE;    xmax = -MAXDOUBLE;
  315.     
  316.     if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) fgets (buffer, 512, fp_in);
  317.     
  318.     while (fgets (buffer, 512, fp_in)) {
  319.     
  320.         /* Check for a blank line  */
  321.         
  322.         if (buffer[0] == '\0') continue;
  323.  
  324.         if (sscanf(buffer, "%f %f", &x[n], &tmp) ) {
  325.             if (second) x[n] = tmp;
  326.             n++;
  327.         }
  328.         if (n == n_alloc) {
  329.             n_alloc += GMT_CHUNK;
  330.             x = (float *) memory ((char *)x, n_alloc , sizeof (float), "pshistogram");
  331.         }
  332.         xmin = MIN (xmin, x[n-1]);
  333.         xmax = MAX (xmax, x[n-1]);
  334.     }
  335.     
  336.     if (gmtdefs.verbose) fprintf (stderr,"N points read:\t%d\n", n);
  337.  
  338.     x = (float *) memory ((char *)x, n , sizeof (float), "pshistogram");
  339.     if (fp_in != stdin) fclose (fp_in);
  340.  
  341.     get_loc_scl (x, n, stats);
  342.     
  343.     return( !(n) );
  344. }
  345.  
  346. int    fill_boxes () {
  347.  
  348.     double    add_half = 0.0;
  349.     
  350.     int    b0, b1, i, ibox, count_sum;
  351.  
  352.     n_boxes = ceil( (east - west) / box_width);
  353.     
  354.     if (center_box) {
  355.         n_boxes++;
  356.         add_half = 0.5;
  357.     }
  358.     
  359.     if (n_boxes <= 0) return (-1);
  360.  
  361.     boxh = (int *) memory (CNULL, n_boxes , sizeof (int), "pshistogram");
  362.  
  363.     n_counted = 0;
  364.     
  365.     /* First fill boxes with counts  */
  366.     
  367.     for (i = 0; i < n; i++) {
  368.         ibox = floor( ( (x[i] - west) / box_width) + add_half);
  369.         if (ibox < 0 || ibox >= n_boxes) continue;
  370.         boxh[ibox] ++;
  371.         n_counted++;
  372.     }
  373.  
  374.     if (cumulative) {
  375.         count_sum = 0;
  376.         for (ibox = 0; ibox < n_boxes; ibox++) {
  377.             count_sum += boxh[ibox];
  378.             boxh[ibox] = count_sum;
  379.         }
  380.         b0 = 0;
  381.         b1 = count_sum;
  382.     }
  383.     else {
  384.         b0 = n_counted;
  385.         b1 = 0;
  386.         for (ibox = 0; ibox < n_boxes; ibox++) {
  387.             if (b0 > boxh[ibox]) b0 = boxh[ibox];
  388.             if (b1 < boxh[ibox]) b1 = boxh[ibox];
  389.         }
  390.     }
  391.  
  392.     /* Now find out what the min max y will be  */
  393.     
  394.     if (b0) {
  395.         if (hist_type == LOG_COUNTS)
  396.             yy0 = d_log1p( (double)b0);
  397.         else if (hist_type == FREQ_PCT)
  398.             yy0 = (100.0 * b0) / n_counted;
  399.         else if (hist_type == LOG_FREQ_PCT)
  400.             yy0 = d_log1p( 100.0 * b0 / n_counted );
  401.         else
  402.             yy0 = b0;
  403.     }
  404.     else {
  405.         yy0 = 0;
  406.     }
  407.     if (b1) {
  408.         if (hist_type == LOG_COUNTS)
  409.             yy1 = d_log1p( (double)b1);
  410.         else if (hist_type == FREQ_PCT)
  411.             yy1 = (100.0 * b1) / n_counted;
  412.         else if (hist_type == LOG_FREQ_PCT)
  413.             yy1 = d_log1p( 100.0 * b1 / n_counted );
  414.         else
  415.             yy1 = b1;
  416.     }
  417.     else {
  418.         yy1 = 0;
  419.     }
  420.     return (0);
  421. }
  422.  
  423. int    plot_boxes (stairs)
  424. BOOLEAN stairs;
  425. {
  426.  
  427.     int    i, ibox, first = FALSE;
  428.     
  429.     double    x[4], y[4], xx, yy;
  430.     
  431.     if (draw_outline) {        
  432.         ps_setline (pen.width);
  433.         ps_setpaint (pen.r, pen.g, pen.b);
  434.         if (pen.texture) ps_setdash (pen.texture, pen.offset);
  435.     }
  436.     
  437.     for (ibox = 0; ibox < n_boxes; ibox++) {
  438.         if (stairs || boxh[ibox]) {
  439.             x[0] = west + ibox * box_width;
  440.             if (center_box) x[0] -= (0.5 * box_width);
  441.             x[1] = x[0] + box_width;
  442.             if (x[0] < west) x[0] = west;
  443.             if (x[1] > east) x[1] = east;
  444.             x[2] = x[1];
  445.             x[3] = x[0];
  446.             y[0] = y[1] = south;
  447.             if (hist_type == LOG_COUNTS)
  448.                 y[2] = d_log1p( (double)boxh[ibox]);
  449.             else if (hist_type == FREQ_PCT)
  450.                 y[2] = (100.0 * boxh[ibox]) / n_counted;
  451.             else if (hist_type == LOG_FREQ_PCT)
  452.                 y[2] = d_log1p( 100.0 * boxh[ibox] / n_counted );
  453.             else
  454.                 y[2] = boxh[ibox];
  455.             y[3] = y[2];
  456.             
  457.             for (i = 0; i < 4; i++) {
  458.                 geo_to_xy (x[i], y[i], &xx, &yy);
  459.                 if (project_info.three_D) xyz_to_xy (xx, yy, project_info.z_level, &xx, &yy);
  460.                 x[i] = xx;    y[i] = yy;
  461.             }
  462.             
  463.             if (stairs) {
  464.                 if (first) {
  465.                     first = FALSE;
  466.                     ps_plot (x[0], y[0], 2);
  467.                 }
  468.                 ps_plot (x[3], y[3], 3);
  469.                 ps_plot (x[2], y[2], 3);
  470.             }
  471.             else if (fill.use_pattern)
  472.                 ps_imagefill (x, y, 4, fill.pattern_no, fill.pattern, fill.inverse, fill.icon_size, draw_outline);
  473.             else
  474.                 ps_polygon (x, y, 4, fill.r, fill.g, fill.b, draw_outline);
  475.         }
  476.     }
  477.  
  478.     if (stairs) ps_plot (x[1], y[1], 3);
  479.     
  480.     if (draw_outline && pen.texture) ps_setdash (CNULL, 0);
  481.     
  482.     return(0);
  483. }
  484.  
  485. int    get_loc_scl (x, n, stats)
  486. float    x[], stats[];    /* L2, L1, LMS location, L2, L1, LMS scale  */
  487. int    n;
  488. {
  489.  
  490.     int    i, j, istop, multiplicity = 1, compx();
  491.     double    xsum = 0.0, x2sum = 0.0, mid_point_sum = 0.0, length, short_length = 1.0e+30;
  492.  
  493.     if (n < 3) return(-1);
  494.     
  495.     qsort((char *) x, n, sizeof(float), compx);
  496.  
  497.     j = n/2;
  498.     
  499.     if (n%2)
  500.         stats[1] = x[j];
  501.     else
  502.         stats[1] = 0.5 * (x[j] + x[j-1]);
  503.  
  504.     istop = n - j;
  505.     for (i = 0; i < istop; i++) {
  506.         length = x[i + j] - x[i];
  507.         if (length == short_length) {
  508.             multiplicity++;
  509.             mid_point_sum += (0.5 * (x[i + j]  + x[i]) );
  510.         }
  511.         else if (length < short_length) {
  512.             multiplicity = 1;
  513.             mid_point_sum = (0.5 * (x[i + j]  + x[i]) );
  514.             short_length = length;
  515.         }
  516.     }
  517.     if (multiplicity > 1) mid_point_sum /=  multiplicity;
  518.     stats[2] = mid_point_sum;
  519.     
  520.     getmad (x, n, &stats[1], &stats[4]);
  521.     getmad (x, n, &stats[2], &stats[5]);
  522.  
  523.     for (i = 0; i < n; i++) {
  524.         xsum += x[i];
  525.         x2sum += (x[i] * x[i]);
  526.     }
  527.     stats[0] = xsum / n;
  528.     stats[3] = sqrt( (n * x2sum - xsum * xsum) / (n * (n - 1) ) );
  529.  
  530.     return(0);
  531. }
  532.  
  533. int    getmad (x, n, location, scale)
  534. float    x[], *location, *scale;
  535. int    n;
  536. {
  537.     double    e_low, e_high, error, last_error;
  538.     int    i_low, i_high, n_dev, n_dev_stop;
  539.     
  540.     i_low = 0;
  541.     while (i_low < n && x[i_low] <= *location) i_low++;
  542.     i_low--;
  543.     
  544.     i_high = n - 1;
  545.     while (i_high >= 0 && x[i_high] >= *location) i_high--;
  546.     i_high++;
  547.  
  548.     n_dev_stop = n / 2;
  549.     error = 0.0;
  550.     n_dev = 0;
  551.  
  552.  
  553.     while (n_dev < n_dev_stop) {
  554.     
  555.         last_error = error;
  556.  
  557.         if (i_low < 0) {
  558.             error = x[i_high] - *location;
  559.             i_high++;
  560.             n_dev++;
  561.         }
  562.         
  563.         else if (i_high == n) {
  564.             error = *location - x[i_low];
  565.             i_low--;
  566.             n_dev++;
  567.         }
  568.         else {
  569.             e_low = *location - x[i_low];
  570.             e_high = x[i_high] - *location;
  571.             
  572.             if (e_low < e_high) {
  573.                 error = e_low;
  574.                 i_low--;
  575.                 n_dev++;
  576.             }
  577.             else if (e_high < e_low) {
  578.                 error = e_high;
  579.                 i_high++;
  580.                 n_dev++;
  581.             }
  582.             else {
  583.                 error = e_high;
  584.                 i_low--;
  585.                 i_high++;
  586.                 if ( !(n_dev)) n_dev--; /* First time count only once  */
  587.                 n_dev += 2;
  588.             }
  589.         }
  590.     }
  591.     
  592.     if (n%2) 
  593.         *scale = 1.4826 * error;
  594.     else
  595.         *scale = 0.7413 * (error + last_error);
  596.  
  597.     return (0);
  598. }
  599.  
  600. int    compx (p1, p2)
  601. float    *p1, *p2;
  602. {
  603.     if ( (*p1) < (*p2) )
  604.         return(-1);
  605.     else if ( (*p1) > (*p2) )
  606.         return(1);
  607.     else
  608.         return(0);
  609. }
  610.  
  611.