#include <fstream>
#include "iohelp.h"
#include "plotting.h"

// *********************************************************************

// useful global variables

const double symbol_size=0.6;
int line_thickness=1,Rplot_type;
double xmin,xmax,ymin,ymax;
double x_major_tick_spacing,y_major_tick_spacing;
int x_minor_tick_num,y_minor_tick_num;
charString title,subtitle,xlabel,ylabel;
plot_labels label_strings(1);

int Nchannels;
Svector channel_match,channel_label;
ifstream in;
ofstream html_out;
xy_curves goldstone;
xy_curves zero_line(1);

// ********************************************************************************

void plot_label::assign(const charString& label_string, 
                        double xloc, double yloc)
{
 label=label_string;
 label_xloc=xloc;
 label_yloc=yloc;
}

// *********************************************************************

plot_labels::plot_labels()
{n=0;}

plot_labels::plot_labels(int dim)
{n=dim;
 if (dim<1){
    cerr << "negative or zero dimension in plot_labels"<<endl;
    exit(1);}
 m=new plot_label[n];
 if (!m){
    cerr << " problem allocating memory"<<endl;
    exit(1);}
 --m;
}

plot_labels::plot_labels(const plot_labels& pvec)
{int i;
 n=pvec.n;
 if (n>0){
    m=new plot_label[n];
    if (!m){
       cerr << " problem allocating memory"<<endl;
       exit(1);}
    m--;
    for (i=1;i<=n;i++) m[i]=pvec.m[i];
    }
}


int plot_labels::dim() const
{return n;}

plot_labels::~plot_labels()
{if (n>0) delete [] ++m;}

void plot_labels::assign(int index, const charString& label_string,
                         double xloc, double yloc)
{
 if ((index>=1)&&(index<=n)){
    m[index].assign(label_string,xloc,yloc);
    }
 else{
    cerr << " bad index; no assignment made in plot_labels"<<endl;
    }
}

void plot_labels::reset(int dim)
{
 if (n!=dim){
    if (n>0) delete [] ++m;
    n=dim;
    if (n>0){
       m=new plot_label[n];
       if (!m){
          cerr << " problem allocating memory"<<endl;
          exit(1);}
       m--;
    }}
}

plot_label &plot_labels::operator()(int i) const
 {return ((i>=1)&&(i<=n)) ? m[i]
    : (cerr << "plot_labels out of bounds"<<endl,exit(-1),m[1]);}



 // ********************************************************************

   //  This routine takes a '\0'-terminated character string in "ff" and returns
   //  all of the tokens.  A token is either a sequence on non-blank characters
   //  of a sequence of any characters between parentheses.
   
void get_tokens(char *ff, Svector& tokens)
{
 int i,count;
 char ch,ch2;
 char *fptr;

 count=i=0;
 ch=ff[i];
 while (ch==' ') ch=ff[++i];
 while (ch!='\0'){
    count++;
    if (ch=='(') ch2=')';
    else ch2=' ';
    while ((ch!=ch2)&&(ch!='\0')) ch=ff[++i];
    if (ch==')') ch=ff[++i];
    while (ch==' ') ch=ff[++i];
    }
 if (count==0){
    tokens(); return;}
 tokens.reset(count);
 i=count=0;
 ch=ff[i];
 while (ch==' ') ch=ff[++i];
 while (ch!='\0'){
    count++;
    fptr=&ff[i];
    if (ch=='(') ch2=')';
    else ch2=' ';
    while ((ch!=ch2)&&(ch!='\0')) ch=ff[++i];
    if (ch==')') ch=ff[++i];
    ch2=ff[i];
    ff[i]='\0';
    tokens(count)=charString(fptr);
    ff[i]=ch2;
    while (ch==' ') ch=ff[++i];
    }
}

// ****************************************************************

  //  linked list for results (used before storing in arrays)

struct node
{
 double R_value;
 double Rinv_value;
 double mean_value;
 double stdev_value;
 node *next;
};

node *head=NULL;
node *tail=NULL;

void add_node(double R, double Rinv, double ans, double err)
{
 if (head==NULL){
    head=new node;
    head->next=NULL;
    head->R_value=R;
    head->Rinv_value=Rinv;
    head->mean_value=ans;
    head->stdev_value=err;
    tail=head;
    }
 else{
    tail->next=new node;
    tail=tail->next;
    tail->R_value=R;
    tail->Rinv_value=Rinv;
    tail->mean_value=ans;
    tail->stdev_value=err;
    tail->next=NULL;
    }
}

void delete_nodes()
{
 if (head==NULL) return;
 node *pnode=head;
 node *pnext=head->next;
 while (pnext != NULL){
    delete pnode;
    pnode=pnext;
    pnext=pnext->next;
    }
 head=NULL;
 tail=NULL;
}

// ****************************************************************

 // Rplot_type:  1 => x-axis is R,  0 => x-axis is 1/R

int xydy_curve::setup_from_stream(ifstream& in,
                                  const charString& channel_string,
          int Rplot_type, int type_of_line, int color_of_line,
          int thickness_of_line, int type_of_symbol, 
          int color_of_symbol, int fill_of_symbol, 
          double size_of_symbol, const charString& label_string,
          double xloc_of_label, int yloc_of_label  )
{
 const int maxlinelength=132;
 static char ff[maxlinelength];
 Svector tokens;
 int flag;
 double Rvalue,Rinvvalue,ans,err;

 in.clear();  // in case eof reached
 in.seekg(0,ios::beg); // go back to beginning of file
 flag=Npoints=0;

   // first, search sequentially through file to find start
   // of relevant channel

 while (!flag){
    in.getline(ff,maxlinelength);
    get_tokens(ff,tokens);
    if (tokens.dim()==5){
       if (tokens(4).locate(channel_string)>=0) flag=1;
       }
    if (in.eof()){
       //cerr << " could not find requested channel"<<endl;
       wipeout();
       return 1;
       }
    }

   // next, loop through any blank lines to get to the data

 in.getline(ff,maxlinelength);
 get_tokens(ff,tokens);
 while (tokens.dim()==0){
    in.getline(ff,maxlinelength);
    get_tokens(ff,tokens);
    }
 
   // now begin extracting the data (until blank line)

 while (tokens.dim()==5){
    if (string_to_double(tokens(2),Rvalue)) break;
    if (string_to_double(tokens(3),Rinvvalue)) break;
    if (string_to_double(tokens(4),ans)) break;
    if (string_to_double(tokens(5),err)) break;
    Npoints++;
    add_node(Rvalue,Rinvvalue,ans,err);
    in.getline(ff,maxlinelength);
    get_tokens(ff,tokens);
    }
 xvals.reset(Npoints);
 yvals.reset(Npoints);
 dyvals.reset(Npoints);
 flag=0;
 node *pnode=head;
 while (pnode!=NULL){
   yvals(++flag)=pnode->mean_value;
   dyvals(flag)=pnode->stdev_value;
   if (Rplot_type) xvals(flag)=pnode->R_value;
   else xvals(flag)=pnode->Rinv_value;
   pnode=pnode->next;
   }
 delete_nodes();

 line_color=color_of_line;
 line_type=type_of_line;
 line_thickness=thickness_of_line;
 symbol_type=type_of_symbol;
 symbol_color=color_of_symbol;
 symbol_fill=fill_of_symbol;
 symbol_size=size_of_symbol;
 label=label_string;
 horizontal_label_loc=xloc_of_label;
 vertical_label_loc=yloc_of_label;
 return 0;
}


void xydy_curve::setup(int number_of_points, double *xvalues, 
                       double *yvalues, double *dyvalues, 
          int type_of_line, int color_of_line,
          int thickness_of_line, int type_of_symbol, 
          int color_of_symbol, int fill_of_symbol, 
          double size_of_symbol, const charString& label_string,
          double xloc_of_label, int yloc_of_label )
{
 int i;
 Npoints=number_of_points;
 xvals.reset(Npoints);
 yvals.reset(Npoints);
 dyvals.reset(Npoints);
 for (i=1;i<=Npoints;i++){
    xvals(i)=xvalues[i-1];
    yvals(i)=yvalues[i-1];
    dyvals(i)=dyvalues[i-1];
    }
 line_color=color_of_line;
 line_type=type_of_line;
 line_thickness=thickness_of_line;
 symbol_type=type_of_symbol;
 symbol_color=color_of_symbol;
 symbol_fill=fill_of_symbol;
 symbol_size=size_of_symbol;
 label=label_string;
 horizontal_label_loc=xloc_of_label;
 vertical_label_loc=yloc_of_label;
}


void xydy_curve::wipeout()
{
 Npoints=0;
 xvals();
 yvals();
 dyvals();
 label="";
 horizontal_label_loc=0.5;
 vertical_label_loc=0;
}


 // ******************************************************************

xydy_curves::xydy_curves()
{n=0;}

xydy_curves::xydy_curves(int dim)
{n=dim;
 if (dim<1){
    cerr << "negative or zero dimension in xydy_curves"<<endl;
    exit(1);}
 m=new xydy_curve[n];
 if (!m){
    cerr << " problem allocating memory"<<endl;
    exit(1);}
 --m;
}

xydy_curves::xydy_curves(const xydy_curves& pvec)
{int i;
 n=pvec.n;
 if (n>0){
    m=new xydy_curve[n];
    if (!m){
       cerr << " problem allocating memory"<<endl;
       exit(1);}
    m--;
    for (i=1;i<=n;i++) m[i]=pvec.m[i];
    }
}


int xydy_curves::dim() const
{return n;}

xydy_curves::~xydy_curves()
{if (n>0) delete [] ++m;}


int xydy_curves::setup_from_stream(int index, 
                                   ifstream& in, 
            const charString& channel_string,
            int Rplot_type, int type_of_line, int color_of_line,
            int thickness_of_line, int type_of_symbol, 
            int color_of_symbol, int fill_of_symbol, 
            double size_of_symbol, const charString& label_string,
            double xloc_of_label, int yloc_of_label  )
{
 if ((index>=1)&&(index<=n)){
    return m[index].setup_from_stream(in,channel_string,Rplot_type,
                               type_of_line,color_of_line,
                               thickness_of_line,type_of_symbol, 
                               color_of_symbol,fill_of_symbol, 
                               size_of_symbol,label_string,
                               xloc_of_label,yloc_of_label);
    }
 else{
    cerr << " bad index; no assignment made in xydy_curves"<<endl;
    return 1;
    }
}

void xydy_curves::setup(int index, 
                        int number_of_points, 
             double *xvalues, double *yvalues,
             double *dyvalues, int type_of_line, int color_of_line,
             int thickness_of_line, int type_of_symbol, 
             int color_of_symbol, int fill_of_symbol,
             double size_of_symbol, const charString& label_string,
             double xloc_of_label, int yloc_of_label)
{
 if ((index>=1)&&(index<=n)){
    m[index].setup(number_of_points, xvalues,yvalues,
                   dyvalues,type_of_line,color_of_line,
                   thickness_of_line,type_of_symbol, 
                   color_of_symbol,fill_of_symbol,
                   size_of_symbol,label_string,
                   xloc_of_label,yloc_of_label);
    }
 else{
    cerr << " bad index; no assignment made in xydy_curves"<<endl;
    }
}


void xydy_curves::reset(int dim)
{
 if (n!=dim){
    if (n>0) delete [] ++m;
    n=dim;
    if (n>0){
       m=new xydy_curve[n];
       if (!m){
          cerr << " problem allocating memory"<<endl;
          exit(1);}
       m--;
    }}
}

xydy_curve &xydy_curves::operator()(int i) const
 {return ((i>=1)&&(i<=n)) ? m[i]
    : (cerr << "xydy_curves out of bounds"<<endl,exit(-1),m[1]);}


 // ******************************************************************


void xy_curve::setup(int number_of_points, double *xvalues, 
                     double *yvalues, 
          int type_of_line, int color_of_line,
          int thickness_of_line, int type_of_symbol, 
          int color_of_symbol, int fill_of_symbol, 
          double size_of_symbol, const charString& label_string,
          double xloc_of_label, int yloc_of_label )
{
 int i;
 Npoints=number_of_points;
 xvals.reset(Npoints);
 yvals.reset(Npoints);
 for (i=1;i<=Npoints;i++){
    xvals(i)=xvalues[i-1];
    yvals(i)=yvalues[i-1];
    }
 line_color=color_of_line;
 line_type=type_of_line;
 line_thickness=thickness_of_line;
 symbol_type=type_of_symbol;
 symbol_color=color_of_symbol;
 symbol_fill=fill_of_symbol;
 symbol_size=size_of_symbol;
 label=label_string;
 horizontal_label_loc=xloc_of_label;
 vertical_label_loc=yloc_of_label;
}


void xy_curve::wipeout()
{
 Npoints=0;
 xvals();
 yvals();
 label="";
 horizontal_label_loc=0.5;
 vertical_label_loc=0;
}


 // ******************************************************************

xy_curves::xy_curves()
{n=0;}

xy_curves::xy_curves(int dim)
{n=dim;
 if (dim<1){
    cerr << "negative or zero dimension in xy_curves"<<endl;
    exit(1);}
 m=new xy_curve[n];
 if (!m){
    cerr << " problem allocating memory"<<endl;
    exit(1);}
 --m;
}

xy_curves::xy_curves(const xy_curves& pvec)
{int i;
 n=pvec.n;
 if (n>0){
    m=new xy_curve[n];
    if (!m){
       cerr << " problem allocating memory"<<endl;
       exit(1);}
    m--;
    for (i=1;i<=n;i++) m[i]=pvec.m[i];
    }
}


int xy_curves::dim() const
{return n;}

xy_curves::~xy_curves()
{if (n>0) delete [] ++m;}


void xy_curves::setup(int index, 
                      int number_of_points, 
             double *xvalues, double *yvalues,
             int type_of_line, int color_of_line,
             int thickness_of_line, int type_of_symbol, 
             int color_of_symbol, int fill_of_symbol,
             double size_of_symbol, const charString& label_string,
             double xloc_of_label, int yloc_of_label)
{
 if ((index>=1)&&(index<=n)){
    m[index].setup(number_of_points,xvalues,yvalues,
                   type_of_line,color_of_line,
                   thickness_of_line,type_of_symbol, 
                   color_of_symbol,fill_of_symbol,
                   size_of_symbol,label_string,
                   xloc_of_label,yloc_of_label);
    }
 else{
    cerr << " bad index; no assignment made in xy_curves"<<endl;
    }
}


void xy_curves::reset(int dim)
{
 if (n!=dim){
    if (n>0) delete [] ++m;
    n=dim;
    if (n>0){
       m=new xy_curve[n];
       if (!m){
          cerr << " problem allocating memory"<<endl;
          exit(1);}
       m--;
    }}
}

xy_curve &xy_curves::operator()(int i) const
 {return ((i>=1)&&(i<=n)) ? m[i]
    : (cerr << "xy_curves out of bounds"<<endl,exit(-1),m[1]);}

 // *******************************************************************
 //
 //    Computes the location of the curve label.
 //    The location of the label is specified as follows:
 //
 //      HORIZONTAL SPECIFICATION
 //
 //          xloc_of_label < 0   =>  left of smallest-x point
 //          xloc_of_label > 1   =>  right of largest-x point
 //      0 < xloc_of_label < 1   =>  fractional distance between
 //                                  smallest and largest x-values
 //
 //      VERTICAL SPECIFICATION
 //
 //          yloc_of_label > 0   =>  above the curve
 //          yloc_of_label = 0   =>  on the curve
 //          yloc_of_label < 0   =>  below the curve
 //

 
void compute_label_location(const dvector& xvals, const dvector& yvals,
                            const dvector& dyvals, 
                            double horizontal_label_loc, int vertical_label_loc,
                            double xmin, double xmax, double ymin, double ymax,
                            double& x_pos, double& y_pos)
{
 int Npoints,Nerror,i;
 double xstart,xstop,dy,dx;
 
 Npoints=xvals.dim(); 
 Nerror=dyvals.dim();
 
 if ((horizontal_label_loc<=0.0)||(Npoints<2)){   // put label to left of smallest-x point
    x_pos=xmax; y_pos=ymax;
    for (i=1;i<=Npoints;i++){
       if ((xvals(i)<x_pos)&&(xvals(i)>=xmin)){
          x_pos=xvals(i);
          y_pos=yvals(i);
          }
       }
    x_pos-=0.05*(xmax-xmin);
    if (vertical_label_loc>0) y_pos+=0.06*(ymax-ymin);
    else if (vertical_label_loc<0) y_pos-=0.06*(ymax-ymin);
    }
 else if (horizontal_label_loc>=1.0){   // put label to right of largest-x point
    x_pos=xmin; y_pos=ymin;
    for (i=1;i<=Npoints;i++){
       if ((xvals(i)>x_pos)&&(xvals(i)<=xmax)){
          x_pos=xvals(i);
          y_pos=yvals(i);
          }
       }
    x_pos+=0.05*(xmax-xmin);
    if (vertical_label_loc>0) y_pos+=0.06*(ymax-ymin);
    else if (vertical_label_loc<0) y_pos-=0.06*(ymax-ymin);
    }
 else{                  // put label fractional distance horizontal between smallest
    xstart=xmax; xstop=xmin;                   // and largest-x values
    for (i=1;i<=Npoints;i++){
       if ((xvals(i)<xstart)&&(xvals(i)>=xmin)) xstart=xvals(i);
       if ((xvals(i)>xstop)&&(xvals(i)<=xmax)) xstop=xvals(i);
       }
    x_pos=xstart+(xstop-xstart)*horizontal_label_loc;
    int k1=0,k2=0;
    double x1=xmin,x2=xmax,y1=ymin,y2=ymax;
    for (i=1;i<=Npoints;i++){
       if ((xvals(i)<x_pos)&&(xvals(i)>x1)){
          x1=xvals(i); y1=yvals(i); k1=i;}
       if ((xvals(i)>x_pos)&&(xvals(i)<x2)){
          x2=xvals(i); y2=yvals(i); k2=i;}
       }
    if ((k1==0)||(k2==0)||(k1==k2)) return; 
    y_pos=(x1*y2-y1*x2+x_pos*(y1-y2))/(x1-x2);   // linear interpolation
    if (vertical_label_loc==0) return;
    dy=(Nerror==0)?0.0:((dyvals(k1)>dyvals(k2))?dyvals(k1):dyvals(k2));
    double r1=1.0/(xmax-xmin), r2=1.0/(ymax-ymin);
    double dh,dhmax=0.0,h=0.07+dy*r2;
    if (vertical_label_loc<0) h=-h;
    y_pos+=(ymax-ymin)*h;
    for (i=1;i<=Npoints;i++){
       dx=r1*(xvals(i)-x_pos); dy=fabs(r2*(yvals(i)-y_pos));
       if ((dx*dx+dy*dy)<h*h){
          dh=sqrt(h*h-dx*dx)-dy;
          if (dh>dhmax) dhmax=dh;
          }
       }
    if (vertical_label_loc<0) dhmax=-dhmax;
    y_pos+=dhmax*(ymax-ymin); 
    }
}


void compute_label_location(const dvector& xvals, const dvector& yvals,
                            double horizontal_label_loc, int vertical_label_loc,
                            double xmin, double xmax, double ymin, double ymax,
                            double& x_pos, double& y_pos)
{
 dvector dyvals;
 compute_label_location(xvals,yvals,dyvals,horizontal_label_loc,vertical_label_loc,
                        xmin,xmax,ymin,ymax,x_pos,y_pos);
}

 //  *****************************************************************
 //
 //   This procedure produces a plot
 //
 // *******************************************************************

void produce_plot(const charString& filename,
                  const charString& title, const charString& subtitle,
                  const charString& xlabel, const charString& ylabel,
                  double xmin, double xmax, double ymin, double ymax,
                  double x_major_tick_spacing, int x_minor_tick_num,
                  double y_major_tick_spacing, int y_minor_tick_num,
                  const xydy_curves& xydy_results, 
                  const xy_curves& xy_curves_to_plot,
                  const plot_labels& label_strings)
{
 int i,j;
 charString str;
 double x_pos,y_pos;
 
   // do some checks

 if ((xmax<=xmin)||(ymax<=ymin)){
    cerr << "ERROR in produce_plot: problem with plotting limits"<<endl;
    return;}

 ofstream out(convert_to_charstar(filename));
 if (!out){
    cerr << " could not open file for output"<<endl;
    return;}

 out << "@version 50106"<<endl;

    // Send some initialization commands to Grace

 out << "@ view xmin 0.200000"<<endl;
 out << "@ view xmax 0.950000"<<endl;
 out << "@ view ymin 0.150000"<<endl;
 out << "@ view ymax 0.850000"<<endl;

 out << "@ world xmin "<<xmin<<endl;
 out << "@ world xmax "<<xmax<<endl;
 out << "@ world ymin "<<ymin<<endl;
 out << "@ world ymax "<<ymax<<endl;
 out << "@ xaxes scale Normal"<<endl;
 out << "@ yaxes scale Normal"<<endl;
 out << "@ xaxes invert off"<<endl;
 out << "@ yaxes invert off"<<endl;

 if (title.non_blank()){
    out << "@ title \""<<title<<"\""<<endl;
    out << "@ title font 0"<<endl;
    out << "@ title size 1.250000"<<endl;}
 if (subtitle.non_blank()){
    out << "@ subtitle \""<<subtitle<<"\""<<endl;
    out << "@ subtitle font 0"<<endl;
    out << "@ subtitle size 1.040000"<<endl;}

 out << "@ xaxis  on"<<endl;
 out << "@ xaxis  type zero false"<<endl;
 out << "@ xaxis  offset 0.000000 , 0.000000"<<endl;
 out << "@ xaxis  bar on"<<endl;
 out << "@ xaxis  bar color 1"<<endl;
 out << "@ xaxis  bar linewidth 2.0"<<endl;
 if (xlabel.non_blank()){
    out << "@ xaxis  label \""<<xlabel<<"\""<<endl;
    out << "@ xaxis  label layout para"<<endl;
    out << "@ xaxis  label place auto"<<endl;
    out << "@ xaxis  label char size 1.830000"<<endl;
    out << "@ xaxis  label font 1"<<endl;
    out << "@ xaxis  label color 1"<<endl;
    out << "@ xaxis  label place normal"<<endl;}
 out << "@ xaxis  tick on"<<endl;
 out << "@ xaxis  tick major "<<x_major_tick_spacing<<endl;
 out << "@ xaxis  tick minor ticks "<<x_minor_tick_num<<endl;
 out << "@ xaxis  tick default 6"<<endl;
 out << "@ xaxis  tick place rounded true"<<endl;
 out << "@ xaxis  tick in"<<endl;
 out << "@ xaxis  tick major size 1.440000"<<endl;
 out << "@ xaxis  tick major color 1"<<endl;
 out << "@ xaxis  tick major linewidth 2.0"<<endl;
 out << "@ xaxis  tick major linestyle 1"<<endl;
 out << "@ xaxis  tick major grid off"<<endl;
 out << "@ xaxis  tick minor color 1"<<endl;
 out << "@ xaxis  tick minor linewidth 2.0"<<endl;
 out << "@ xaxis  tick minor linestyle 1"<<endl;
 out << "@ xaxis  tick minor grid off"<<endl;
 out << "@ xaxis  tick minor size 0.720000"<<endl;
 out << "@ xaxis  ticklabel on"<<endl;
 out << "@ xaxis  ticklabel format general"<<endl;
 out << "@ xaxis  ticklabel prec 5"<<endl;
 out << "@ xaxis  ticklabel formula \"\""<<endl;
 out << "@ xaxis  ticklabel append \"\""<<endl;
 out << "@ xaxis  ticklabel prepend \"\""<<endl;
 out << "@ xaxis  ticklabel angle 0"<<endl;
 out << "@ xaxis  ticklabel skip 0"<<endl;
 out << "@ xaxis  ticklabel stagger 0"<<endl;
 out << "@ xaxis  ticklabel place normal"<<endl;
 out << "@ xaxis  ticklabel offset auto"<<endl;
 out << "@ xaxis  ticklabel offset 0.000000 , 0.010000"<<endl;
 out << "@ xaxis  ticklabel start type auto"<<endl;
 out << "@ xaxis  ticklabel start 0.000000"<<endl;
 out << "@ xaxis  ticklabel stop type auto"<<endl;
 out << "@ xaxis  ticklabel stop 0.000000"<<endl;
 out << "@ xaxis  ticklabel char size 1.480000"<<endl;
 out << "@ xaxis  ticklabel font 0"<<endl;
 out << "@ xaxis  ticklabel color 1"<<endl;
 out << "@ xaxis  tick place both"<<endl;
 out << "@ xaxis  tick spec type none"<<endl;


 out << "@ yaxis  on"<<endl;
 out << "@ yaxis  type zero false"<<endl;
 out << "@ yaxis  offset 0.000000 , 0.000000"<<endl;
 out << "@ yaxis  bar on"<<endl;
 out << "@ yaxis  bar color 1"<<endl;
 out << "@ yaxis  bar linestyle 1"<<endl;
 out << "@ yaxis  bar linewidth 2.0"<<endl;
 if (ylabel.non_blank()){  
    out << "@ yaxis  label \""<<ylabel<<"\""<<endl;
    out << "@ yaxis  label layout para"<<endl;
    out << "@ yaxis  label place auto"<<endl;
    out << "@ yaxis  label char size 1.830000"<<endl;
    out << "@ yaxis  label font 1"<<endl;
    out << "@ yaxis  label color 1"<<endl;
    out << "@ yaxis  label place normal"<<endl;}
 out << "@ yaxis  tick on"<<endl;
 out << "@ yaxis  tick major "<<y_major_tick_spacing<<endl;
 out << "@ yaxis  tick minor ticks "<<y_minor_tick_num<<endl;
 out << "@ yaxis  tick default 6"<<endl;
 out << "@ yaxis  tick place rounded true"<<endl;
 out << "@ yaxis  tick in"<<endl;
 out << "@ yaxis  tick major size 1.440000"<<endl;
 out << "@ yaxis  tick major color 1"<<endl;
 out << "@ yaxis  tick major linewidth 2.0"<<endl;
 out << "@ yaxis  tick major linestyle 1"<<endl;
 out << "@ yaxis  tick major grid off"<<endl;
 out << "@ yaxis  tick minor color 1"<<endl;
 out << "@ yaxis  tick minor linewidth 2.0"<<endl;
 out << "@ yaxis  tick minor linestyle 1"<<endl;
 out << "@ yaxis  tick minor grid off"<<endl;
 out << "@ yaxis  tick minor size 0.720000"<<endl;
 out << "@ yaxis  ticklabel on"<<endl;
 out << "@ yaxis  ticklabel format general"<<endl;
 out << "@ yaxis  ticklabel prec 5"<<endl;
 out << "@ yaxis  ticklabel formula \"\""<<endl;
 out << "@ yaxis  ticklabel append \"\""<<endl;
 out << "@ yaxis  ticklabel prepend \"\""<<endl;
 out << "@ yaxis  ticklabel angle 0"<<endl;
 out << "@ yaxis  ticklabel skip 0"<<endl;
 out << "@ yaxis  ticklabel stagger 0"<<endl;
 out << "@ yaxis  ticklabel place normal"<<endl;
 out << "@ yaxis  ticklabel offset auto"<<endl;
 out << "@ yaxis  ticklabel offset 0.000000 , 0.010000"<<endl;
 out << "@ yaxis  ticklabel start type auto"<<endl;
 out << "@ yaxis  ticklabel start 0.000000"<<endl;
 out << "@ yaxis  ticklabel stop type auto"<<endl;
 out << "@ yaxis  ticklabel stop 0.000000"<<endl;
 out << "@ yaxis  ticklabel char size 1.480000"<<endl;
 out << "@ yaxis  ticklabel font 0"<<endl;
 out << "@ yaxis  ticklabel color 1"<<endl;
 out << "@ yaxis  tick place both"<<endl;
 out << "@ yaxis  tick spec type none"<<endl;


 for (i=1;i<=xydy_results.dim();i++){

    str="@ s"+int_to_string(i)+" ";
    out << str << "on"<<endl;
    out << str << "type xydy"<<endl;
    out << str << "symbol "<<xydy_results(i).symbol_type<<endl;
    out << str << "symbol size "<<xydy_results(i).symbol_size<<endl;
    out << str << "symbol color "<<xydy_results(i).symbol_color<<endl;
    out << str << "symbol pattern 1"<<endl;
    out << str << "symbol fill color "<<xydy_results(i).symbol_color<<endl;
    out << str << "symbol fill pattern "<<xydy_results(i).symbol_fill<<endl;
    out << str << "symbol linewidth 1.0"<<endl;
    out << str << "symbol linestyle 1"<<endl;
    out << str << "symbol char 65"<<endl;
    out << str << "symbol char font 0"<<endl;
    out << str << "symbol skip 0"<<endl;
    out << str << "line type 1"<<endl;
    out << str << "line linestyle "<<xydy_results(i).line_type<<endl;
    out << str << "line linewidth "<<xydy_results(i).line_thickness<<endl;
    out << str << "line color "<<xydy_results(i).line_color<<endl;
    out << str << "line pattern 1"<<endl;
    out << str << "baseline type 0"<<endl;
    out << str << "baseline off"<<endl;
    out << str << "dropline off"<<endl;
    out << str << "fill type 0"<<endl;
    out << str << "fill rule 0"<<endl;
    out << str << "fill color 1"<<endl;
    out << str << "fill pattern 1"<<endl;
    out << str << "avalue off"<<endl;
    out << str << "avalue type 2"<<endl;
    out << str << "avalue char size 1.0"<<endl;
    out << str << "avalue font 0"<<endl;
    out << str << "avalue color 1"<<endl;
    out << str << "avalue rot 0"<<endl;
    out << str << "avalue format general"<<endl;
    out << str << "avalue prec 3"<<endl;
    out << str << "avalue prepend \"\""<<endl;
    out << str << "avalue append \"\""<<endl;
    out << str << "avalue offset 0.0, 0.0"<<endl;
    out << str << "errorbar on"<<endl;
    out << str << "errorbar place both"<<endl;
    out << str << "errorbar color "<<xydy_results(i).symbol_color<<endl;
    out << str << "errorbar pattern 1"<<endl;
    out << str << "errorbar size 0.62"<<endl;
    out << str << "errorbar linewidth 2.0"<<endl;
    out << str << "errorbar linestyle 1"<<endl;
    out << str << "errorbar riser linewidth 2.0"<<endl;
    out << str << "errorbar riser linestyle 1"<<endl;
    out << str << "errorbar riser clip off"<<endl;
    out << str << "errorbar riser clip length 0.10"<<endl;

    out << "@target g0.s"<<i<<endl;
    out << "@type xydy"<<endl;
    for (j=1;j<=xydy_results(i).Npoints;j++){
       out <<" "<<xydy_results(i).xvals(j)
           <<" "<<xydy_results(i).yvals(j)
           <<" "<<xydy_results(i).dyvals(j)<<endl;
       }
    out << "&"<<endl;

    if (xydy_results(i).label.non_blank()){
       compute_label_location(xydy_results(i).xvals,xydy_results(i).yvals,
                              xydy_results(i).dyvals,
                              xydy_results(i).horizontal_label_loc,
                              xydy_results(i).vertical_label_loc,
                              xmin,xmax,ymin,ymax,x_pos,y_pos);
       out << "@with string"<<endl;
       out << "@ string on"<<endl;
       out << "@ string loctype world"<<endl;
       out << "@ string "<<x_pos<< ", "<<y_pos<<endl;
       out << "@ string color "<<xydy_results(i).symbol_color<<endl;
       out << "@ string rot 0"<<endl;
       out << "@ string font 3"<<endl;
       out << "@ string just 14"<<endl;
       out << "@ string char size 1.51"<<endl;
       out << "@ string def \""<<xydy_results(i).label<<"\""<<endl;
       }
    }

 for (i=1;i<=xy_curves_to_plot.dim();i++){

    str="@ s"+int_to_string(i+xydy_results.dim())+" ";
    out << str << "on"<<endl;
    out << str << "type xydy"<<endl;
    out << str << "symbol "<<xy_curves_to_plot(i).symbol_type<<endl;
    out << str << "symbol size "<<xy_curves_to_plot(i).symbol_size<<endl;
    out << str << "symbol color "<<xy_curves_to_plot(i).symbol_color<<endl;
    out << str << "symbol pattern 1"<<endl;
    out << str << "symbol fill color "<<xy_curves_to_plot(i).symbol_color<<endl;
    out << str << "symbol fill pattern "<<xy_curves_to_plot(i).symbol_fill<<endl;
    out << str << "symbol linewidth 1.0"<<endl;
    out << str << "symbol linestyle 1"<<endl;
    out << str << "symbol char 65"<<endl;
    out << str << "symbol char font 0"<<endl;
    out << str << "symbol skip 0"<<endl;
    out << str << "line type 1"<<endl;
    out << str << "line linestyle "<<xy_curves_to_plot(i).line_type<<endl;
    out << str << "line linewidth "<<xy_curves_to_plot(i).line_thickness<<endl;
    out << str << "line color "<<xy_curves_to_plot(i).line_color<<endl;
    out << str << "line pattern 1"<<endl;
    out << str << "baseline type 0"<<endl;
    out << str << "baseline off"<<endl;
    out << str << "dropline off"<<endl;
    out << str << "fill type 0"<<endl;
    out << str << "fill rule 0"<<endl;
    out << str << "fill color 1"<<endl;
    out << str << "fill pattern 1"<<endl;
    out << str << "avalue off"<<endl;
    out << str << "avalue type 2"<<endl;
    out << str << "avalue char size 1.0"<<endl;
    out << str << "avalue font 0"<<endl;
    out << str << "avalue color 1"<<endl;
    out << str << "avalue rot 0"<<endl;
    out << str << "avalue format general"<<endl;
    out << str << "avalue prec 3"<<endl;
    out << str << "avalue prepend \"\""<<endl;
    out << str << "avalue append \"\""<<endl;
    out << str << "avalue offset 0.0, 0.0"<<endl;

    out << "@target g0.s"<<i+xydy_results.dim()<<endl;
    out << "@type xy"<<endl;
    for (j=1;j<=xy_curves_to_plot(i).Npoints;j++){
       out <<" "<<xy_curves_to_plot(i).xvals(j)
           <<" "<<xy_curves_to_plot(i).yvals(j)<<endl;
       }
    out << "&"<<endl;

    if (xy_curves_to_plot(i).label.non_blank()){
       compute_label_location(xy_curves_to_plot(i).xvals,xy_curves_to_plot(i).yvals,
                              xy_curves_to_plot(i).horizontal_label_loc,
                              xy_curves_to_plot(i).vertical_label_loc,
                              xmin,xmax,ymin,ymax,x_pos,y_pos);
       out << "@with string"<<endl;
       out << "@ string on"<<endl;
       out << "@ string loctype world"<<endl;
       out << "@ string "<<x_pos<< ", "<<y_pos<<endl;
       out << "@ string color "<<xy_curves_to_plot(i).symbol_color<<endl;
       out << "@ string rot 0"<<endl;
       out << "@ string font 3"<<endl;
       out << "@ string just 14"<<endl;
       out << "@ string char size 1.51"<<endl;
       out << "@ string def \""<<xy_curves_to_plot(i).label<<"\""<<endl;
       }
    }

 for (i=1;i<=label_strings.dim();i++){
    if (label_strings(i).label.non_blank()){
       out << "@with string"<<endl;
       out << "@ string on"<<endl;
       out << "@ string loctype world"<<endl;
       out << "@ string "<<xmin+(xmax-xmin)*label_strings(i).label_xloc
           << ", "<<ymin+(ymax-ymin)*label_strings(i).label_yloc<<endl;
       out << "@ string color 1"<<endl;
       out << "@ string rot 0"<<endl;
       out << "@ string font 0"<<endl;
       out << "@ string just 0"<<endl;
       out << "@ string char size 1.0"<<endl;
       out << "@ string def \""<<label_strings(i).label<<"\""<<endl;
       }
    }

 out.close();
}

void produce_plot(const charString& filename,
                  const charString& title, const charString& subtitle,
                  const charString& xlabel, const charString& ylabel,
                  double xmin, double xmax, double ymin, double ymax,
                  double x_major_tick_spacing, int x_minor_tick_num,
                  double y_major_tick_spacing, int y_minor_tick_num,
                  const xydy_curves& xydy_results,
                  const plot_labels& label_strings)
{
 xy_curves xy_curves_to_plot;
 produce_plot(filename,title,subtitle,xlabel,ylabel,
              xmin,xmax,ymin,ymax,x_major_tick_spacing,x_minor_tick_num,
              y_major_tick_spacing,y_minor_tick_num,
              xydy_results,xy_curves_to_plot,label_strings);
}

void produce_plot(const charString& filename,
                  const charString& title, const charString& subtitle,
                  const charString& xlabel, const charString& ylabel,
                  double xmin, double xmax, double ymin, double ymax,
                  double x_major_tick_spacing, int x_minor_tick_num,
                  double y_major_tick_spacing, int y_minor_tick_num,
                  const xy_curves& xy_curves_to_plot,
                  const plot_labels& label_strings)
{
 xydy_curves xydy_results;
 produce_plot(filename,title,subtitle,xlabel,ylabel,
              xmin,xmax,ymin,ymax,x_major_tick_spacing,x_minor_tick_num,
              y_major_tick_spacing,y_minor_tick_num,
              xydy_results,xy_curves_to_plot,label_strings);
}


// ***********************************************************************


void setup_goldstone(int Rplot_type, double xmin, double xmax, double xi, double Zxi)
{
 const double pi=3.14159265358979324;
 const double eps=1e-10;
 const int Ngold_points=300;
 double xgold[Ngold_points];
 double ygold[Ngold_points];
 double step=(xmax-xmin)/Ngold_points;
 const int Ngoldstone=4;
 double rr;
 int i,j;
 
 goldstone.reset(Ngoldstone);
 rr=xmin;
 if ((fabs(rr)<eps)&&(Rplot_type)) rr+=step;
 for (j=0;j<Ngold_points;j++)
    xgold[j]=rr+j*step;
 if (Rplot_type){
    for (i=1;i<=Ngoldstone;i++){
       rr=pi*i/(xi*Zxi);
       for (j=0;j<Ngold_points;j++)
          ygold[j]=rr/xgold[j];
       goldstone.setup(i,Ngold_points,xgold,ygold,solid_line,black,line_thickness,
                       no_symbol,white,none,0.8,"",0.0,0);
       }
    }
 else{
    for (i=1;i<=Ngoldstone;i++){
       rr=pi*i/(xi*Zxi);
       for (j=0;j<Ngold_points;j++)
          ygold[j]=rr*xgold[j];
       goldstone.setup(i,Ngold_points,xgold,ygold,solid_line,black,line_thickness,
                       no_symbol,white,none,0.8,"",0.0,0);
    }
 }
}
 
// *********************************************************************************

void pump_out_plot(const charString& outfile,
                   int level1, double label1_xloc, int label1_yloc)
{
 int Nlevels=1;
 xydy_curves levels(Nlevels);
 if (levels.setup_from_stream(1,in,channel_match(level1),Rplot_type,
                              no_line,white,line_thickness,circle,blue, 
                              solid_fill,symbol_size,channel_label(level1),
                              label1_xloc,label1_yloc)) return;
 produce_plot(outfile+".agr",title,subtitle,xlabel,ylabel,xmin,xmax,ymin,ymax,
              x_major_tick_spacing,x_minor_tick_num,
              y_major_tick_spacing,y_minor_tick_num,
              levels,goldstone,label_strings);
#ifdef WEBPAGE
 charString str="xmgrace -hdevice EPS -printfile "+outfile+".eps -hardcopy "+outfile+".agr";
 system(convert_to_charstar(str));
 str="rm -f "+outfile+".gif";
 system(convert_to_charstar(str));
 str="ps2gif "+outfile;
 system(convert_to_charstar(str));
 str=outfile;
 while (str.locate('/')>=0){
    str=str.right_relax_substring('/');}
 str=str+".gif";
 html_out << "<IMG SRC=\"" << str <<"\"><br><br><br><br>" <<endl;
#endif
}

//  Parameter:  missing_curve_flag   1 = do not do a plot if any of its curves are missing
//                                   0 = do plot but leave out missing curves 

void pump_out_plot(const charString& outfile,
                   int level1, double label1_xloc, int label1_yloc,
                   int level2, double label2_xloc, int label2_yloc,
                   int missing_curve_flag)
{
 int Nlevels=2;
 int flag1,flag2;
 xydy_curves levels(Nlevels);
 flag1=levels.setup_from_stream(1,in,channel_match(level1),Rplot_type,
                                no_line,white,line_thickness,circle,blue, 
                                solid_fill,symbol_size,channel_label(level1),
                                label1_xloc,label1_yloc);
 flag2=levels.setup_from_stream(2,in,channel_match(level2),Rplot_type,
                                no_line,white,line_thickness,square,red, 
                                solid_fill,symbol_size,channel_label(level2),
                                label2_xloc,label2_yloc);
 Nlevels-=flag1+flag2;
 if ((missing_curve_flag)&&(Nlevels<2)) Nlevels=0;
 if (Nlevels==0){
    return;
    }
 if (flag1){
    pump_out_plot(outfile,level2,label2_xloc,label2_yloc);
    return;
    }
 else if (flag2){
    pump_out_plot(outfile,level1,label1_xloc,label1_yloc);
    return;
    }

 produce_plot(outfile+".agr",title,subtitle,xlabel,ylabel,xmin,xmax,ymin,ymax,
              x_major_tick_spacing,x_minor_tick_num,
              y_major_tick_spacing,y_minor_tick_num,
              levels,goldstone,label_strings);
#ifdef WEBPAGE
 charString str="xmgrace -hdevice EPS -printfile "+outfile+".eps -hardcopy "+outfile+".agr";
 system(convert_to_charstar(str));
 str="rm -f "+outfile+".gif";
 system(convert_to_charstar(str));
 str="ps2gif "+outfile;
 system(convert_to_charstar(str));
 str=outfile;
 while (str.locate('/')>=0){
    str=str.right_relax_substring('/');}
 str=str+".gif";
 html_out << "<IMG SRC=\"" << str <<"\"><br><br><br><br>" <<endl;
#endif
}


void pump_out_plot(const charString& outfile,
                   int level1, double label1_xloc, int label1_yloc,
                   int level2, double label2_xloc, int label2_yloc,
                   int level3, double label3_xloc, int label3_yloc,
                   int missing_curve_flag)
{
 int flag1,flag2,flag3,Nlevels=3;
 xydy_curves levels(Nlevels);
 flag1=levels.setup_from_stream(1,in,channel_match(level1),Rplot_type,
                                no_line,white,line_thickness,circle,blue,
                                solid_fill,symbol_size,channel_label(level1),
                                label1_xloc,label1_yloc);
 flag2=levels.setup_from_stream(2,in,channel_match(level2),Rplot_type,
                                no_line,white,line_thickness,square,red,
                                solid_fill,symbol_size,channel_label(level2),
                                label2_xloc,label2_yloc);
 flag3=levels.setup_from_stream(3,in,channel_match(level3),Rplot_type,
                                no_line,white,line_thickness,triangle_up,green4,
                                solid_fill,symbol_size,channel_label(level3),
                                label3_xloc,label3_yloc);
 Nlevels-=flag1+flag2+flag3;
 if ((missing_curve_flag)&&(Nlevels<3)) Nlevels=0;
 if (Nlevels==0){
     return;
     }
 else if (Nlevels==1){
     if (!flag1) pump_out_plot(outfile,level1,label1_xloc,label1_yloc);
     else if (!flag2) pump_out_plot(outfile,level2,label2_xloc,label2_yloc);
     else pump_out_plot(outfile,level3,label3_xloc,label3_yloc);
     return;
     }
 else if (Nlevels==2){
     if (flag1) pump_out_plot(outfile,level2,label2_xloc,label2_yloc, 
                                      level3,label3_xloc,label3_yloc,missing_curve_flag);
     else if (flag2) pump_out_plot(outfile,level1,label1_xloc,label1_yloc, 
                                           level3,label3_xloc,label3_yloc,missing_curve_flag);
     else pump_out_plot(outfile,level1,label1_xloc,label1_yloc, 
                                level2,label2_xloc,label2_yloc,missing_curve_flag);
     return;
     }

 produce_plot(outfile+".agr",title,subtitle,xlabel,ylabel,xmin,xmax,ymin,ymax,
              x_major_tick_spacing,x_minor_tick_num,
              y_major_tick_spacing,y_minor_tick_num,
              levels,goldstone,label_strings);
#ifdef WEBPAGE
 charString str="xmgrace -hdevice EPS -printfile "+outfile+".eps -hardcopy "+outfile+".agr";
 system(convert_to_charstar(str));
 str="rm -f "+outfile+".gif";
 system(convert_to_charstar(str));
 str="ps2gif "+outfile;
 system(convert_to_charstar(str));
 str=outfile;
 while (str.locate('/')>=0){
    str=str.right_relax_substring('/');}
 str=str+".gif";
 html_out << "<IMG SRC=\"" << str <<"\"><br><br><br><br>" <<endl;
#endif
}


void pump_out_plot(const charString& outfile,
                   int level1, double label1_xloc, int label1_yloc,
                   int level2, double label2_xloc, int label2_yloc,
                   int level3, double label3_xloc, int label3_yloc,
                   int level4, double label4_xloc, int label4_yloc,
                   int missing_curve_flag)
{
 int flag1,flag2,flag3,flag4,Nlevels=4;
 xydy_curves levels(Nlevels);
 flag1=levels.setup_from_stream(1,in,channel_match(level1),Rplot_type,
                                no_line,white,line_thickness,circle,blue,
                                solid_fill,symbol_size,channel_label(level1),
                                label1_xloc,label1_yloc);
 flag2=levels.setup_from_stream(2,in,channel_match(level2),Rplot_type,
                                no_line,white,line_thickness,square,red,
                                solid_fill,symbol_size,channel_label(level2),
                                label2_xloc,label2_yloc);
 flag3=levels.setup_from_stream(3,in,channel_match(level3),Rplot_type,
                                no_line,white,line_thickness,triangle_up,green4,
                                solid_fill,symbol_size,channel_label(level3),
                                label3_xloc,label3_yloc);
 flag4=levels.setup_from_stream(4,in,channel_match(level4),Rplot_type,
                                no_line,white,line_thickness,diamond,magenta,
                                solid_fill,symbol_size,channel_label(level4),
                                label4_xloc,label4_yloc);
 Nlevels-=flag1+flag2+flag3+flag4;
 if ((missing_curve_flag)&&(Nlevels<4)) Nlevels=0;
 if (Nlevels==0){
     return;
     }
 else if (Nlevels==1){
     if (!flag1) pump_out_plot(outfile,level1,label1_xloc,label1_yloc);
     else if (!flag2) pump_out_plot(outfile,level2,label2_xloc,label2_yloc);
     else if (!flag3) pump_out_plot(outfile,level3,label3_xloc,label3_yloc);
     else pump_out_plot(outfile,level4,label4_xloc,label4_yloc);
     return;
     }
 else if (Nlevels==2){
     if ((flag1)&&(flag2)) pump_out_plot(outfile,level3,label3_xloc,label3_yloc,
                                                 level4,label4_xloc,label4_yloc,missing_curve_flag);
     else if ((flag1)&&(flag3)) pump_out_plot(outfile,level2,label2_xloc,label2_yloc,
                                                      level4,label4_xloc,label4_yloc,missing_curve_flag);
     else if ((flag1)&&(flag4)) pump_out_plot(outfile,level2,label2_xloc,label2_yloc,
                                                      level3,label3_xloc,label3_yloc,missing_curve_flag);
     else if ((flag2)&&(flag3)) pump_out_plot(outfile,level1,label1_xloc,label1_yloc,
                                                      level4,label4_xloc,label4_yloc,missing_curve_flag);
     else if ((flag2)&&(flag4)) pump_out_plot(outfile,level1,label1_xloc,label1_yloc,
                                                      level3,label3_xloc,label3_yloc,missing_curve_flag);
     else pump_out_plot(outfile,level1,label1_xloc,label1_yloc,
                                level2,label2_xloc,label2_yloc,missing_curve_flag);
     return;
     }
 else if (Nlevels==3){
     if (flag1) pump_out_plot(outfile,level2,label2_xloc,label2_yloc, 
                                      level3,label3_xloc,label3_yloc,
                                      level4,label4_xloc,label4_yloc,missing_curve_flag);
     else if (flag2) pump_out_plot(outfile,level1,label1_xloc,label1_yloc, 
                                           level3,label3_xloc,label3_yloc,
                                           level4,label4_xloc,label4_yloc,missing_curve_flag);
     else if (flag3) pump_out_plot(outfile,level1,label1_xloc,label1_yloc, 
                                           level2,label2_xloc,label2_yloc,
                                           level4,label4_xloc,label4_yloc,missing_curve_flag);
     else pump_out_plot(outfile,level1,label1_xloc,label1_yloc, 
                                level2,label2_xloc,label2_yloc,
                                level3,label3_xloc,label3_yloc,missing_curve_flag);
     return;
     }
 
 produce_plot(outfile+".agr",title,subtitle,xlabel,ylabel,xmin,xmax,ymin,ymax,
              x_major_tick_spacing,x_minor_tick_num,
              y_major_tick_spacing,y_minor_tick_num,
              levels,goldstone,label_strings);
#ifdef WEBPAGE
 charString str="xmgrace -hdevice EPS -printfile "+outfile+".eps -hardcopy "+outfile+".agr";
 system(convert_to_charstar(str));
 str="rm -f "+outfile+".gif";
 system(convert_to_charstar(str));
 str="ps2gif "+outfile;
 system(convert_to_charstar(str));
 str=outfile;
 while (str.locate('/')>=0){
    str=str.right_relax_substring('/');}
 str=str+".gif";
 html_out << "<IMG SRC=\"" << str <<"\"><br><br><br><br>" <<endl;
#endif
}

 // ***************************************************************************

void html_header(const charString filename, const charString& up_filename)
{
 html_out.open(convert_to_charstar(filename));
 if (!html_out){
     cerr << " could not open html file for output"<<endl;
     exit(1);}

 html_out << "<HTML>"<<endl;
 html_out << "<HEAD>"<<endl;
 html_out << "<TITLE> Plots</TITLE>"<<endl;
 html_out << "</HEAD>"<<endl;
 html_out << "<BODY BGCOLOR=\"#99FFCC\">"<<endl;
 html_out << "<a href=\"" << up_filename << "\">"<<endl;
 html_out << "<img src=\"../../../images/back.gif\" align=middle border=0></a>"<<endl;
 html_out << "<center>"<<endl;
}


void html_footer(const charString& up_filename)
{
 html_out << "</center>"<<endl;
 html_out << "<BR><BR><BR>"<<endl;
 html_out << "<a href=\"" << up_filename << "\">"<<endl;
 html_out << "<img src=\"../../../images/back.gif\" align=middle border=0></a>"<<endl;
 html_out << "<a href=\"#top\">"<<endl;
 html_out << "<img src=\"../../../images/to_top.gif\" align=middle border=0></a>"<<endl;
 html_out << "</BODY>"<<endl;
 html_out << "</HTML>"<<endl;
 html_out.close();
}

// *******************************************************************************

void pump_out_string_plot(const charString& outfile,
                          int level1, double label1_xloc, int label1_yloc,
                          int string_level, double xi, double Zxi, double Zxi_err)
{
 const double pi=3.141592653589793238;
 int i,Nlevels=1;
 double y,dy,r;
 xydy_curves levels(Nlevels);
 if (levels.setup_from_stream(1,in,channel_match(level1),1,
                              no_line,white,line_thickness,circle,blue, 
                              solid_fill,symbol_size,channel_label(level1),
                              label1_xloc,label1_yloc)) return;
 for (i=1;i<=levels(1).Npoints;i++){
    r=levels(1).xvals(i)*xi/(pi*string_level);
    y=levels(1).yvals(i);
    dy=levels(1).dyvals(i);
    levels(1).dyvals(i)=r*sqrt(dy*dy*Zxi*Zxi+Zxi_err*Zxi_err*y*y);
    levels(1).yvals(i)=y*Zxi*r-1.0;
    }
 double xx[2],yy[2];
 xx[0]=xmin; xx[1]=xmax;
 yy[0]=yy[1]=0.0;
 zero_line.setup(1,2,xx,yy,solid_line,black,2*line_thickness,
                       no_symbol,white,none,0.8,"",0.0,0);
 produce_plot(outfile+".agr",title,subtitle,xlabel,ylabel,xmin,xmax,ymin,ymax,
              x_major_tick_spacing,x_minor_tick_num,
              y_major_tick_spacing,y_minor_tick_num,
              levels,zero_line,label_strings);
#ifdef WEBPAGE
 charString str="xmgrace -hdevice EPS -printfile "+outfile+".eps -hardcopy "+outfile+".agr";
 system(convert_to_charstar(str));
 str="rm -f "+outfile+".gif";
 system(convert_to_charstar(str));
 str="ps2gif "+outfile;
 system(convert_to_charstar(str));
 str=outfile;
 while (str.locate('/')>=0){
    str=str.right_relax_substring('/');}
 str=str+".gif";
 html_out << "<IMG SRC=\"" << str <<"\"><br><br><br><br>" <<endl;
#endif
}

// *********************************************************************************