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


//#define XERRORS
#define COMPACT_LEGEND

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

// useful global variables

const double symbol_size=0.6;
int line_thickness=2,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);
double legend_xpos,legend_ypos;

int Nchannels;
Svector channel_match,channel_label;
ifstream in;
ofstream html_out;
xy_curves goldstone;
xy_curves comparison_curves(3);


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

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]);}

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

level_info::level_info()
{Ncurves=0;}

void level_info::assign(int string_N,  const charString& rundir, 
     const charString& stub1, const charString& patt1, const charString &legend1)
{
 int count=0;
 string_level=string_N;
 Ncurves=1;
 datafiles.reset(Ncurves);
 gapdatafiles.reset(Ncurves);
 pattern_search.reset(Ncurves);
 legend_strings.reset(Ncurves);
 count++;
 datafiles(count)=rundir+stub1+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub1+"/scaled_gap_results.dat";
 pattern_search(count)=patt1;
 legend_strings(count)=legend1;
}

void level_info::assign(int string_N,  const charString& rundir, 
     const charString& stub1, const charString& patt1, const charString &legend1,
     const charString& stub2, const charString& patt2, const charString &legend2)
{
 int count=0;
 string_level=string_N;
 Ncurves=2;
 datafiles.reset(Ncurves);
 gapdatafiles.reset(Ncurves);
 pattern_search.reset(Ncurves);
 legend_strings.reset(Ncurves);
 count++;
 datafiles(count)=rundir+stub1+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub1+"/scaled_gap_results.dat";
 pattern_search(count)=patt1;
 legend_strings(count)=legend1;
 count++;
 datafiles(count)=rundir+stub2+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub2+"/scaled_gap_results.dat";
 pattern_search(count)=patt2;
 legend_strings(count)=legend2;
}

void level_info::assign(int string_N,  const charString& rundir, 
     const charString& stub1, const charString& patt1, const charString &legend1,
     const charString& stub2, const charString& patt2, const charString &legend2,
     const charString& stub3, const charString& patt3, const charString &legend3)
{
 int count=0;
 string_level=string_N;
 Ncurves=3;
 datafiles.reset(Ncurves);
 gapdatafiles.reset(Ncurves);
 pattern_search.reset(Ncurves);
 legend_strings.reset(Ncurves);
 count++;
 datafiles(count)=rundir+stub1+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub1+"/scaled_gap_results.dat";
 pattern_search(count)=patt1;
 legend_strings(count)=legend1;
 count++;
 datafiles(count)=rundir+stub2+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub2+"/scaled_gap_results.dat";
 pattern_search(count)=patt2;
 legend_strings(count)=legend2;
 count++;
 datafiles(count)=rundir+stub3+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub3+"/scaled_gap_results.dat";
 pattern_search(count)=patt3;
 legend_strings(count)=legend3;
}

void level_info::assign(int string_N,  const charString& rundir, 
     const charString& stub1, const charString& patt1, const charString &legend1,
     const charString& stub2, const charString& patt2, const charString &legend2,
     const charString& stub3, const charString& patt3, const charString &legend3,
     const charString& stub4, const charString& patt4, const charString &legend4)
{
 int count=0;
 string_level=string_N;
 Ncurves=4;
 datafiles.reset(Ncurves);
 gapdatafiles.reset(Ncurves);
 pattern_search.reset(Ncurves);
 legend_strings.reset(Ncurves);
 count++;
 datafiles(count)=rundir+stub1+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub1+"/scaled_gap_results.dat";
 pattern_search(count)=patt1;
 legend_strings(count)=legend1;
 count++;
 datafiles(count)=rundir+stub2+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub2+"/scaled_gap_results.dat";
 pattern_search(count)=patt2;
 legend_strings(count)=legend2;
 count++;
 datafiles(count)=rundir+stub3+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub3+"/scaled_gap_results.dat";
 pattern_search(count)=patt3;
 legend_strings(count)=legend3;
 count++;
 datafiles(count)=rundir+stub4+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub4+"/scaled_gap_results.dat";
 pattern_search(count)=patt4;
 legend_strings(count)=legend4;
}

void level_info::assign(int string_N,  const charString& rundir, 
     const charString& stub1, const charString& patt1, const charString &legend1,
     const charString& stub2, const charString& patt2, const charString &legend2,
     const charString& stub3, const charString& patt3, const charString &legend3,
     const charString& stub4, const charString& patt4, const charString &legend4,
     const charString& stub5, const charString& patt5, const charString &legend5)
{
 int count=0;
 string_level=string_N;
 Ncurves=5;
 datafiles.reset(Ncurves);
 gapdatafiles.reset(Ncurves);
 pattern_search.reset(Ncurves);
 legend_strings.reset(Ncurves);
 count++;
 datafiles(count)=rundir+stub1+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub1+"/scaled_gap_results.dat";
 pattern_search(count)=patt1;
 legend_strings(count)=legend1;
 count++;
 datafiles(count)=rundir+stub2+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub2+"/scaled_gap_results.dat";
 pattern_search(count)=patt2;
 legend_strings(count)=legend2;
 count++;
 datafiles(count)=rundir+stub3+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub3+"/scaled_gap_results.dat";
 pattern_search(count)=patt3;
 legend_strings(count)=legend3;
 count++;
 datafiles(count)=rundir+stub4+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub4+"/scaled_gap_results.dat";
 pattern_search(count)=patt4;
 legend_strings(count)=legend4;
 count++;
 datafiles(count)=rundir+stub5+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub5+"/scaled_gap_results.dat";
 pattern_search(count)=patt5;
 legend_strings(count)=legend5;
}

void level_info::assign(int string_N,  const charString& rundir, 
     const charString& stub1, const charString& patt1, const charString &legend1,
     const charString& stub2, const charString& patt2, const charString &legend2,
     const charString& stub3, const charString& patt3, const charString &legend3,
     const charString& stub4, const charString& patt4, const charString &legend4,
     const charString& stub5, const charString& patt5, const charString &legend5,
     const charString& stub6, const charString& patt6, const charString &legend6)
{
 int count=0;
 string_level=string_N;
 Ncurves=6;
 datafiles.reset(Ncurves);
 gapdatafiles.reset(Ncurves);
 pattern_search.reset(Ncurves);
 legend_strings.reset(Ncurves);
 count++;
 datafiles(count)=rundir+stub1+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub1+"/scaled_gap_results.dat";
 pattern_search(count)=patt1;
 legend_strings(count)=legend1;
 count++;
 datafiles(count)=rundir+stub2+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub2+"/scaled_gap_results.dat";
 pattern_search(count)=patt2;
 legend_strings(count)=legend2;
 count++;
 datafiles(count)=rundir+stub3+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub3+"/scaled_gap_results.dat";
 pattern_search(count)=patt3;
 legend_strings(count)=legend3;
 count++;
 datafiles(count)=rundir+stub4+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub4+"/scaled_gap_results.dat";
 pattern_search(count)=patt4;
 legend_strings(count)=legend4;
 count++;
 datafiles(count)=rundir+stub5+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub5+"/scaled_gap_results.dat";
 pattern_search(count)=patt5;
 legend_strings(count)=legend5;
 count++;
 datafiles(count)=rundir+stub6+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub6+"/scaled_gap_results.dat";
 pattern_search(count)=patt6;
 legend_strings(count)=legend6;
}

void level_info::assign(int string_N,  const charString& rundir, 
     const charString& stub1, const charString& patt1, const charString &legend1,
     const charString& stub2, const charString& patt2, const charString &legend2,
     const charString& stub3, const charString& patt3, const charString &legend3,
     const charString& stub4, const charString& patt4, const charString &legend4,
     const charString& stub5, const charString& patt5, const charString &legend5,
     const charString& stub6, const charString& patt6, const charString &legend6,
     const charString& stub7, const charString& patt7, const charString &legend7)
{
 int count=0;
 string_level=string_N;
 Ncurves=7;
 datafiles.reset(Ncurves);
 gapdatafiles.reset(Ncurves);
 pattern_search.reset(Ncurves);
 legend_strings.reset(Ncurves);
 count++;
 datafiles(count)=rundir+stub1+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub1+"/scaled_gap_results.dat";
 pattern_search(count)=patt1;
 legend_strings(count)=legend1;
 count++;
 datafiles(count)=rundir+stub2+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub2+"/scaled_gap_results.dat";
 pattern_search(count)=patt2;
 legend_strings(count)=legend2;
 count++;
 datafiles(count)=rundir+stub3+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub3+"/scaled_gap_results.dat";
 pattern_search(count)=patt3;
 legend_strings(count)=legend3;
 count++;
 datafiles(count)=rundir+stub4+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub4+"/scaled_gap_results.dat";
 pattern_search(count)=patt4;
 legend_strings(count)=legend4;
 count++;
 datafiles(count)=rundir+stub5+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub5+"/scaled_gap_results.dat";
 pattern_search(count)=patt5;
 legend_strings(count)=legend5;
 count++;
 datafiles(count)=rundir+stub6+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub6+"/scaled_gap_results.dat";
 pattern_search(count)=patt6;
 legend_strings(count)=legend6;
 count++;
 datafiles(count)=rundir+stub7+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub7+"/scaled_gap_results.dat";
 pattern_search(count)=patt7;
 legend_strings(count)=legend7;
}

void level_info::assign(int string_N,  const charString& rundir, 
     const charString& stub1, const charString& patt1, const charString &legend1,
     const charString& stub2, const charString& patt2, const charString &legend2,
     const charString& stub3, const charString& patt3, const charString &legend3,
     const charString& stub4, const charString& patt4, const charString &legend4,
     const charString& stub5, const charString& patt5, const charString &legend5,
     const charString& stub6, const charString& patt6, const charString &legend6,
     const charString& stub7, const charString& patt7, const charString &legend7,
     const charString& stub8, const charString& patt8, const charString &legend8)
{
 int count=0;
 string_level=string_N;
 Ncurves=8;
 datafiles.reset(Ncurves);
 gapdatafiles.reset(Ncurves);
 pattern_search.reset(Ncurves);
 legend_strings.reset(Ncurves);
 count++;
 datafiles(count)=rundir+stub1+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub1+"/scaled_gap_results.dat";
 pattern_search(count)=patt1;
 legend_strings(count)=legend1;
 count++;
 datafiles(count)=rundir+stub2+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub2+"/scaled_gap_results.dat";
 pattern_search(count)=patt2;
 legend_strings(count)=legend2;
 count++;
 datafiles(count)=rundir+stub3+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub3+"/scaled_gap_results.dat";
 pattern_search(count)=patt3;
 legend_strings(count)=legend3;
 count++;
 datafiles(count)=rundir+stub4+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub4+"/scaled_gap_results.dat";
 pattern_search(count)=patt4;
 legend_strings(count)=legend4;
 count++;
 datafiles(count)=rundir+stub5+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub5+"/scaled_gap_results.dat";
 pattern_search(count)=patt5;
 legend_strings(count)=legend5;
 count++;
 datafiles(count)=rundir+stub6+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub6+"/scaled_gap_results.dat";
 pattern_search(count)=patt6;
 legend_strings(count)=legend6;
 count++;
 datafiles(count)=rundir+stub7+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub7+"/scaled_gap_results.dat";
 pattern_search(count)=patt7;
 legend_strings(count)=legend7;
 count++;
 datafiles(count)=rundir+stub8+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub8+"/scaled_gap_results.dat";
 pattern_search(count)=patt8;
 legend_strings(count)=legend8;
}

void level_info::assign(int string_N,  const charString& rundir, 
     const charString& stub1, const charString& patt1, const charString &legend1,
     const charString& stub2, const charString& patt2, const charString &legend2,
     const charString& stub3, const charString& patt3, const charString &legend3,
     const charString& stub4, const charString& patt4, const charString &legend4,
     const charString& stub5, const charString& patt5, const charString &legend5,
     const charString& stub6, const charString& patt6, const charString &legend6,
     const charString& stub7, const charString& patt7, const charString &legend7,
     const charString& stub8, const charString& patt8, const charString &legend8,
     const charString& stub9, const charString& patt9, const charString &legend9)
{
 int count=0;
 string_level=string_N;
 Ncurves=9;
 datafiles.reset(Ncurves);
 gapdatafiles.reset(Ncurves);
 pattern_search.reset(Ncurves);
 legend_strings.reset(Ncurves);
 count++;
 datafiles(count)=rundir+stub1+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub1+"/scaled_gap_results.dat";
 pattern_search(count)=patt1;
 legend_strings(count)=legend1;
 count++;
 datafiles(count)=rundir+stub2+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub2+"/scaled_gap_results.dat";
 pattern_search(count)=patt2;
 legend_strings(count)=legend2;
 count++;
 datafiles(count)=rundir+stub3+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub3+"/scaled_gap_results.dat";
 pattern_search(count)=patt3;
 legend_strings(count)=legend3;
 count++;
 datafiles(count)=rundir+stub4+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub4+"/scaled_gap_results.dat";
 pattern_search(count)=patt4;
 legend_strings(count)=legend4;
 count++;
 datafiles(count)=rundir+stub5+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub5+"/scaled_gap_results.dat";
 pattern_search(count)=patt5;
 legend_strings(count)=legend5;
 count++;
 datafiles(count)=rundir+stub6+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub6+"/scaled_gap_results.dat";
 pattern_search(count)=patt6;
 legend_strings(count)=legend6;
 count++;
 datafiles(count)=rundir+stub7+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub7+"/scaled_gap_results.dat";
 pattern_search(count)=patt7;
 legend_strings(count)=legend7;
 count++;
 datafiles(count)=rundir+stub8+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub8+"/scaled_gap_results.dat";
 pattern_search(count)=patt8;
 legend_strings(count)=legend8;
 count++;
 datafiles(count)=rundir+stub9+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub9+"/scaled_gap_results.dat";
 pattern_search(count)=patt9;
 legend_strings(count)=legend9;
}

void level_info::assign(int string_N,  const charString& rundir, 
     const charString& stub1, const charString& patt1, const charString &legend1,
     const charString& stub2, const charString& patt2, const charString &legend2,
     const charString& stub3, const charString& patt3, const charString &legend3,
     const charString& stub4, const charString& patt4, const charString &legend4,
     const charString& stub5, const charString& patt5, const charString &legend5,
     const charString& stub6, const charString& patt6, const charString &legend6,
     const charString& stub7, const charString& patt7, const charString &legend7,
     const charString& stub8, const charString& patt8, const charString &legend8,
     const charString& stub9, const charString& patt9, const charString &legend9,
     const charString& stub10,const charString& patt10,const charString &legend10)
{
 int count=0;
 string_level=string_N;
 Ncurves=10;
 datafiles.reset(Ncurves);
 gapdatafiles.reset(Ncurves);
 pattern_search.reset(Ncurves);
 legend_strings.reset(Ncurves);
 count++;
 datafiles(count)=rundir+stub1+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub1+"/scaled_gap_results.dat";
 pattern_search(count)=patt1;
 legend_strings(count)=legend1;
 count++;
 datafiles(count)=rundir+stub2+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub2+"/scaled_gap_results.dat";
 pattern_search(count)=patt2;
 legend_strings(count)=legend2;
 count++;
 datafiles(count)=rundir+stub3+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub3+"/scaled_gap_results.dat";
 pattern_search(count)=patt3;
 legend_strings(count)=legend3;
 count++;
 datafiles(count)=rundir+stub4+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub4+"/scaled_gap_results.dat";
 pattern_search(count)=patt4;
 legend_strings(count)=legend4;
 count++;
 datafiles(count)=rundir+stub5+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub5+"/scaled_gap_results.dat";
 pattern_search(count)=patt5;
 legend_strings(count)=legend5;
 count++;
 datafiles(count)=rundir+stub6+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub6+"/scaled_gap_results.dat";
 pattern_search(count)=patt6;
 legend_strings(count)=legend6;
 count++;
 datafiles(count)=rundir+stub7+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub7+"/scaled_gap_results.dat";
 pattern_search(count)=patt7;
 legend_strings(count)=legend7;
 count++;
 datafiles(count)=rundir+stub8+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub8+"/scaled_gap_results.dat";
 pattern_search(count)=patt8;
 legend_strings(count)=legend8;
 count++;
 datafiles(count)=rundir+stub9+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub9+"/scaled_gap_results.dat";
 pattern_search(count)=patt9;
 legend_strings(count)=legend9;
 count++;
 datafiles(count)=rundir+stub10+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub10+"/scaled_gap_results.dat";
 pattern_search(count)=patt10;
 legend_strings(count)=legend10;
}

void level_info::assign(int string_N,  const charString& rundir, 
     const charString& stub1, const charString& patt1, const charString &legend1,
     const charString& stub2, const charString& patt2, const charString &legend2,
     const charString& stub3, const charString& patt3, const charString &legend3,
     const charString& stub4, const charString& patt4, const charString &legend4,
     const charString& stub5, const charString& patt5, const charString &legend5,
     const charString& stub6, const charString& patt6, const charString &legend6,
     const charString& stub7, const charString& patt7, const charString &legend7,
     const charString& stub8, const charString& patt8, const charString &legend8,
     const charString& stub9, const charString& patt9, const charString &legend9,
     const charString& stub10,const charString& patt10,const charString &legend10,
     const charString& stub11,const charString& patt11,const charString &legend11)
{
 int count=0;
 string_level=string_N;
 Ncurves=11;
 datafiles.reset(Ncurves);
 gapdatafiles.reset(Ncurves);
 pattern_search.reset(Ncurves);
 legend_strings.reset(Ncurves);
 count++;
 datafiles(count)=rundir+stub1+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub1+"/scaled_gap_results.dat";
 pattern_search(count)=patt1;
 legend_strings(count)=legend1;
 count++;
 datafiles(count)=rundir+stub2+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub2+"/scaled_gap_results.dat";
 pattern_search(count)=patt2;
 legend_strings(count)=legend2;
 count++;
 datafiles(count)=rundir+stub3+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub3+"/scaled_gap_results.dat";
 pattern_search(count)=patt3;
 legend_strings(count)=legend3;
 count++;
 datafiles(count)=rundir+stub4+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub4+"/scaled_gap_results.dat";
 pattern_search(count)=patt4;
 legend_strings(count)=legend4;
 count++;
 datafiles(count)=rundir+stub5+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub5+"/scaled_gap_results.dat";
 pattern_search(count)=patt5;
 legend_strings(count)=legend5;
 count++;
 datafiles(count)=rundir+stub6+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub6+"/scaled_gap_results.dat";
 pattern_search(count)=patt6;
 legend_strings(count)=legend6;
 count++;
 datafiles(count)=rundir+stub7+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub7+"/scaled_gap_results.dat";
 pattern_search(count)=patt7;
 legend_strings(count)=legend7;
 count++;
 datafiles(count)=rundir+stub8+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub8+"/scaled_gap_results.dat";
 pattern_search(count)=patt8;
 legend_strings(count)=legend8;
 count++;
 datafiles(count)=rundir+stub9+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub9+"/scaled_gap_results.dat";
 pattern_search(count)=patt9;
 legend_strings(count)=legend9;
 count++;
 datafiles(count)=rundir+stub10+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub10+"/scaled_gap_results.dat";
 pattern_search(count)=patt10;
 legend_strings(count)=legend10;
 count++;
 datafiles(count)=rundir+stub11+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub11+"/scaled_gap_results.dat";
 pattern_search(count)=patt11;
 legend_strings(count)=legend11;
}

void level_info::assign(int string_N,  const charString& rundir, 
     const charString& stub1, const charString& patt1, const charString &legend1,
     const charString& stub2, const charString& patt2, const charString &legend2,
     const charString& stub3, const charString& patt3, const charString &legend3,
     const charString& stub4, const charString& patt4, const charString &legend4,
     const charString& stub5, const charString& patt5, const charString &legend5,
     const charString& stub6, const charString& patt6, const charString &legend6,
     const charString& stub7, const charString& patt7, const charString &legend7,
     const charString& stub8, const charString& patt8, const charString &legend8,
     const charString& stub9, const charString& patt9, const charString &legend9,
     const charString& stub10,const charString& patt10,const charString &legend10,
     const charString& stub11,const charString& patt11,const charString &legend11,
     const charString& stub12,const charString& patt12,const charString &legend12)
{
 int count=0;
 string_level=string_N;
 Ncurves=12;
 datafiles.reset(Ncurves);
 gapdatafiles.reset(Ncurves);
 pattern_search.reset(Ncurves);
 legend_strings.reset(Ncurves);
 count++;
 datafiles(count)=rundir+stub1+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub1+"/scaled_gap_results.dat";
 pattern_search(count)=patt1;
 legend_strings(count)=legend1;
 count++;
 datafiles(count)=rundir+stub2+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub2+"/scaled_gap_results.dat";
 pattern_search(count)=patt2;
 legend_strings(count)=legend2;
 count++;
 datafiles(count)=rundir+stub3+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub3+"/scaled_gap_results.dat";
 pattern_search(count)=patt3;
 legend_strings(count)=legend3;
 count++;
 datafiles(count)=rundir+stub4+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub4+"/scaled_gap_results.dat";
 pattern_search(count)=patt4;
 legend_strings(count)=legend4;
 count++;
 datafiles(count)=rundir+stub5+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub5+"/scaled_gap_results.dat";
 pattern_search(count)=patt5;
 legend_strings(count)=legend5;
 count++;
 datafiles(count)=rundir+stub6+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub6+"/scaled_gap_results.dat";
 pattern_search(count)=patt6;
 legend_strings(count)=legend6;
 count++;
 datafiles(count)=rundir+stub7+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub7+"/scaled_gap_results.dat";
 pattern_search(count)=patt7;
 legend_strings(count)=legend7;
 count++;
 datafiles(count)=rundir+stub8+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub8+"/scaled_gap_results.dat";
 pattern_search(count)=patt8;
 legend_strings(count)=legend8;
 count++;
 datafiles(count)=rundir+stub9+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub9+"/scaled_gap_results.dat";
 pattern_search(count)=patt9;
 legend_strings(count)=legend9;
 count++;
 datafiles(count)=rundir+stub10+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub10+"/scaled_gap_results.dat";
 pattern_search(count)=patt10;
 legend_strings(count)=legend10;
 count++;
 datafiles(count)=rundir+stub11+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub11+"/scaled_gap_results.dat";
 pattern_search(count)=patt11;
 legend_strings(count)=legend11;
 count++;
 datafiles(count)=rundir+stub12+"/scaled_results.dat"; 
 gapdatafiles(count)=rundir+stub12+"/scaled_gap_results.dat";
 pattern_search(count)=patt12;
 legend_strings(count)=legend12;
}

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

   //  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
   //  or 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 Rerr_value;
 double Rinv_value;
 double Rinverr_value;
 double mean_value;
 double stdev_value;
 node *next;
};

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

void add_node(double R, double Rerr, double Rinv, double Rinverr, double ans, double err)
{
 if (head==NULL){
    head=new node;
    head->next=NULL;
    head->R_value=R;
    head->Rinv_value=Rinv;
    head->Rerr_value=Rerr;
    head->Rinverr_value=Rinverr;
    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->Rerr_value=Rerr;
    tail->Rinverr_value=Rinverr;
    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 plot_curve::setup_from_file(
          const charString& datafile, 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)
{
 const int maxlinelength=132;
 static char ff[maxlinelength];
 Svector tokens;
 int flag;
 double Rvalue,Rerrvalue,Rinvvalue,Rinverrvalue,ans,err;

 ifstream in(convert_to_charstar(datafile));
 if (!in){
    cerr << " could not open data file";
    in.close();
    return 1;}
 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()==11){
       if (tokens(9).locate(channel_string)>=0) flag=1;
       }
    if (in.eof()){
       cerr << " could not find requested channel"<<endl;
       wipeout();
       in.close();
       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()==7){
    if (string_to_double(tokens(2),Rvalue)) break;
    if (string_to_double(tokens(3),Rerrvalue)) break;
    if (string_to_double(tokens(4),Rinvvalue)) break;
    if (string_to_double(tokens(5),Rinverrvalue)) break;
    if (string_to_double(tokens(6),ans)) break;
    if (string_to_double(tokens(7),err)) break;
    Npoints++;
    add_node(Rvalue,Rerrvalue,Rinvvalue,Rinverrvalue,ans,err);
    in.getline(ff,maxlinelength);
    get_tokens(ff,tokens);
    }
 xvals.reset(Npoints);
 yvals.reset(Npoints);
 dxvals.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;
      dxvals(flag)=pnode->Rerr_value;}
   else{
      xvals(flag)=pnode->Rinv_value;
      dxvals(flag)=pnode->Rinverr_value;}
   pnode=pnode->next;
   }
 delete_nodes();
 in.close();
 
 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;
 return 0;
}


void plot_curve::wipeout()
{
 Npoints=0;
 xvals();
 dxvals();
 yvals();
 dyvals();
 label="";
}


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

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

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

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


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

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


int plot_curves::setup_from_file(
            int index,  const charString& datafile, 
            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)
{
 if ((index>=1)&&(index<=n)){
    return m[index].setup_from_file(datafile,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);
    }
 else{
    cerr << " bad index; no assignment made in plot_curves"<<endl;
    return 1;
    }
}

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

plot_curve &plot_curves::operator()(int i) const
 {return ((i>=1)&&(i<=n)) ? m[i]
    : (cerr << "plot_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)
{
 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;
}


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)
{
 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="";
}

void xy_curve::wipeout()
{
 Npoints=0;
 xvals();
 yvals();
}


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

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)
{
 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);
    }
 else{
    cerr << " bad index; no assignment made in xy_curves"<<endl;
    }
}

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)
{
 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);
    }
 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]);}

 // *******************************************************************
 //
 //   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 plot_curves& plot_results, 
                  const xy_curves& xy_curves_to_plot,
                  const plot_labels& label_strings,
                  double legend_xloc, double legend_yloc)
{
 int i,j;
 charString str;
 
   // 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;

 out << "@ legend on"<<endl;
 out << "@ legend loctype view"<<endl;
 out << "@ legend "<<legend_xloc<<", "<<legend_yloc<<endl;
 out << "@ legend box color 1"<<endl;
 out << "@ legend box pattern 1"<<endl;
 out << "@ legend box linewidth 1.0"<<endl;
 out << "@ legend box linestyle 1"<<endl;
 out << "@ legend box fill color 0"<<endl;
 out << "@ legend box fill pattern 1"<<endl;
 out << "@ legend font 0"<<endl;
#ifdef COMPACT_LEGEND
 out << "@ legend char size 0.88"<<endl;
 out << "@ legend vgap 0"<<endl;
#else
 out << "@ legend char size 1.0"<<endl;
 out << "@ legend vgap 1"<<endl;
#endif
 out << "@ legend color 1"<<endl;
 out << "@ legend length 4"<<endl;
 out << "@ legend hgap 1"<<endl;
 out << "@ legend invert false"<<endl;

 for (i=1;i<=plot_results.dim();i++){
    if (plot_results(i).Npoints>0){
    
    str="@ s"+int_to_string(i)+" ";
    out << str << "on"<<endl;
#ifdef XERRORS
    out << str << "type xydxdy"<<endl;
#else
    out << str << "type xydy"<<endl;
#endif
    out << str << "symbol "<<plot_results(i).symbol_type<<endl;
    out << str << "symbol size "<<plot_results(i).symbol_size<<endl;
    out << str << "symbol color "<<plot_results(i).symbol_color<<endl;
    out << str << "symbol pattern 1"<<endl;
    out << str << "symbol fill color "<<plot_results(i).symbol_color<<endl;
    out << str << "symbol fill pattern "<<plot_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 "<<plot_results(i).line_type<<endl;
    out << str << "line linewidth "<<plot_results(i).line_thickness<<endl;
    out << str << "line color "<<plot_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 "<<plot_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 << str << "legend \""<<plot_results(i).label<<"\""<<endl;

    out << "@target g0.s"<<i<<endl;
#ifdef XERRORS
    out << "@type xydxdy"<<endl;
#else
    out << "@type xydy"<<endl;
#endif
    for (j=1;j<=plot_results(i).Npoints;j++){
       out <<" "<<plot_results(i).xvals(j)
           <<" "<<plot_results(i).yvals(j)
#ifdef XERRORS
           <<" "<<plot_results(i).dxvals(j)
#endif
           <<" "<<plot_results(i).dyvals(j)<<endl;
       }
    out << "&"<<endl;

    }}

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

    str="@ s"+int_to_string(i+plot_results.dim())+" ";
    out << str << "on"<<endl;
    out << str << "type xy"<<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;
    if (xy_curves_to_plot(i).label.length()>0)
       out << str << "legend \""<<xy_curves_to_plot(i).label<<"\""<<endl;

    out << "@target g0.s"<<i+plot_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;

    }

 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 plot_curves& plot_results,
                  const plot_labels& label_strings,
                  double legend_xloc, double legend_yloc)
{
 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,
              plot_results,xy_curves_to_plot,label_strings,
              legend_xloc,legend_yloc);
}

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,
                  double legend_xloc, double legend_yloc)
{
 plot_curves plot_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,
              plot_results,xy_curves_to_plot,label_strings,
              legend_xloc,legend_yloc);
}


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

void setup_goldstone(int Rplot_type, double xmin, double xmax, int string_level, int dim)
{
 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;
 double rr,An,A0,En,E0;
 int j;
 
 goldstone.reset(3);
 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){
    rr=pi*string_level;
    for (j=0;j<Ngold_points;j++)
       ygold[j]=rr/xgold[j];
    goldstone.setup(1,Ngold_points,xgold,ygold,solid_line,black,line_thickness,
                    no_symbol,white,none,0.8,"Goldstone");

    An=pi*(2.0*string_level-(dim-2.0)/12.0);
    A0=-pi*(dim-2.0)/12.0;
    rr=sqrt(-A0)+0.01; // minimum value allowed
    if (xmin>rr) rr=xmin;
    step=(xmax-rr)/Ngold_points;
    for (j=0;j<Ngold_points;j++)
       xgold[j]=rr+j*step;
    for (j=0;j<Ngold_points;j++){
       rr=xgold[j]*xgold[j]; 
       En=sqrt(1.0+An/rr); 
       E0=sqrt(1.0+A0/rr);
       ygold[j]=xgold[j]*(En-E0);
       }
    goldstone.setup(2,Ngold_points,xgold,ygold,dashed_line,blue,line_thickness,
                    no_symbol,white,none,0.8,"Full Nambu-Goto");
    A0=string_level*pi;
    An=0.5*pi*((dim-2.0)/12.0-string_level);
    for (j=0;j<Ngold_points;j++){
       ygold[j]=A0/xgold[j]*(1.0+An/(xgold[j]*xgold[j]));
       }
    goldstone.setup(3,Ngold_points,xgold,ygold,dot_dashed_line,red,line_thickness,
                    no_symbol,white,none,0.8,"One-loop Nambu-Goto");
    }
 else{
    rr=pi*string_level;
    for (j=0;j<Ngold_points;j++)
       ygold[j]=rr*xgold[j];
    goldstone.setup(1,Ngold_points,xgold,ygold,solid_line,black,line_thickness,
                    no_symbol,white,none,0.8,"Goldstone");
    An=pi*(2*string_level-(dim-2.0)/12.0);
    A0=-pi*(dim-2.0)/12.0;
    En=1.0/sqrt(-A0)-0.04; // maximum value allowed
    if (xmax<En) En=xmax;
    rr=xmin; if (rr<0.01) rr=0.01;  // avoid division by zero
    step=(En-rr)/Ngold_points;
    for (j=0;j<Ngold_points;j++){
       xgold[j]=rr+j*step;}
    for (j=0;j<Ngold_points;j++){
       rr=xgold[j]*xgold[j];
       En=sqrt(1.0+An*rr);
       E0=sqrt(1.0+A0*rr);
       ygold[j]=(En-E0)/xgold[j];}
    goldstone.setup(2,Ngold_points,xgold,ygold,dashed_line,blue,line_thickness,
                    no_symbol,white,none,0.8,"Full Nambu-Goto");
    A0=string_level*pi;
    An=0.5*pi*((dim-2.0)/12.0-string_level);
    for (j=0;j<Ngold_points;j++){
       ygold[j]=A0*xgold[j]*(1.0+An*xgold[j]*xgold[j]);
       }
    goldstone.setup(3,Ngold_points,xgold,ygold,dot_dashed_line,red,line_thickness,
                    no_symbol,white,none,0.8,"One-loop Nambu-Goto");
 }
}

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

void setup_full_nambu_goto(double xmin, double xmax, int string_level, int dim)
{
 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;
 double rr,An,A0,En,E0;
 int j;
 
 // comparison_curves(2) = one_loop NambuGoto
 // comparison_curves(3) = full NambuGoto

 rr=xmin;
 if (fabs(rr)<eps) rr+=step;
 for (j=0;j<Ngold_points;j++)
    xgold[j]=rr+j*step;
 An=pi*(2*string_level-(dim-2.0)/12.0);
 A0=-pi*(dim-2.0)/12.0;
 for (j=0;j<Ngold_points;j++){
    rr=xgold[j]*xgold[j];
    En=sqrt(1.0+An/rr);
    E0=sqrt(1.0+A0/rr);
    ygold[j]=rr*(En-E0)/(string_level*pi)-1.0;}
 comparison_curves.setup(3,Ngold_points,xgold,ygold,dashed_line,blue,line_thickness,
                         no_symbol,white,none,0.8,"Full Nambu-Goto");
}

void setup_one_loop_nambu_goto(double xmin, double xmax, int string_level, int dim)
{
 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;
 double rr,A0;
 int j;
 
 // comparison_curves(2) = one_loop NambuGoto
 // comparison_curves(3) = full NambuGoto

 rr=xmin;
 if (fabs(rr)<eps) rr+=step;
 for (j=0;j<Ngold_points;j++)
    xgold[j]=rr+j*step;
 A0=-pi/2.0*(string_level-1.0/12.0);
 for (j=0;j<Ngold_points;j++){
    rr=xgold[j]*xgold[j];
    ygold[j]=A0/rr;}
 comparison_curves.setup(2,Ngold_points,xgold,ygold,dot_dashed_line,red,line_thickness,
                         no_symbol,white,none,0.8,"One-loop Nambu-Goto");
}

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

void pump_out_plot(const charString& outfile, int Nlevels,
                   const Svector& datafiles, const Svector& channel_match,
                   const Svector& legend_strings, int ymin_flag)
{
 const int maxsymbols=6;
 const int maxcolors=8;
 int i,j,flag;
 double r,yymin;
 plot_curves levels(Nlevels);
 int colors[maxcolors]={blue,red,green4,black,magenta,orange,cyan,green};
 int symbols[maxsymbols]={circle,square,triangle_up,triangle_down,triangle_left,triangle_right};

 ymax=0.0;
 yymin=1e4;
 for (i=1;i<=Nlevels;i++){
    flag=levels.setup_from_file(i,datafiles(i),channel_match(i),Rplot_type,
                                no_line,white,line_thickness,symbols[(i-1)%maxsymbols],colors[(i-1)%maxcolors],
                                solid_fill,symbol_size,legend_strings(i));
//    if (flag){
//       cerr << " problem setting up plots ... "<<i<<endl;
//       return;}
    if (!flag){
       for (j=1;j<=levels(i).Npoints;j++){
          r=levels(i).xvals(j);
          if ((r>=xmin)&&(r<=xmax)){
             r=levels(i).yvals(j);
             if (r>ymax) ymax=r;
             if (r<yymin) yymin=r;
             }
          }
        }
    }
 if (ymin_flag) ymin=ymax-1.3*(ymax-yymin);
 else ymin=0.0;
 ymax=1.3*(ymax-ymin)+ymin;
 
 r=0.333*(ymax-ymin);
 if (r<=0.15){      y_major_tick_spacing=0.1;  y_minor_tick_num=4;}
 else if (r<=0.35){ y_major_tick_spacing=0.25; y_minor_tick_num=4;}
 else if (r<=0.65){ y_major_tick_spacing=0.5;  y_minor_tick_num=4;}
 else if (r<=1.35){ y_major_tick_spacing=1.0;  y_minor_tick_num=4;}
 else {             y_major_tick_spacing=2.0;  y_minor_tick_num=4;}

 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,
              legend_xpos,legend_ypos);

#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="/usr/local/bin/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 Nlevels,
                   const Svector& datafiles, const Svector& channel_match,
                   const Svector& legend_strings)
{
 pump_out_plot(outfile,Nlevels,datafiles,channel_match,legend_strings,0);
}

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

void html_header(const charString filename, const charString& up_filename,
                 const charString title1, const charString title2)
{
 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></td>"<<endl;
 html_out << "<center>"<<endl;
 html_out << "<h2><font COLOR=\"#990000\">"<<title1<<"</font></h2>"<<endl;
 html_out << "<h2><font COLOR=\"#000099\">"<<title2<<"</font></h2>"<<endl;
 html_out << "<br><br><br>"<<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 Nlevels,
                          const Svector& datafiles, const Svector& channel_match,
                          const Svector& legend_strings, int string_level, int dim)
{
 const double pi=3.141592653589793238;
 const int maxsymbols=6;
 const int maxcolors=8;
 int i,j,flag;
 double r,dr,y,dy;
 plot_curves levels(Nlevels);
 int colors[maxcolors]={blue,red,green4,black,magenta,orange,cyan,green};
 int symbols[maxsymbols]={circle,square,triangle_up,triangle_down,triangle_left,triangle_right};
 
 ymax=-1e4;
 ymin=1e4;
 for (i=1;i<=Nlevels;i++){
    flag=levels.setup_from_file(i,datafiles(i),channel_match(i),1,
                                no_line,white,line_thickness,symbols[(i-1)%maxsymbols],colors[(i-1)%maxcolors],
                                solid_fill,symbol_size,legend_strings(i));
//    if (flag){
//       cerr << " problem setting up plots";
//       return;}
    if (!flag){
       for (j=1;j<=levels(i).Npoints;j++){
       r=levels(i).xvals(j);
       dr=levels(i).dxvals(j);
       y=levels(i).yvals(j);
       dy=levels(i).dyvals(j);
       levels(i).yvals(j)=y/string_level-1.0;
       levels(i).dyvals(j)=dy/string_level;
       }
    for (j=1;j<=levels(i).Npoints;j++){
       r=levels(i).xvals(j);
       if ((r>=xmin)&&(r<=xmax)){
          r=levels(i).yvals(j);
          if (r>ymax) ymax=r;
          if (r<ymin) ymin=r;
          }
       }}
    }
 ymin=ymax-1.3*(ymax-ymin);
 ymax=1.3*(ymax-ymin)+ymin;
 
 r=0.333*(ymax-ymin);
 if (r<=0.15){      y_major_tick_spacing=0.1;  y_minor_tick_num=4;}
 else if (r<=0.35){ y_major_tick_spacing=0.25; y_minor_tick_num=4;}
 else if (r<=0.65){ y_major_tick_spacing=0.5;  y_minor_tick_num=4;}
 else if (r<=1.35){ y_major_tick_spacing=1.0;  y_minor_tick_num=4;}
 else {             y_major_tick_spacing=2.0;  y_minor_tick_num=4;}

 double xx[2],yy[2];
 xx[0]=xmin; xx[1]=xmax;
 yy[0]=yy[1]=0.0;
 comparison_curves.setup(1,2,xx,yy,solid_line,black,2*line_thickness,
                         no_symbol,white,none,0.8);
 setup_one_loop_nambu_goto(2.5,xmax,string_level,dim);
 setup_full_nambu_goto(1.0,xmax,string_level,dim);
 
 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,comparison_curves,
              label_strings,legend_xpos,legend_ypos);

#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
}

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