#define PROG_NAME "makefigs_mscale"
#define PROG_DESC "generate figures for the multiscale talk, UTFPR 2012-12-12"
#define PROG_VERS "1.0"

/* Copyright © 2012 by the State University of Campinas (UNICAMP). */
/* See the copyright, authorship, and warranty notice at end of file. */
/* Last edited on 2024-12-21 14:02:25 by stolfi */

/* Created dec/2012 by Jorge Stolfi, UNICAMP */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <values.h>

#include <interval.h>
#include <affirm.h>
#include <pswr.h>
#include <frgb.h>
#include <rn.h>
#include <jsrandom.h>
#include <jsmath.h>
#include <jsfile.h>

#include <figtools_mscale.h>

/* INTERNAL PROTOTYPES */

int main(int argc, char **argv);

void figs_poisson(PSStream **psP, plot_options_t *o);
  /* Outputs figures reated to Poisson problem. */
  
void figs_poisson_metbar(PSStream **psP, plot_options_t *o, int nst);
  /* Metal bar with ice baths at ends and heat sources/sinks along length. */

void figs_poisson_disbar(PSStream **psP, plot_options_t *o, int nst);
  /* Discrete version of Poisson heat problem. */
  
void figs_poisson_bar_heaters(PSStream **psP, plot_options_t *o, int nst, double xmin, double xmax, double ymin);
  /* Used by {figs_poisson_metbar,figs_poisson_disbar}. */

void figs_poisson_graphs
  ( PSStream **psP, 
    plot_options_t *o, 
    int nst, 
    int level, 
    int freq, 
    int itplot[], 
    int nplot
  );
  /* Generates figures of one-dimensional heat equilibrium problem. 
    Uses the discrete formulation with {nst} intervals and {nst+1}
    equally spaced sample points.  Plots initial state, and also after
    iteration {itplot[k]} for {k=0..nplot-1}.
    
    If {freq} is zero, starts with a random initial guess. If {freq} is
    positive, starts with a guess that is the true solution plus a
    sinusoid with {freq} half-periods.
    
    The file names are "{tag}-v{LL}-f{FFF}-i{IIIIII}.eps" where {tag}
    is "pwr" (input-output heat power), "tru" (true solution), "cur"
    (current approximation to the solution), or "err" (error, that is,
    current minus true); {LL} is the 2-digit {level}; {FFF} is the 3-digit
    {freq}; and {IIIIII} is the 6-digit iteration number {it}.
    The {level} is used only for the file name. */
  
void figs_poisson_set_power(double P[], int nst, int level);
void figs_poisson_get_true (double P[], double S[], int nst, double kappa);
void figs_poisson_graphs_get_guess(double P[], double kappa, double S[], double T[], int nst, int freq);
void figs_poisson_graphs_iterate(double T[], double P[], int nst, double kappa);
void figs_poisson_graphs_compute_error(double T[], double S[], double E[], int nst);  
void figs_poisson_graphs_downscale_temp(double U[], int nst, double UR[], int nstR);
void figs_poisson_graphs_upscale_temp(double UR[], int nstR, double U[], int nst);
void figs_poisson_downscale_power(double P[], int nst, double PR[], int nstR);

void figs_poisson_graphs_plot
  ( PSStream **psP, 
    plot_options_t *o, 
    double Y[], 
    double Ymin,
    double Ymax,
    int nst, 
    frgb_t dotcolor,
    char *tag, 
    int level, 
    int freq, 
    int it
  );
  /* Plots {(i,Y[i])} for {i} in {0..nst} as stick-and-ball. */

void parse_options(int argc, char **argv, plot_options_t **op);
void get_arg_double(double *varp, int *argnp, int argc, char **argv, char *usage);
void get_arg_float(float *varp, int *argnp, int argc, char **argv, char *usage);
void get_arg_string(char **varp, int *argnp, int argc, char **argv, char *usage);
void arg_error(char *msg, char *arg, char *pname, char *usage);

/* IMPLEMENTATIONS */

int main(int argc, char **argv)
  {   
    PSStream *ps = NULL;

    plot_options_t *o;
    parse_options(argc, argv, &o);
    start_album(&ps,o);
    
    figs_poisson(&ps, o);
    
    fprintf(stderr, "OK so far...\n");
    finish_album(&ps, o);
    return(0);
  }

void figs_poisson(PSStream **psP, plot_options_t *o)
  {
    int nst = 32; /* Number of steps along bar. */
    
    figs_poisson_metbar(psP, o, nst);
    figs_poisson_disbar(psP, o, nst);

    /* Plots of fuctions and errors: */
    int maxlevel = 2;
    int level;
    int rfact = 1;
    for (level = 0; level <= maxlevel; level++)
      { int itplot[20] = { 0,1,2,3,4,5,6,7,8,9,10,20,30,40,50, 100,200,400,800,1600 };
        int nplot = 20;
        int rnst = nst/rfact;
        figs_poisson_graphs(psP, o, rnst, level,   0, itplot, nplot);
        figs_poisson_graphs(psP, o, rnst, level, 999, itplot, nplot);
        if (rnst >=  2) { figs_poisson_graphs(psP, o, rnst, level,   1, itplot, nplot); }
        if (rnst >=  4) { figs_poisson_graphs(psP, o, rnst, level,   2, itplot, nplot); }
        if (rnst >=  8) { figs_poisson_graphs(psP, o, rnst, level,   4, itplot, nplot); }
        if (rnst >= 16) { figs_poisson_graphs(psP, o, rnst, level,   8, itplot, nplot); }
        if (rnst >= 32) { figs_poisson_graphs(psP, o, rnst, level,  16, itplot, nplot); }
        rfact = rfact * 2;
      }
  }

void figs_poisson_metbar(PSStream **psP, plot_options_t *o, int nst)
  {
    int xst = 4;  /* Sampling step. */
    
    double xsz = nst*xst;          /* Length of bar. */
    double xmrg = 12;              /* Margin at sides of bar. */
    double xmin = xmrg;            /* Start of bar. */
    double xmax = xmin + xsz;      /* End of bar. */
    
    double ysz = 2;                /* Height of bar. */
    double ymrg = 8;               /* Space at top and bottom. */
    double ymin = ymrg;
    double ymax = ymin + ysz;
    
    double xtot = xsz + 2*xmrg;  /* Width of figure. */
    double ytot = ysz + 2*ymrg;  /* height of figure. */

    interval_t bbox[2] = { (interval_t){{ 0, xtot }}, (interval_t){{ 0, ytot }} };
    start_figure(psP, o, "metbar", NULL, bbox, 1.0);
    
    /* Ice boxes: */
    double xbmin = -7.5;
    double xbmax = +0.5;
    double ybd = 4.00;
    pswr_set_fill_color(*psP, 0.300, 0.800, 1.000);
    pswr_set_pen(*psP, 0,0,0, 0.15, 0,0);
    pswr_rectangle(*psP, xmin+xbmin, xmin+xbmax, ymin-ybd, ymax+ybd, TRUE, TRUE);
    pswr_rectangle(*psP, xmax-xbmax, xmax-xbmin, ymin-ybd, ymax+ybd, TRUE, TRUE);
    
    /* Metal bar: */
    pswr_set_fill_color(*psP, 0.800, 0.500, 0.400);
    pswr_set_pen(*psP, 0,0,0, 0.15, 0,0);
    pswr_rectangle(*psP, xmin, xmax, ymin, ymax, TRUE, TRUE);
    
    /* Heaters/coolers: */
    figs_poisson_bar_heaters(psP, o, nst, xmin, xmax, ymin);
    
    finish_figure(psP, o);
  }

void figs_poisson_disbar(PSStream **psP, plot_options_t *o, int nst)
  {
    int xst = 4;  /* Sampling step. */
    
    double dotr = 1.50; /* Bar element radius. */
    
    double xsz = nst*xst;          /* Length of bar. */
    double xmrg = 12;              /* Margin at sides of bar. */
    double xmin = xmrg;            /* Start of bar. */
    double xmax = xmin + xsz;      /* End of bar. */
    
    double ysz = 2*dotr;           /* Height of bar. */
    double ymrg_bot = 6;           /* Space at bottom. */
    double ymrg_top = 10;          /* Space at top. */
    double ymin = ymrg_bot;
    double ymed = ymin + dotr;
    double ymax = ymin + ysz;
    double ylab = ymax + 1.00*dotr;
    
    double fontsz = 14.0;
    
    double xtot = xsz + 2*xmrg;  /* Width of figure. */
    double ytot = ysz + ymrg_bot + ymrg_top;  /* height of figure. */

    interval_t bbox[2] = { (interval_t){{ 0, xtot }}, (interval_t){{ 0, ytot }} };
    start_figure(psP, o, "disbar", NULL, bbox, 1.0);
    
    /* Connection between dots: */
    pswr_set_pen(*psP, 0,0,0, 0.50, 0,0);
    pswr_segment(*psP, xmin, ymed, xmax, ymed);
    
    /* Sampling points: */
    int i;
    pswr_set_pen(*psP, 0,0,0, 0.15, 0,0);
    
    for (i = 0; i <= nst; i++)
      { double xt = xmin + i*xst;
      
        /* Element label: */
        if (i == 0)
          { pswr_set_label_font(*psP, "TimesRoman", fontsz);
            pswr_label(*psP, "0", xt, ylab, 0.0, 0.5, 0.0);
          }
        else if (i == nst)
          { pswr_set_label_font(*psP, "TimesItalic", fontsz);
            pswr_label(*psP, "n", xt, ylab, 0.0, 0.5, 0.0);
          }

        /* Bar element: */
        if ((i == 0) || ( i == nst))
          { pswr_set_fill_color(*psP, 0.300, 0.800, 1.000); }
        else
          { pswr_set_fill_color(*psP, 0.800, 0.500, 0.400); }
        pswr_circle(*psP, xt, ymed, dotr, TRUE, TRUE);
      }
    
    /* Heaters/coolers: */
    figs_poisson_bar_heaters(psP, o, nst, xmin, xmax, ymin);
    
    finish_figure(psP, o);
  }    
    
void figs_poisson_bar_heaters(PSStream **psP, plot_options_t *o, int nst, double xmin, double xmax, double ymin)
  {
    double xst = (xmax - xmin)/nst;  /* Spacing between arrows. */
    srandom(344615);
    int i;
    for (i = 0; i < 1000; i++) { (void)drandom(); } /* Random generator warmup. */
    double xsr = 0.25;  /* Half-width of stem. */
    double xtr = 0.80;  /* Half-width of triangle. */
    double yth = 2.00;  /* Height of triangle. */
    double ybas = ymin - 2.00*yth;  /* Y of triangle base. */
    double yalo = ybas - yth;  /* Y of bottomof arrow. */
    double yahi = ybas + yth;  /* Y of bottomof arrow. */
    for (i = 1; i < nst; i++)
      { double xt = xmin + i*xst;
        if (drandom() < 0.400)
          { /* Choose between heat or cold source: */
            if (drandom() < 0.500) 
              { /* Heat source, up arrow. */
                pswr_set_fill_color(*psP, 1.000, 0.200, 0.000); 
                pswr_rectangle(*psP, xt-xsr, xt+xsr, yalo, ybas+0.5*yth, TRUE, FALSE);
                pswr_triangle(*psP, xt-xtr, ybas, xt, yahi, xt+xtr, ybas, TRUE, FALSE);
              }
            else
              { /* Cold source, down arrow. */
                pswr_set_fill_color(*psP, 0.300, 0.800, 1.000);
                pswr_rectangle(*psP, xt-xsr, xt+xsr, ybas-0.5*yth, yahi, TRUE, FALSE);
                pswr_triangle(*psP, xt-xtr, ybas, xt, yalo, xt+xtr, ybas, TRUE, FALSE);
              }
          }
      }
  }

void figs_poisson_graphs
  ( PSStream **psP, 
    plot_options_t *o, 
    int nst, 
    int level, 
    int freq, 
    int itplot[], 
    int nplot
  )
  {
    int nit = itplot[nplot-1] + 1; /* Number of iterations. */
    double T[nst+1]; /* Temperatures of elements computed by Gauss-Seidel. */
    double P[nst+1]; /* Heat input-output power at each element. */
    double S[nst+1]; /* True solution for temperatures. */
    double E[nst+1]; /* Error {T-S}. */
    
    double kappa = 3.2/nst;
    
    frgb_t red = (frgb_t){{ 1.000f, 0.000f, 0.000f }};
    frgb_t grn = (frgb_t){{ 0.200f, 0.800f, 0.000f }};
    frgb_t org = (frgb_t){{ 1.000f, 0.500f, 0.000f }};
    frgb_t blu = (frgb_t){{ 0.000f, 0.500f, 1.000f }};
    
    /* Get and plot the input-output power term: */
    figs_poisson_set_power(P, nst, level);
    if (freq == 0) 
      { double Pmin = -0.25*(double)ipow(2,level);
        double Pmax = +1.25*(double)ipow(2,level);
        figs_poisson_graphs_plot(psP, o, P, Pmin, Pmax, nst, org, "pwr", level, 0, 0);
      }

    /* Get and plot the true solution: */
    figs_poisson_get_true(P, S, nst, kappa);
    if (freq == 0) 
      { double Tmin = -0.25;
        double Tmax = +1.25;
        figs_poisson_graphs_plot(psP, o, S, Tmin, Tmax, nst, blu, "tru", level, 0, 0);
      }
    
    /* File of error per iteration: */
    char *errfname = jsprintf("%s-eit-v%02d-f%03d.txt", o->outName, level, freq);
    FILE *errwr = open_write(errfname, TRUE);
    free(errfname);

    /* Gauss-Seidel iteration: */
    int kplot = 0;
    int it;
    for (it = 0; it < nit; it++)
      { /* Get initial solution or update current one: */
        if (it == 0)
          { figs_poisson_graphs_get_guess(P, kappa, S, T, nst, freq); }
        else
          { figs_poisson_graphs_iterate(T, P, nst, kappa); }
        /* Compute the error: */
        figs_poisson_graphs_compute_error(T, S, E, nst);
        
        /* Save the root-means-square error: */
        double emag = sqrt(rn_norm_sqr(nst+1, E)/(nst-1));
        fprintf(errwr, "%8d %10.5e\n", it+1, emag);
        
        /* Plot current solution and error: */
        if (it == itplot[kplot])
          { figs_poisson_graphs_plot(psP, o, T, -1.25, +2.25, nst, grn, "cur", level, freq, it);
            figs_poisson_graphs_plot(psP, o, E, -1.25, +1.25, nst, red, "err", level, freq, it);
            kplot++;
          }
      }
    
    fclose(errwr);
  }
  
void figs_poisson_set_power(double P[], int nst, int level)  
  {
    if (level == 0)
      { demand(nst % 2 == 0, "num steps must be even");
        int i;
        for (i = 0; i <= nst; i++) { P[i] = 0.0; }
        P[nst/2+1] = 1.00;
      }
    else
      { int nstU = 2*nst; /* Size of {level-1} problem */
        double PU[nstU+1];
        figs_poisson_set_power(PU, nstU, level-1);
        figs_poisson_downscale_power(PU, nstU, P, nst);
      }
  }
 
void figs_poisson_get_true (double P[], double S[], int nst, double kappa)   
  {
    /* Compute a solution with zero flow on the left: */
    double F = 0; /* Heat flow from {i-1} to {i}. */
    S[0] = 0;
    int i;
    for (i = 1; i <= nst; i++) { S[i] = S[i-1] - kappa*F; F += P[i]; }
    /* Adjust boundary conditions to 0: */
    double B = S[0];
    double A = (S[nst] - S[0])/nst;
    for (i = 0; i <= nst; i++) { S[i] = S[i] - (A*i + B); }
  }

void figs_poisson_graphs_get_guess(double P[], double kappa, double S[], double T[], int nst, int freq)
  {
    int i;
    T[0] = 0.0;
    if (freq == 0)
      { 
        srandom(417417417);
        for (i = 0; i < 1000; i++) { (void)drandom(); }
        for (i = 1; i < nst; i++) { T[i] = 2*drandom() - 1; }
      }
    else if (freq == 999)
      { 
        int nstR = nst/2;
        double PR[nstR+1];
        figs_poisson_downscale_power(P, nst, PR, nstR);
        double SR[nstR+1];
        figs_poisson_get_true (PR, SR, nstR, 2*kappa);   
        figs_poisson_graphs_upscale_temp(SR, nstR, T, nst);
      }
    else
      { for (i = 1; i < nst; i++) 
          { double phase = (M_PI*freq*i)/nst;
            T[i] = S[i] + sin(phase);
          }
      }
    T[nst] = 0;
  }

void figs_poisson_graphs_iterate(double T[], double P[], int nst, double kappa)
  {
    double U[nst+1];
    U[0] = 0.0;
    int i;
    for (i = 1; i < nst; i++) { U[i] = 0.5*(T[i-1] + T[i+1]) + 0.5*kappa*P[i]; }
    U[nst] = 0;
    for (i = 0; i <= nst; i++) { T[i] = U[i]; }
  }

void figs_poisson_graphs_compute_error(double T[], double S[], double E[], int nst)
  {
    int i;
    for (i = 0; i <= nst; i++) { E[i] = T[i] - S[i]; }
  }

void figs_poisson_graphs_downscale_temp(double U[], int nst, double UR[], int nstR)
  {
    demand(nst == 2*nstR, "bad scaling");
    UR[0] = U[0];
    int i;
    for (i = 1; i < nstR; i++) { UR[i] = 0.5*U[2*i] + 0.25*(U[2*i-1] + U[2*i+1]); }
    UR[nstR] = U[nst];
  }

void figs_poisson_downscale_power(double P[], int nst, double PR[], int nstR)
  {
    demand(nst == 2*nstR, "bad scaling");
    demand(P[0] == 0, "bad power at 0");
    demand(P[nst] == 0, "bad power at {nst}");
    PR[0] = P[0];
    int i;
    for (i = 1; i < nstR; i++) { PR[i] = P[2*i] + 0.5*(P[2*i-1] + P[2*i+1]); }
    PR[nstR] = P[nst];
  }

void figs_poisson_graphs_upscale_temp(double UR[], int nstR, double U[], int nst)
  {
    demand(nst == 2*nstR, "bad scaling");
    U[0] = UR[0];
    int i;
    for (i = 1; i < nst; i++) 
      { U[i] = ((i % 2) == 0 ? UR[i/2] : 0.5*(UR[(i-1)/2] + UR[(i+1)/2])); }
    U[nst] = UR[nstR];
  }

void figs_poisson_graphs_plot
  ( PSStream **psP, 
    plot_options_t *o, 
    double Y[], 
    double Ymin,
    double Ymax,
    int nst, 
    frgb_t dotcolor,
    char *tag, 
    int level, 
    int freq, 
    int it
  )
  {
    int xst = 4;  /* Sampling step. */
    
    double xsz = nst*xst;          /* Length of bar. */
    double xmrg = 8;               /* Margin at sides of bar. */
    double xmin = xmrg;            /* Start of bar. */
    double xmax = xmin + xsz;      /* End of bar. */
    
    double Yscale = 10;
    
    double ysz = Yscale*(Ymax-Ymin); /* Height of graph. */
    double ymrg_bot = 4;             /* Space at bottom. */
    double ymrg_top = 4;             /* Space at top. */
    double ymin = ymrg_bot;
    /* double ymax = ymin + ysz; */
    double yzer = ymin - Yscale*Ymin; /* Y position of zero axis. */
    
    double dotr = 0.75; /* Dot radius (mm). */
    
    double xtot = xsz + 2*xmrg;  /* Width of figure. */
    double ytot = ysz + ymrg_bot + ymrg_top;  /* height of figure. */
    interval_t bbox[2] = { (interval_t){{ 0, xtot }}, (interval_t){{ 0, ytot }} };
 
    char *fulltag = jsprintf("%s-v%02d-f%03d-i%06d", tag, level, freq, it);
    start_figure(psP, o, fulltag, NULL, bbox, 1.0);
    free(fulltag);
    
    if ((Ymin < 0) && (Ymax > 0))
      { /* Zero axis: */
        pswr_set_pen(*psP, 0,0,0, 0.50, 0,0);
        pswr_segment(*psP, xmin, yzer, xmax, yzer);
      }
    
    /* Bar-and-dot graph: */
    pswr_set_fill_color(*psP, dotcolor.c[0], dotcolor.c[1], dotcolor.c[2]); 
    pswr_set_pen(*psP, 0,0,0, 0.15, 0,0);
    int i;
    for (i = 0; i <= nst; i++)
      { double xt = xmin + i*xst;
        double yt = yzer+Y[i]*Yscale;
        pswr_segment(*psP, xt, yzer, xt, yt);
        pswr_dot(*psP, xt, yt, dotr, TRUE, TRUE);
      }
    
    finish_figure(psP, o);
  }

#define ARG_ERROR(Msg,Arg) arg_error((Msg),(Arg),argv[0],usage)
  
#define GET_STRING(Var) get_arg_string(&(Var), &argn, argc, argv, usage)
#define GET_DOUBLE(Var) get_arg_double(&(Var), &argn, argc, argv, usage)
#define GET_FLOAT(Var) get_arg_float(&(Var), &argn, argc, argv, usage)

void parse_options(int argc, char **argv, plot_options_t **op)
  {
    plot_options_t *o = (plot_options_t *)malloc(sizeof(plot_options_t));
    char* usage = "\n  [ -help ] [ -outName PREFIX ] [ -eps | -ps ]";
    int argn;

    /* Defaults: */
    o->outName = "dg";
    o->eps = FALSE;
    o->paperSize = "letter";
    o->color = TRUE;
    o->scale = 25.4;
    o->margin[0] = 1.0; o->margin[1] = 1.0;
    o->meshSize = 1.0;

    /* Default depends on {color} option: */
    frgb_t defColorFill = (frgb_t){{ 1.00f, 0.95f, 0.85f }};
    frgb_t defGrayFill  = (frgb_t){{ 0.90f, 0.90f, 0.90f }};
    o->fill = (frgb_t){{-1.0, -1.0, -1.0}};

    argn = 1;

    /* Scan command line options. */
    while ((argn < argc) && (argv[argn][0] == '-') && (argv[argn][1] != '\0'))
      {
        char *key = argv[argn];
        if ((key[0] == '-') && (key[1] == '-') && (key[2] != '\0')) { key++; }
        if (strcmp(key, "-help") == 0)
          { fprintf(stderr, "usage: %s %s\n", argv[0], usage); exit(0); }
        else if (strcmp(key, "-outName") == 0)
          { GET_STRING(o->outName); }
        else if (strcmp(key, "-eps") == 0)
          { o->eps = TRUE; }
        else if (strcmp(key, "-ps") == 0)
          { o->eps = FALSE; }
        else if (strcmp(key, "-paperSize") == 0)
          { GET_STRING(o->paperSize); }
        else if (strcmp(key, "-color") == 0)
          { o->color = TRUE; }
        else if ((strcmp(key, "-gray") == 0) || (strcmp(key, "-grey") == 0))
          { o->color = FALSE; }
        else if (strcmp(key, "-fill") == 0)
          { GET_FLOAT(o->fill.c[0]);
            GET_FLOAT(o->fill.c[1]);
            GET_FLOAT(o->fill.c[2]);
          }
        else if (strcmp(key, "-scale") == 0)
          { GET_DOUBLE(o->scale); }
        else if (strcmp(key, "-margin") == 0)
          { GET_DOUBLE(o->margin[0]);
            GET_DOUBLE(o->margin[1]);
          }
        else if (strcmp(key, "-meshSize") == 0)
          { GET_DOUBLE(o->meshSize);
          }
        else 
          { ARG_ERROR("unknown option", argv[argn]); }
        ++argn;
      }
    if (argn != argc) { ARG_ERROR("extraneous arguments", argv[argn]); }

    if (o->fill.c[0] == -1) 
      { o->fill = (o->color ? defColorFill : defGrayFill); }
      
    (*op) = o;
  }

void get_arg_string(char **varp, int *argnp, int argc, char **argv, char *usage)
   /*
     Stores the next command line argument (as a string) into "*varp" */
  {
    int argn = *argnp;
    if (argn+1 >= argc)
      { ARG_ERROR("missing arg value", argv[argn]); }
    (*varp) = argv[argn+1];
    (*argnp) = argn+1;
  }
  
void get_arg_double(double *varp, int *argnp, int argc, char **argv, char *usage)
   /*
     Stores the next command line argument (as a double) into "*varp" */
  {
    int argn = *argnp;
    char *end;
    if (argn+1 >= argc)
      { ARG_ERROR("missing arg value", argv[argn]); }
    (*varp) = strtod(argv[argn+1], &end);
    if ((*end) != '\0') 
      { ARG_ERROR("invalid numeric argument", argv[argn+1]); }
    (*argnp) = argn+1;
  }
  
void get_arg_float(float *varp, int *argnp, int argc, char **argv, char *usage)
   /*
     Stores the next command line argument (as a double) into "*varp" */
  {
    int argn = *argnp;
    char *end;
    if (argn+1 >= argc)
      { ARG_ERROR("missing arg value", argv[argn]); }
    (*varp) = (float)strtod(argv[argn+1], &end);
    if ((*end) != '\0') 
      { ARG_ERROR("invalid numeric argument", argv[argn+1]); }
    (*argnp) = argn+1;
  }

void arg_error(char *msg, char *arg, char *pname, char *usage)
   /*
     Prints "msg", "arg", the "usage" message, and exits.
     Handy while parsing the command line arguments. */
  {
    fprintf(stderr, "%s %s\n", msg, arg);
    fprintf(stderr, "usage: %s %s\n", pname, usage);
    exit(1);
  }