#include "analysis.h"

int max_level=5;

 //  This program makes the html files for a given simulation.
 //  It takes as input a directory name.  This directory is where
 //  the upper html file will be placed.  In this directory should
 //  be a file named "run_parameters.dat" from which all information
 //  about the run will be extracted.  In this directory should be
 //  several subdirectories: one for each symmetry channel. In each
 //  of these subdirectories should be a file which contains a
 //  summary of the fit results.

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


   //  Read lines from stream "in" starting a beginning of file until a line
   //  is found containing the string "match".  Returns value 1 if successful,
   //  0 otherwise. Returns the line in "str" if substring_flag=0, otherwise
   //  it returns in "str" that part of the string which comes after the matched
   //  pattern.
   
int read_lines_until(ifstream& in, charString& str, const charString& match,
                     int substring_flag)
{
 const int maxline=132;
 static char ff[maxline];
 int j,flag;
 
 in.clear(); // in case eof reached
 in.seekg(0,ios::beg); // go back to beginning of file
 flag=0;
 while (!flag){
    j=0;
    in.get(ff[j]); 
    while ((ff[j] != '\n')&&(j<(maxline-1))&&(!in.eof()))
       in.get(ff[++j]);
    ff[j]='\0';
    str=ff;
    flag=str.locate(match); 
    if (flag>=0){
       if (substring_flag) str=str.right_substring(match);
       return 1;
       }
    else flag=0;
    if (in.eof()){
       in.clear();
       return 0;}
    }
 return 0;
}

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

int error_quit(ifstream& in, ofstream& out, const charString& message)
{
 in.close();
 out.close();
 cerr << " ERROR RETURN:  "<<message<<" not found"<<endl;
 return 1;
}

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

charString trim1(charString& str)
{
 str=str.substring_replace(".0000",".0");
 str=str.substring_replace("000","");
 str=str.substring_replace(" ","");
 return str;
}

charString trim(charString& str)
{
 str=str.substring_replace(" ","");
 return str;
}

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

double calc_Rmagnitude(const charString& Rstring, int Nspace_dim)
{
 int x,y,z;
 double r;
 charString str=Rstring.substring('(',',');
 if (string_to_int(str,x)){
    cerr << "conversion to integer not possible in calc_Rmagnitude"<<endl;
    exit(1);}
 if (Nspace_dim==3){
    str=Rstring.substring(',',',');
    if (string_to_int(str,y)){
       cerr << "conversion to integer not possible in calc_Rmagnitude"<<endl;
       exit(1);}
    }
 else y=0;
 str=Rstring.substring(',',')');
 if (Nspace_dim==3) str=str.right_substring(',');
 if (string_to_int(str,z)){
    cerr << "conversion to integer not possible in calc_Rmagnitude"<<endl;
    exit(1);}
 if ((y==0)&&(z==0)) r=(double)x;
 else r=sqrt((double)(x*x+y*y+z*z));
 return r;
}

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

   // This routine looks in the file "resultsfilename" and proceeds until
   // the string "Rstring" is found.  From the **next** line in this file,
   // the answer is extracted and returned in "answer".  From the **next**
   // line after this, the Q fit value is extracted and returned in "Qval".
   // It also extracts the answer with error from the "Final fit results" lines.
   // If the "Final fit" line contains the string "OMIT", then "-" is returned in
   // answer with Qval=-1 and st_dev is returned with a negative value.  
   // If the string "Rstring" is not found, returns 1; successful returns 0.

int get_results(const charString& resultsfilename, const charString& Rstring,
                charString& answer_string, double& Qval, int orientation,
                double& answer, double& st_dev, int gap_flag)
{
 const int maxline=132;
 static char ff[maxline];
 int j,flag;
 charString str;
 charString match;

#ifdef DEBUG
cout << "resultsfilename = "<<resultsfilename<<endl;
cout << "Rstring = "<<Rstring<<endl;
cout << "orientation = "<<orientation<<endl;
cout << "gap_flag = "<<gap_flag<<endl;
#endif

 if (orientation==0){
    match="(R="+Rstring.substring('(',',')+")";
    str=resultsfilename+"_res.dat";}
 else if (orientation==1){
    match="[R="+Rstring+"]";
    str=resultsfilename+"_PD_res.dat";}
 else if (orientation==2){
    match="[R="+Rstring+"]";
    str=resultsfilename+"_CD_res.dat";}

#ifdef DEBUG
cout << " searching file "<<str<<endl;
#endif

 char *ifname=convert_to_charstar(str);
 ifstream in(ifname);
 delete [] ifname;
 if (!in){
    //cerr << " *** error opening file "<<str<<endl;
    return 1;}
    
 flag=0;
 while (!flag){
    j=0;
    in.get(ff[j]); 
    while ((ff[j] != '\n')&&(j<(maxline-1))&&(!in.eof()))
       in.get(ff[++j]);
    ff[j]='\0';
    str=ff; str=trim(str);
    flag=str.locate(match);
    if (flag<0) flag=0;
    else flag=1;
    if (in.eof()){
       in.clear();
       //cerr << "could not find "<<match<<endl;
       in.close();
       return 1;}
    }

     // get next line
     
 j=0;
 in.get(ff[j]); 
 while ((ff[j] != '\n')&&(j<(maxline-1))&&(!in.eof()))
    in.get(ff[++j]);
 ff[j]='\0';
 str=ff;
 if (!gap_flag) flag=str.locate("E=");
 else flag=str.locate("Delta");
 if (flag<0){
    cerr << "problem finding E or Delta value"<<endl;
    in.close();
    return 1;}
 if (!gap_flag) answer_string=str.right_substring("E=");
 else answer_string=str.right_substring(")=");

     // get next line
     
 j=0;
 in.get(ff[j]); 
 while ((ff[j] != '\n')&&(j<(maxline-1))&&(!in.eof()))
    in.get(ff[++j]);
 ff[j]='\0';
 str=ff; 
 flag=str.locate("Q=");
 if (flag<0){
    cerr << "problem finding Q value"<<endl;
    in.close();
    return 1;}
 if (string_to_double(str.right_substring("Q="),Qval)){
    cerr << "conversion to double not possible"<<endl;
    in.close(); return 1;}

// cout << "answer string = "<<answer_string<<"   Q = "<<Qval<<endl;

   // continue on down the file to find the corresponding line
   // with the result and error on a single line
   
 flag=0;
 while (!flag){
    j=0;
    in.get(ff[j]); 
    while ((ff[j] != '\n')&&(j<(maxline-1))&&(!in.eof()))
       in.get(ff[++j]);
    ff[j]='\0';
    str=ff; //cout << "str = "<<str<<endl;
    if (in.eof()){
       in.clear();
       //cerr << "could not find "<<Rstring<<endl;
       in.close();
       return 1;}
    flag=str.locate(match);// cout << "flag = "<<flag<<endl;
    if (flag<0) flag=0;
    else if (str.locate("Final fit result:")>=0) flag=1;
    else{
       cerr << " problem with file"<<resultsfilename<<endl;
       in.close(); return 1;}
    }
    
 flag=str.locate("OMIT");
 if (flag>=0){
    answer_string="------------";
    Qval=-1.0;
    answer=0.0;
    st_dev=-1.0;
    }
 else{
    str=str.right_substring("Final fit result:");
    // cout << "str = "<<str<<endl;
    j=0; while (str[j]==' ') ++j; 
    str=str.right_substring(j);  // remove leading blanks
    // cout << "str = "<<str<<endl;
    // cout << " string to be converted to double is "<<str.left_substring(' ')<<endl;
    if (string_to_double(str.left_substring(' '),answer)){
        cerr << " convert of answer to double not possible"<<endl;
        in.close();return 1;}
    if (string_to_double(str.right_substring(' '),st_dev)){
        cerr << " convert of st_dev to double not possible"<<endl;
        in.close();return 1;}
    }
 in.close();
#ifdef DEBUG
 cout << "answer string = "<<answer_string<<"   Q = "<<Qval<<" answer = "<<answer<<" st_dev = "<<st_dev<<endl;
 read_string(" continue? ",str);
#endif
 return 0;
}

// *******************************************************************
   
int get_orientations(ifstream& in, Svector& Rstrings, int& Norientations,
                     Svector& Rfilenames, int Nspace_dim, ivector& orientations)
{
 const int maxline=132;
 static char ff[maxline];
 int j,flag;
 charString str;
 
 in.clear(); // in case eof reached
 in.seekg(0,ios::beg); // go back to beginning of file
 flag=Norientations=0;
 while (!flag){
    j=0;
    in.get(ff[j]); 
    while ((ff[j] != '\n')&&(j<(maxline-1))&&(!in.eof()))
       in.get(ff[++j]);
    ff[j]='\0';
    str=ff; 
    flag=str.locate("Orientations studied");
    if (flag<0) flag=0;
    else flag=1;
    if (in.eof()){
       in.clear();
       return 0;}
    }
 flag=0;
 while (!flag){
    j=0;
    in.get(ff[j]); 
    while ((ff[j] != '\n')&&(j<(maxline-1))&&(!in.eof()))
       in.get(ff[++j]);
    ff[j]='\0';
    str=ff;  
    str=trim(str); 
    flag=str.locate("Maximum");
    if (flag<0){
       Norientations++;
       flag=0;}
    else flag=1;
    if (in.eof()){
       in.clear();
       return 0;}
    }

 Rfilenames.reset(Norientations);
 Rstrings.reset(Norientations);
 orientations.reset(Norientations);

 in.clear();
 in.seekg(0,ios::beg); // go back to beginning of file
 flag=Norientations=0;
 while (!flag){
    j=0;
    in.get(ff[j]); 
    while ((ff[j] != '\n')&&(j<(maxline-1))&&(!in.eof()))
       in.get(ff[++j]);
    ff[j]='\0';
    str=ff; 
    flag=str.locate("Orientations studied");
    if (flag<0) flag=0;
    else flag=1;
    if (in.eof()){
       in.clear();
       return 0;}
    }
 flag=0;   
 while (!flag){
    j=0;
    in.get(ff[j]); 
    while ((ff[j] != '\n')&&(j<(maxline-1))&&(!in.eof()))
       in.get(ff[++j]);
    ff[j]='\0';
    str=ff; 
    str=trim(str); 
    flag=str.locate("Maximum");
    if (flag<0){
       Norientations++;
    
       if (Nspace_dim==2){
          if (str.locate(",0)")>=0){
             Rstrings(Norientations)=str;
             Rfilenames(Norientations)="R"+str.substring('(',',');
             orientations(Norientations)=0;
             }
          else{
             Rstrings(Norientations)=str;
             Rfilenames(Norientations)="PD_R"+str.substring('(',',');
             orientations(Norientations)=1;
             }
          }
       else if (Nspace_dim==3){
          if (str.locate(",0,0)")>=0){
             Rstrings(Norientations)=str;
             Rfilenames(Norientations)="R"+str.substring('(',',');
             orientations(Norientations)=0;
             }
          else if (str.locate(",0)")>=0){
             Rstrings(Norientations)=str;
             Rfilenames(Norientations)="PD_R"+str.substring('(',',');
             orientations(Norientations)=1;
             }
          else{
             Rstrings(Norientations)=str;
             Rfilenames(Norientations)="CD_R"+str.substring('(',',');
             orientations(Norientations)=2;
             }
          }
    
        
        flag=0;}
    else flag=1;
    if (in.eof()){
       in.clear();
       return 0;}
    }
 return 0;
}

// *******************************************************************
//
//                DRIVER ROUTINE
//
// *******************************************************************

int write_pages(const charString& directory, const charString& mainoutfile,
                const charString& run_name, int Nspace_dim, int Nirreps,
                const Svector& labels, double glueball_mass)
{
 int i,j,k,l,row,col,Ncol,flag;
 const int maxcol=8;
 double Qval,ans,st_dev;
 charString answer,str,stub;

   // ****************************
   //
   // first, extract run parameters from the run_parameters.dat file,
   // and get scale parameter information
   //
   // ****************************

 char *ifname=convert_to_charstar(directory+"/run_parameters.dat");
 ifstream in(ifname);
 delete [] ifname;
 if (!in){
    cerr << " *** error opening file run_parameters.dat"<<endl;
    return 1;}

 charString sigma,sigma_inv,zxi;
 
 char *sfname=convert_to_charstar(directory+"/scale/scale_setting.html");
 ifstream scalein(sfname);
 delete [] sfname;
 if (!scalein){
    sigma="???"; sigma_inv="???";}
 else{
    if (read_lines_until(scalein,str,"a_s a_t sigma",1)==0) sigma="???";
    else sigma=" 0."+str.substring('.',')')+")";
    if (read_lines_until(scalein,str,"1/(a_s a_t sigma)",1)==0) sigma_inv="???";
    else{ 
       charString tmp=str.left_substring('.'); 
       while (tmp.right_relax_substring('>')!=tmp) tmp=tmp.right_relax_substring('>');
       sigma_inv=tmp+"."+str.substring('.',')')+")";}
    }
 scalein.close();
    
 char *zfname=convert_to_charstar(directory+"/scale/zxi.html");
 ifstream zin(zfname);
 delete [] zfname;
 if (!zin){
    zxi="???";}
 else{
    if (read_lines_until(zin,str,"Z_xi",1)==0) zxi="???";
    else{ 
       charString tmp=str.left_substring('.'); 
       while (tmp.right_relax_substring('>')!=tmp) tmp=tmp.right_relax_substring('>');
       zxi=tmp+"."+str.substring('.',')')+")";}
    }
 zin.close();

 char *fname=convert_to_charstar(mainoutfile);
 ofstream out(fname);
 delete [] fname;
 if (!out){
    cerr << " *** error opening file "<<mainoutfile<<endl;
    return 1;}

 out << "<HTML>" << endl;
 out << "<HEAD>" << endl;
 out << "<TITLE> Run " << run_name <<" </TITLE>"<<endl;
 out << "<!-- code below removes annoying underlines from all links -->"<<endl;
 out << "<style>"<<endl;
 out << "<!--"<<endl;
 out << "a{text-decoration:none}"<<endl;
 out << "//-->"<<endl;
 out << "</style>"<<endl;
 out << "</HEAD>"<<endl<<endl;

 if (read_lines_until(in,str,"static",0)==0) return error_quit(in,out,"static");
 str=str.substring(' ',' ');

 out << "<BODY BGCOLOR=\"#FFFFCC\" LINK=\"#000099\" ALINK=\"#000099\" VLINK=\"#000099\">"<<endl;
 out << "<a href=\"../raw_results.html\">" << endl;
 out << "<img src=\"../../images/back.gif\" align=middle border=0></a>"<<endl;
 out << "<center>"<<endl;
 out << "<H2><FONT COLOR=\"#990000\">"<<endl;
 out << "Static quark potential and its excitations<br>"<<endl;
 out << "in "<<int_to_string(Nspace_dim+1)<<"-dimensional "<<str<<" gauge theory</FONT></H2>"<<endl;
 out << ""<<endl;
 out << "<H2><FONT COLOR=\"#FF6666\">Run "<<run_name<<"</FONT></H2>"<<endl;
 out << "<br>"<<endl;
 out << ""<<endl;
 out << "<H3><font color=green> RUN PARAMETERS</font></H3>"<<endl<<endl;

 out << "<table><td valign=top>"<<endl;
 out << "<table border=8 cellpadding=5 bgcolor=\"#99FFCC\" bordercolor=blue>"<<endl;
 out << "<tr align=center><td bgcolor=\"#CCFFFF\" colspan=2><b> Improved anisotropic action </b></td></tr>"<<endl;
 if (read_lines_until(in,str,"beta =",1)==0) return error_quit(in,out,"beta");
 out << "<tr align=center  ><td bgcolor=\"#66FFFF\"><b> beta                       </b></td><td><b>"<<trim1(str)<<" </b></td></tr>"<<endl;
 if (read_lines_until(in,str,"anisotropy a_s / a_t =",1)==0) return error_quit(in,out,"anisotropy");
 out << "<tr align=center  ><td bgcolor=\"#66FFFF\"><b> input anisotropy a_s / a_t </b></td><td><b>"<<trim1(str)<<" </b></td></tr>"<<endl;
 if (read_lines_until(in,str,"input u_t =",1)==0) return error_quit(in,out,"u_t");
 out << "<tr align=center  ><td bgcolor=\"#66FFFF\"><b> input u_t                  </b></td><td><b>"<<trim(str)<<" </b></td></tr>"<<endl;
 if (read_lines_until(in,str,"input u_s =",1)==0) return error_quit(in,out,"u_s");
 out << "<tr align=center  ><td bgcolor=\"#66FFFF\"><b> input u_s                  </b></td><td><b>"<<trim(str)<<" </b></td></tr>"<<endl;
 out << "</table>"<<endl;
 out << "</td><td width=50></td><td valign=top>"<<endl;
 out << "<table border=8 cellpadding=5 bgcolor=\"#99FFCC\" bordercolor=blue>"<<endl;
 out << "<tr align=center><td bgcolor=\"#CCFFFF\" colspan=2><b> Individual term couplings: </b></td></tr>"<<endl;
 if (read_lines_until(in,str,"spatial plaquette",1)==0) return error_quit(in,out,"sp plaq");
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>      spatial plaquette  </b></td><td><b>"<<trim(str)<<"</b></td></tr>"<<endl;
 if (read_lines_until(in,str,"temporal plaquette",1)==0) return error_quit(in,out,"tp plaq");
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>     temporal plaquette  </b></td><td><b>"<<trim(str)<<"</b></td></tr>"<<endl;
 if (read_lines_until(in,str,"spatial 2x1 rectangle",1)==0) return error_quit(in,out,"sp 2x1");
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>  spatial 2x1 rectangle  </b></td><td><b>"<<trim(str)<<"</b></td></tr>"<<endl;
 if (read_lines_until(in,str,"2s x 1t rectangle",1)==0) return error_quit(in,out,"2s x 1t");
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>      2s x 1t rectangle  </b></td><td><b>"<<trim(str)<<"</b></td></tr>"<<endl;
 if (read_lines_until(in,str,"2t x 1s rectangle",1)==0) return error_quit(in,out,"2t x 1s");
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>      2t x 1s rectangle  </b></td><td><b>"<<trim(str)<<"</b></td></tr>"<<endl;
 out << "</table>"<<endl;
 out << "</td></table>"<<endl;
 out << "<br>"<<endl<<endl;

 out << "<table border=8 cellpadding=5 bgcolor=\"#99FFCC\" bordercolor=blue>"<<endl;
 if (read_lines_until(in,str,"Lattice size",1)==0) return error_quit(in,out,"lattice size");
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>    lattice size  </b></td><td><b>"<<str<<"</b></td></tr>"<<endl;
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>    scale setting a_s a_t sigma  </b></td><td><a href=\"scale/scale_setting.html\"> <b> "<<sigma<<" </b></a></td></tr>"<<endl;
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>    inverse scale 1/(a_s a_t sigma)  </b></td><td><a href=\"scale/scale_setting.html\"> <b> "<<sigma_inv<<" </b></a></td></tr>"<<endl;
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>    anisotropy renormalization Z_xi </b></td><td><a href=\"scale/zxi.html\"><b> "<<zxi<<" </b></a></td></tr>"<<endl;
 out << "</table>"<<endl;
 out << "<br>"<<endl;
 out << ""<<endl;
 out << "<table><td valign=top>"<<endl;
 out << "<table border=8 cellpadding=5 bgcolor=\"#99FFCC\" bordercolor=blue>"<<endl;
 out << "<tr align=center><td bgcolor=\"#CCFFFF\" colspan=2><b> Configuration updating parameters: </b></td></tr>"<<endl;
 if (read_lines_until(in,str,"number of bins",1)==0) return error_quit(in,out,"number of bins");
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>                                number of bins  </b></td><td><b>"<<str<<"</b></td></tr>"<<endl;
 if (read_lines_until(in,str,"number of measurements per bin",1)==0) return error_quit(in,out,"meas per bin");
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>                number of measurements per bin  </b></td><td><b>"<<str<<"</b></td></tr>"<<endl;
 if (read_lines_until(in,str,"number of updates between measurements",1)==0) return error_quit(in,out,"updates per meas");
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>        number of updates between measurements  </b></td><td><b>"<<str<<"</b></td></tr>"<<endl;
 if (read_lines_until(in,str,"number of Cabibbo-Marinari sweeps per update",1)==0) return error_quit(in,out,"cm sweeps");
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>  number of Cabibbo-Marinari sweeps per update  </b></td><td><b>"<<str<<"</b></td></tr>"<<endl;

 if (read_lines_until(in,str,"number of over-relaxation sweeps (all links) per update",1)){
    out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>   number of over-relaxation sweeps (all links) per update  </b></td><td><b>"<<str<<"</b></td></tr>"<<endl;
    if (read_lines_until(in,str,"number of over-relaxation sweeps (spatial links) per update",1)==0) return error_quit(in,out,"OR spatial");
    out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>   number of over-relaxation sweeps (spatial links) per update  </b></td><td><b>"<<str<<"</b></td></tr>"<<endl;
    }
 else if (read_lines_until(in,str,"number of over-relaxation sweeps per update",1)){
    out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>   number of over-relaxation sweeps per update  </b></td><td><b>"<<str<<"</b></td></tr>"<<endl;
    }
 else{ cout << read_lines_until(in,str,"number of over-relaxation sweeps per update",1)<<endl;
    return error_quit(in,out,"OR all");}
 
 out << "</table>"<<endl;
 out << "</td><td width=50></td><td valign=top>"<<endl;
 out << "<table border=8 cellpadding=5 bgcolor=\"#99FFCC\" bordercolor=blue>"<<endl;
 
 if (read_lines_until(in,str,"Blocking parameters",1)!=0){
    out << "<tr align=center><td bgcolor=\"#CCFFFF\" colspan=2><b> Blocking parameters: </b></td></tr>"<<endl;
    if (read_lines_until(in,str,"blocking alpha =",1)==0) return error_quit(in,out,"blocking alpha");
    out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>  blocking alpha  </b></td><td><b>"<<str<<"</b></td></tr>"<<endl;
    if (read_lines_until(in,str,"blocking levels =",1)==0) return error_quit(in,out,"blocking levels");
    out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>  blocking levels      </b></td><td><b>"<<str<<"</b></td></tr>"<<endl;
    }

 out << "<tr align=center><td bgcolor=\"#CCFFFF\" colspan=2><b> Link smearing parameters: </b></td></tr>"<<endl;
 if (read_lines_until(in,str,"      alpha =",1)==0) return error_quit(in,out,"Link smearing alpha");
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>  alpha (staple weight)  </b></td><td><b>"<<str<<"</b></td></tr>"<<endl;
 if (read_lines_until(in,str,"fuzzing levels =",1)==0) return error_quit(in,out,"link smearing levels");
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>  number of levels  </b></td><td><b>"<<str<<"</b></td></tr>"<<endl;

 out << "<tr align=center><td bgcolor=\"#CCFFFF\" colspan=2><b> Thermal averaging of temporal links: </b></td></tr>"<<endl;
 if (read_lines_until(in,str,"number of Cabibbo-Marinari hits =",1)!=0)
    out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>  number of Cabibbo-Marinari hits  </b></td><td><b>"<<str<<"</b></td></tr>"<<endl;
 else
    out << "<tr align=center><td bgcolor=\"#66FFFF\" colspan=2>  done analytically  </td></tr>"<<endl;

 out << "</table>"<<endl;
 out << "</td></table>"<<endl;
 out << "<br>"<<endl<<endl;

 out << "<table><td valign=top>"<<endl;
 out << "<table border=8 cellpadding=5 bgcolor=\"#99FFCC\" bordercolor=blue>"<<endl;
 out << "<tr align=center><td bgcolor=\"#CCFFFF\" colspan=2><b> Correlator information:  </b></td></tr>"<<endl;
 if (read_lines_until(in,str,"Maximum time separation measured",1)==0) return error_quit(in,out,"max time measured");
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>        Maximum time separation measured </b></td><td><b>"<<str<<"</b></td></tr>"<<endl;
 if (read_lines_until(in,str,"perpendicular spatial-translation increment delta-S",1)==0) return error_quit(in,out,"delta-S");
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>           Perpendicular spatial-translation increment </b></td><td><b>"<<str<<"</b></td></tr>"<<endl;
 if (read_lines_until(in,str,"parallel spatial-translation increment delta-S",1)==0) return error_quit(in,out,"delta-S");
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>           Parallel spatial-translation increment </b></td><td><b>"<<str<<"</b></td></tr>"<<endl;
 if (read_lines_until(in,str,"temporal-translation increment delta-T",1)==0) return error_quit(in,out,"delta-T");
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>          Temporal-translation increment </b></td><td><b>"<<str<<"</b></td></tr>"<<endl;
 if (read_lines_until(in,str,"Number of optimized operators retained",1)==0) return error_quit(in,out,"opt opers");
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>  Number of optimized operators retained </b></td><td><b>"<<str<<"</b></td></tr>"<<endl;
 if (read_lines_until(in,str,"Optimized on time slices",1)==0) return error_quit(in,out,"time slice opt");
 out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>                Optimized on time slices </b></td><td><b>"<<str<<"</b></td></tr>"<<endl;
 out << "</table>"<<endl;
 out << "</td><td width=50></td><td valign=top>"<<endl;
 out << "<table border=8 cellpadding=5 bgcolor=\"#99FFCC\" bordercolor=blue>"<<endl;
 out << "<tr align=center><td bgcolor=\"#CCFFFF\" colspan=2><b> Number of operators before fuzzing: </b></td></tr>"<<endl;
 
 for (i=1;i<=Nirreps;i++){
    if (read_lines_until(in,str,"          "+labels(i),1)==0) return error_quit(in,out,labels(i));
    out << "<tr align=center><td bgcolor=\"#66FFFF\"><b>"<<labels(i)<<"</b></td><td><b>"<<str<<"</b></td></tr>"<<endl;
    }
    
 out << "</table>"<<endl;
 out << "</td></table>"<<endl;
 out << "<br><br><br>"<<endl<<endl;

 
      // ***********************************
      //
      //   Next, extract all of the orientations from the run_parameters.dat file.
      //
      // ***********************************
      
 int Norientations;
 Svector Rfilenames,Rstrings;
 ivector orientations;
 get_orientations(in,Rstrings,Norientations,Rfilenames,Nspace_dim,orientations);
 cout << "Norientations="<<Norientations<<endl;
 for (i=1;i<=Norientations;i++)
    cout << Rstrings(i)<<"  file ending = "<<Rfilenames(i)
         <<"  orientation = "<<orientations(i)<<endl;


    // ******************************************
    //
    //   Now make the table of results  
    //     determine which symmetry channels, levels to include
    //     determine which quark-antiquark orientations to include
    //
    // ******************************************
    
 int Nlevels;
 Svector datafilenames,levelnames,levelstubs;
 ifstream indata;
 dmatrix gapans_mean,gapans_stdev,gapans_Qval;
 dmatrix ans_mean,ans_stdev,ans_Qval;
 int Ngaplevels;
 Svector gapdatafilenames,gaplevelnames,gaplevelstubs;

 Nlevels=0;
 for (i=1;i<=Nirreps;i++){
    for (j=0;j<=max_level;j++){
       for (l=1;l<=Norientations;l++){
          str=directory+"/"+labels(i)+"/"+labels(i);
          for (k=1;k<=j;k++) str=str+"p";
          if (orientations(l)==1) str=str+"_PD";
          else if (orientations(l)==2) str=str+"_CD";
          str=str+"_res.dat"; //cout << " looking for results in file "<<str<<endl;
          indata.clear();
          indata.open(str.convert_to_charstar()); 
          if (indata){
             Nlevels++; cout << "Nlevels found so far = "<<Nlevels<<endl;
             indata.close();
             l=Norientations+1;
             break;}
          indata.close();
          }
       }}
       
 if (Nlevels>0){
 datafilenames.reset(Nlevels);
 levelnames.reset(Nlevels);
 levelstubs.reset(Nlevels);
 Nlevels=0;
 for (i=1;i<=Nirreps;i++){
    for (j=0;j<=max_level;j++){
       for (l=1;l<=Norientations;l++){
          stub=directory+"/"+labels(i)+"/"+labels(i);
          for (k=1;k<=j;k++) stub=stub+"p";
          str=stub;
          if (orientations(l)==1) str=str+"_PD";
          else if (orientations(l)==2) str=str+"_CD";
          str=str+"_res.dat"; // cout << " looking for results in file "<<str<<endl;
          indata.clear();
          indata.open(str.convert_to_charstar());
          if (indata){
             Nlevels++;
             datafilenames(Nlevels)=stub;
             levelnames(Nlevels)=labels(i);
             for (k=1;k<=j;k++) levelnames(Nlevels)=levelnames(Nlevels)+"'";
             levelstubs(Nlevels)=labels(i)+"/"+labels(i);
             for (k=1;k<=j;k++) levelstubs(Nlevels)=levelstubs(Nlevels)+"p";
             l=Norientations+1;
             indata.close();
             break;
             }
          indata.close();
          }
       }}
 cout << "Nlevels="<<Nlevels<<endl;
 for (i=1;i<=Nlevels;i++)
    cout << "  "<<levelnames(i)<<"   stub = "<<levelstubs(i)<<"   data file = "<<datafilenames(i)<<endl;

 out << "<H3> <FONT COLOR=\"#990000\"> Remarks concerning effective mass plots:</FONT></H3>"<<endl<<endl;

 out << "<table><td width=400 align=center><b>"<<endl;
 out << "Results below are obtained by fitting individual, but optimized, correlators"<<endl;
 out << "to a sum of two exponentials.  The fit results are indicated in the"<<endl;
 out << "effective mass plots by horizontal lines.  In the first set of raw results below,"<<endl;
 out << "the solid circles shown with error bars in the plots are standard effective mass"<<endl;
 out << "points associated with the optimized correlators.  The open symbols correspond to"<<endl;
 out << "<font color=red>principal effective mass </font>values as defined by "<<endl;
 out << "<font color=blue> Eq. (5.18)</font>"<<endl;
 out << "in <br><font color=blue>"<<endl;
 out << "  Luscher and Wolff, Nucl.Phys.B339, 222 (1990)</font> <br>"<<endl;
 out << "such that the lowest <i>N</i> effective"<<endl;
 out << "masses tend to the lowest <i>N</i> physical masses as the temporal separation"<<endl;
 out << "becomes large.  Note that the correlation functions themselves are used in"<<endl;
 out << "the fitting, not the effective mass points which are included for illustrative"<<endl;
 out << "purposes and as checks of the fitting results."<<endl;
 out << "In the second set of raw gap results below, all points with error bars correspond to"<<endl;
 out << "principal effective mass points.<br><br>"<<endl;
 out << "Entries below shown in <font color=red>red</font> indicate energies lying higher than one glueball mass"<<endl;
 out << "above the ground state; entries shown in <font color=amber>amber</font> indicate fits with quality factor"<<endl;
 out << "<font color=blue>Q<0.1</font>.<br><br>"<<endl;
 out << "In this run, the glueball mass in inverse temporal lattice spacings is approximately"<<endl;
 out << "<font color=blue> "<<glueball_mass<<"</font>.</b></td></table><br><br><br>"<<endl;

 out << "<H3><font color=green> RAW RESULTS<br>(in terms of temporal lattice spacing)</font></H3>"<<endl;
 out << "<br>"<<endl;

   // store results in array for later manipulation and output
   
 ans_mean.reset(Nlevels,Norientations);
 ans_stdev.reset(Nlevels,Norientations);
 ans_Qval.reset(Nlevels,Norientations);

 k=0;
 while (k<Nlevels){
 
         // output new table
         // first, determine number of columns
         
    Ncol=min(maxcol,Nlevels-k);
    i=Nlevels-(k+Ncol);  // leftover number of columns
    if ((i==1)&&(Ncol>2)) Ncol--;    // reduce by one so last table has at least 2 columns
    
    out << "<table border=8 cellpadding=5 bgcolor=\"#99FFCC\" bordercolor=blue>"
        <<endl;
    out << "<tr align=center><td bgcolor=\"#66FFFF\"><b> R </b></td>"<<endl;
    for (col=1;col<=Ncol;col++) 
       out << "<td><a href=\"" << levelstubs(k+col)<<".html\"><b>"
           <<levelnames(k+col)<<"</b></a></td>"<<endl;
    out << "</tr>"<<endl;

    for (row=1;row<=Norientations;row++){
       out << "<tr align=left><td bgcolor=\"#66FFFF\"><a href=\"Rlinks/"
           <<Rfilenames(row)<<".html\"><b>"
           << Rstrings(row)<<" </b></a></td>"<<endl;
       for (col=1;col<=Ncol;col++){      
            flag=get_results(datafilenames(k+col),Rstrings(row),answer,Qval,
                             orientations(row),ans,st_dev,0);
            if (!flag){
               ans_mean(k+col,row)=ans;
               ans_stdev(k+col,row)=st_dev;
               ans_Qval(k+col,row)=Qval;
               if (ans>= (ans_mean(1,row)+glueball_mass)){
                  out <<"<td><a href=\""<<levelstubs(k+col)<<"_"<<Rfilenames(row)
                      <<".gif\"><font color=red><b>"
                      << answer <<"</b></font></a></td>"<<endl;
                  ans_stdev(k+col,row)=-1.0; // negative error means omit
                  }
               else if (Qval>=0.1)
                  out <<"<td><a href=\""<<levelstubs(k+col)<<"_"<<Rfilenames(row)
                      <<".gif\"><b>"
                      << answer <<"</b></a></td>"<<endl;
               else if (Qval>=0.0)
                  out <<"<td><a href=\""<<levelstubs(k+col)<<"_"<<Rfilenames(row)
                      <<".gif\"><font color=amber><b>"
                      << answer <<"</b></font></a></td>"<<endl;
               else
                  out <<"<td align=center><a href=\""<<levelstubs(k+col)
                      <<"_"<<Rfilenames(row)<<".gif\">"
                      << answer <<"</a></td>"<<endl;
               }
            else{
               out <<"<td align=center> -------------</td>"<<endl;
               ans_stdev(k+col,row)=-1.0; // signals that no such answer
               }
            }
       out << "</tr>"<<endl;
       }
    out <<"</table>"<<endl;
    out <<"<BR><BR><BR>"<<endl;

    k+=Ncol;
    }
 }

    //   Now make the table of results for the gaps....  
    //     determine which symmetry channels, levels to include
    //     determine which quark-antiquark orientations to include
 
 Ngaplevels=0;
 for (i=1;i<=Nirreps;i++){
    for (j=0;j<=max_level;j++){
       for (l=1;l<=Norientations;l++){
          str=directory+"/"+labels(i)+"/"+labels(i);
          for (k=1;k<=j;k++) str=str+"p";
          str=str+"_gap";
          if (orientations(l)==1) str=str+"_PD";
          else if (orientations(l)==2) str=str+"_CD";
          str=str+"_res.dat"; //cout << " trying "<<str<<endl;
          indata.clear();
          indata.open(str.convert_to_charstar());
          if (indata){
              Ngaplevels++;
              indata.close();
              l=Norientations+1;
              break;}
          indata.close();
          }
       }}
 if (Ngaplevels>0){
 gapdatafilenames.reset(Ngaplevels);
 gaplevelnames.reset(Ngaplevels);
 gaplevelstubs.reset(Ngaplevels);
 Ngaplevels=0;
 for (i=1;i<=Nirreps;i++){
    for (j=0;j<=max_level;j++){
       for (l=1;l<=Norientations;l++){
          stub=directory+"/"+labels(i)+"/"+labels(i);
          for (k=1;k<=j;k++) stub=stub+"p";
          stub=stub+"_gap";
          str=stub;
          if (orientations(l)==1) str=str+"_PD";
          else if (orientations(l)==2) str=str+"_CD";
          str=str+"_res.dat";
          indata.clear();
          indata.open(str.convert_to_charstar());
          if (indata){
             Ngaplevels++;
             gapdatafilenames(Ngaplevels)=stub;
             gaplevelnames(Ngaplevels)=labels(i);
             for (k=1;k<=j;k++) 
                gaplevelnames(Ngaplevels)=gaplevelnames(Ngaplevels)+"'";
             gaplevelstubs(Ngaplevels)=labels(i)+"/"+labels(i);
             for (k=1;k<=j;k++)
                gaplevelstubs(Ngaplevels)=gaplevelstubs(Ngaplevels)+"p";
             gaplevelstubs(Ngaplevels)=gaplevelstubs(Ngaplevels)+"_gap";
             indata.close();
             l=Norientations+1;
             break;
             }
          indata.close();
          }
       }}
 cout << "Ngaplevels="<<Ngaplevels<<endl;
 for (i=1;i<=Ngaplevels;i++)
    cout << "  "<<gaplevelnames(i)<<"   stub = "<<gaplevelstubs(i)
         <<"   data file = "<<gapdatafilenames(i)<<endl;

 out << endl <<"<H3><font color=green> RAW GAP (above "<< labels(1) <<") RESULTS<br>(in terms of temporal lattice spacing)</font></H3>"<<endl;
 out << "<br>"<<endl;

   // store results in array for later manipulation and output
   
 gapans_mean.reset(Ngaplevels,Norientations);
 gapans_stdev.reset(Ngaplevels,Norientations);
 gapans_Qval.reset(Ngaplevels,Norientations);

 k=0;
 while (k<Ngaplevels){
 
         // output new table
         // first, determine number of columns
         
    Ncol=min(maxcol,Ngaplevels-k);
    i=Ngaplevels-(k+Ncol);  // leftover number of columns
    if ((i==1)&&(Ncol>2)) Ncol--;    // reduce by one so last table has at least 2 columns
    
    out << "<table border=8 cellpadding=5 bgcolor=\"#99FFCC\" bordercolor=blue>"
        <<endl;
    out << "<tr align=center><td bgcolor=\"#66FFFF\"><b> R </b></td>"<<endl;
    for (col=1;col<=Ncol;col++) 
       out << "<td><a href=\"" << gaplevelstubs(k+col)<<".html\"><b>"
           <<gaplevelnames(k+col)<<"</b></a></td>"<<endl;
    out << "</tr>"<<endl;

    for (row=1;row<=Norientations;row++){
       out << "<tr align=left><td bgcolor=\"#66FFFF\"><a href=\""
           <<"Rlinks/gap_"<<Rfilenames(row)<<".html\"><b>"
           << Rstrings(row)<<" </b></a></td>"<<endl;
       for (col=1;col<=Ncol;col++){      
            flag=get_results(gapdatafilenames(k+col),Rstrings(row),answer,Qval,
                             orientations(row),ans,st_dev,1);
            if (!flag){
               gapans_mean(k+col,row)=ans;
               gapans_stdev(k+col,row)=st_dev;
               gapans_Qval(k+col,row)=Qval;
               if (ans>= glueball_mass){
                  out <<"<td><a href=\""<<gaplevelstubs(k+col)
                      <<"_"<<Rfilenames(row)<<".gif\"><font color=red><b>"
                      << answer <<"</b></font></a></td>"<<endl;
                  gapans_stdev(k+col,row)=-1.0;  // negative error means omit
                  }
               else if (Qval>=0.1)
                  out <<"<td><a href=\""<<gaplevelstubs(k+col)
                      <<"_"<<Rfilenames(row)<<".gif\"><b>"
                      << answer <<"</b></a></td>"<<endl;
               else if (Qval>=0.0)
                  out <<"<td><a href=\""<<gaplevelstubs(k+col)
                      <<"_"<<Rfilenames(row)<<".gif\"><font color=amber><b>"
                      << answer <<"</b></font></a></td>"<<endl;
               else
                  out <<"<td align=center><a href=\""<<gaplevelstubs(k+col)
                      <<"_"<<Rfilenames(row)<<".gif\">"
                      << answer <<"</a></td>"<<endl;
               }
            else{
               out <<"<td align=center> -------------</td>"<<endl;
               gapans_stdev(k+col,row)=-1.0; // signals that no such answer
               }
            }
       out << "</tr>"<<endl;
       }
    out <<"</table>"<<endl;
    out <<"<BR><BR><BR>"<<endl;

    k+=Ncol;
    }
 }
 
 out << endl;
 out << "<H3><font color=green> RAW RESULTS PLOTS </font></H3>"<<endl;
 out << "<table border=8 cellpadding=5 bgcolor=\"#99FFCC\" bordercolor=blue>"<<endl;
 out << "<tr align=left><td><a href=\"plots/Rplots.html\">"<<endl;
 out << "<font COLOR=\"#990000\"><b>Energies plotted against R</b></font></a></td></tr>"<<endl;
 out << "<tr align=left><td><a href=\"plots/gap_Rplots.html\">"<<endl;
 out << "<font COLOR=\"#990000\"><b>Energy gaps above "<<labels(1)<<" plotted against R</b></font></a></td></tr>"<<endl;
 out << "<tr align=left><td><a href=\"plots/gap_Rinv_plots.html\">"<<endl;
 out << "<font COLOR=\"#990000\"><b>Energy gaps above "<<labels(1)<<" plotted against 1/R</b></font></a></td></tr>"<<endl;
 out << "<tr align=left><td><a href=\"plots/string_plots.html\">"<<endl;
 out << "<font COLOR=\"#990000\"><b>Energy gaps above "<<labels(1)<<" compared to <i>N</i> Pi/R</b></font></a></td></tr>"<<endl;
 out << "</table>"<<endl;
 out << "<BR><BR><BR>"<<endl;
 out << endl;
 out << "<H3><font color=green> TEXTUAL DATA FILES</font></H3>"<<endl;
 out << "<table border=8 cellpadding=5 bgcolor=\"#99FFCC\" bordercolor=blue>"<<endl;
 out << "<tr align=left><td><a href=\"all_results.dat\">"<<endl;
 out << "<font COLOR=\"#990000\"><b>Energies</b></font></a></td></tr>"<<endl;
 out << "<tr align=left><td><a href=\"all_gap_results.dat\">"<<endl;
 out << "<font COLOR=\"#990000\"><b>Energy gaps above "<<labels(1)<<"</b></font></a></td></tr>"<<endl;
 out << "</table>"<<endl;
 out << "<BR><BR><BR>"<<endl;
 out << "</center>"<<endl;
 out << endl;
 out << "<a href=\"../raw_results.html\">"<<endl;
 out << "<img src=\"../../images/back.gif\" align=middle border=0></a>"<<endl;
 out << "<a href=\"#top\">"<<endl;
 out << "<img src=\"../../images/to_top.gif\" align=middle border=0></a>"<<endl;
 out << "<br><br>"<<endl;
 out << "</BODY>"<<endl;
 out << "</HTML>"<<endl;

 in.close();
 out.close();
 
 
  // ***********************
  //
  //  Now write the R**.html files
  //
  // ***********************
  
 for (k=1;k<=Norientations;k++){

    fname=convert_to_charstar("Rlinks/"+Rfilenames(k)+".html");
    out.clear();
    out.open(fname);
    delete [] fname;
    if (!out){
       cerr << " *** error opening file Rlinks/"<<Rfilenames(k)<<".html"<<endl;
       return 1;}

    out << "<HTML>"<<endl;
    out << "<HEAD>"<<endl;
    out << "<TITLE> Fit results</TITLE>"<<endl;
    out << "</HEAD>"<<endl;
    out << "<BODY BGCOLOR=\"#99FFCC\">"<<endl;
    out << "<a href=\"../"<<mainoutfile<<"\">"<<endl;
    out << "<img src=\"../../../images/back.gif\" align=middle border=0></a>"<<endl;
    out << "<center>"<<endl;
    out << "<H2><FONT COLOR=\"#990000\">Fit results for the energy of the</FONT><br>"<<endl;
    out << "<FONT COLOR=\"#990000\"> stationary states of the gauge field</FONT><br>"<<endl;
    out << "<FONT COLOR=\"#990000\">in the presence of a static source-antisource pair</FONT><br>"<<endl;
    out << "<FONT COLOR=\"#990000\"> for </FONT>"<<endl;
    out << "<FONT COLOR=black> R = "<<Rstrings(k)<<" </FONT>"<<endl;
    out << "<FONT COLOR=\"990000\"> for Run "<<run_name<<"</FONT><br><br>"<<endl;
    out << "</h2>"<<endl;
    out << "<p><br>"<<endl<<endl;

    for (i=1;i<=Nlevels;i++){
       if (ans_stdev(i,k)>0.0){
          out << "<IMG SRC=\"../"<< levelstubs(i)<<"_"<<Rfilenames(k) <<".gif\"><br><br><br><br>"<<endl;}}

    out << endl<<"</center>"<<endl;
    out << "<BR><BR><BR> "<<endl;
    out << "<a href=\"../"<<mainoutfile<<"\">"<<endl;
    out << "<img src=\"../../../images/back.gif\" align=middle border=0></a>"<<endl;
    out << "<a href=\"#top\">"<<endl;
    out << "<img src=\"../../../images/to_top.gif\" align=middle border=0></a>"<<endl;
    out << "</BODY>"<<endl;
    out << "</HTML>"<<endl;
    
    out.close();

    } 
    
 //  Now write the column level.html files

 for (k=1;k<=Nlevels;k++){

    fname=convert_to_charstar(levelstubs(k)+".html");
    out.clear();
    out.open(fname);
    delete [] fname;
    if (!out){
       cerr << " *** error opening file "<<levelstubs(k)<<".html"<<endl;
       return 1;}

    out << "<HTML>"<<endl;
    out << "<HEAD>"<<endl;
    out << "<TITLE> Fit results</TITLE>"<<endl;
    out << "</HEAD>"<<endl;
    out << "<BODY BGCOLOR=\"#99FFCC\">"<<endl;
    out << "<a href=\"../"<<mainoutfile<<"\">"<<endl;
    out << "<img src=\"../../../images/back.gif\" align=middle border=0></a>"<<endl;
    out << "<center>"<<endl;
    out << "<H2><FONT COLOR=\"#990000\">Fit results for the energy of the</FONT><br>"<<endl;
    out << "<FONT COLOR=black> "<<levelnames(k)<<"</FONT>"<<endl;
    out << "<FONT COLOR=\"#990000\"> stationary state of the gauge field</FONT><br>"<<endl;
    out << "<FONT COLOR=\"#990000\">in the presence of a static source-antisource pair</FONT><br>"<<endl;
    out << "<FONT COLOR=\"#990000\"> for Run "<<run_name<<"</FONT><br><br>"<<endl;
    out << "</h2>"<<endl;
    out << "<p><br>"<<endl<<endl;

    for (i=1;i<=Norientations;i++){
       if (ans_stdev(k,i)>0.0){
          out << "<IMG SRC=\""<< levelstubs(k).right_substring('/')
              <<"_"<<Rfilenames(i) <<".gif\"><br><br><br><br>"<<endl;}}
       
    out << endl<<"</center>"<<endl;
    out << "<BR><BR><BR>"<<endl;
    out << "<a href=\"../"<<mainoutfile<<"\">"<<endl;
    out << "<img src=\"../../../images/back.gif\" align=middle border=0></a>"<<endl;
    out << "<a href=\"#top\">"<<endl;
    out << "<img src=\"../../../images/to_top.gif\" align=middle border=0></a>"<<endl;
    out << "</BODY>"<<endl;
    out << "</HTML>"<<endl;
    
    out.close();

    } 


 //  Now write the gap_R**.html files

 for (k=1;k<=Norientations;k++){

    fname=convert_to_charstar("Rlinks/gap_"+Rfilenames(k)+".html");
    out.clear();
    out.open(fname);
    delete [] fname;
    if (!out){
       cerr << " *** error opening file Rlinks/gap_"<<Rfilenames(k)<<".html"<<endl;
       return 1;}

    out << "<HTML>"<<endl;
    out << "<HEAD>"<<endl;
    out << "<TITLE> Fit results</TITLE>"<<endl;
    out << "</HEAD>"<<endl;
    out << "<BODY BGCOLOR=\"#99FFCC\">"<<endl;
    out << "<a href=\"../"<<mainoutfile<<"\">"<<endl;
    out << "<img src=\"../../../images/back.gif\" align=middle border=0></a>"<<endl;
    out << "<center>"<<endl;
    out << "<H2><FONT COLOR=\"#990000\">Fit results for the energy gap of the</FONT><br>"<<endl;
    out << "<FONT COLOR=\"#990000\"> stationary states of the gauge field</FONT><br>"<<endl;
    out << "<FONT COLOR=\"#990000\">in the presence of a static source-antisource pair</FONT><br>"<<endl;
    out << "<FONT COLOR=\"#990000\"> above the ground state for </FONT>"<<endl;
    out << "<FONT COLOR=black> R = "<<Rstrings(k)<<" </FONT>"<<endl;
    out << "<FONT COLOR=\"990000\"> for Run "<<run_name<<"</FONT><br><br>"<<endl;
    out << "</h2>"<<endl;
    out << "<p><br>"<<endl<<endl;

    for (i=1;i<=Ngaplevels;i++){
       if (gapans_stdev(i,k)>0.0){
          out << "<IMG SRC=\"../"<< gaplevelstubs(i)<<"_"<<Rfilenames(k) <<".gif\"><br><br><br><br>"<<endl;}}

    out << endl<<"</center>"<<endl;
    out << "<BR><BR><BR> "<<endl;
    out << "<a href=\"../"<<mainoutfile<<"\">"<<endl;
    out << "<img src=\"../../../images/back.gif\" align=middle border=0></a>"<<endl;
    out << "<a href=\"#top\">"<<endl;
    out << "<img src=\"../../../images/to_top.gif\" align=middle border=0></a>"<<endl;
    out << "</BODY>"<<endl;
    out << "</HTML>"<<endl;
    
    out.close();

    } 
    
 //  Now write the column gap-level.html files

 for (k=1;k<=Ngaplevels;k++){

    fname=convert_to_charstar(gaplevelstubs(k)+".html");
    out.clear();
    out.open(fname);
    delete [] fname;
    if (!out){
       cerr << " *** error opening file "<<gaplevelstubs(k)<<".html"<<endl;
       return 1;}

    out << "<HTML>"<<endl;
    out << "<HEAD>"<<endl;
    out << "<TITLE> Fit results</TITLE>"<<endl;
    out << "</HEAD>"<<endl;
    out << "<BODY BGCOLOR=\"#99FFCC\">"<<endl;
    out << "<a href=\"../"<<mainoutfile<<"\">"<<endl;
    out << "<img src=\"../../../images/back.gif\" align=middle border=0></a>"<<endl;
    out << "<center>"<<endl;
    out << "<H2><FONT COLOR=\"#990000\">Fit results for the energy gap of the</FONT><br>"<<endl;
    out << "<FONT COLOR=black> "<<gaplevelnames(k)<<"</FONT>"<<endl;
    out << "<FONT COLOR=\"#990000\"> stationary state of the gauge field</FONT><br>"<<endl;
    out << "<FONT COLOR=\"#990000\">in the presence of a static source-antisource pair</FONT><br>"<<endl;
    out << "<FONT COLOR=\"#990000\"> above the ground state for Run "<<run_name<<"</FONT><br><br>"<<endl;
    out << "</h2>"<<endl;
    out << "<p><br>"<<endl<<endl;

    for (i=1;i<=Norientations;i++){
       if (gapans_stdev(k,i)>0.0){
          out << "<IMG SRC=\""<< gaplevelstubs(k).right_substring('/')
              <<"_"<<Rfilenames(i) <<".gif\"><br><br><br><br>"<<endl;}}
       
    out << endl<<"</center>"<<endl;
    out << "<BR><BR><BR>"<<endl;
    out << "<a href=\"../"<<mainoutfile<<"\">"<<endl;
    out << "<img src=\"../../../images/back.gif\" align=middle border=0></a>"<<endl;
    out << "<a href=\"#top\">"<<endl;
    out << "<img src=\"../../../images/to_top.gif\" align=middle border=0></a>"<<endl;
    out << "</BODY>"<<endl;
    out << "</HTML>"<<endl;
    
    out.close();

    } 

  // ******************************
  //
  //  Now write the data files (useful for input to other plotting routines)
  //
  // ******************************
  
 Svector Rmags(Norientations),Rmagsinv(Norientations);
 for (i=1;i<=Norientations;i++){
    ans=calc_Rmagnitude(Rstrings(i),Nspace_dim);
    Rmags(i)=double_to_string(ans,6,10);
    Rmagsinv(i)=double_to_string(1.0/ans,6,10);
    }
 Ncol=0;
 for (i=1;i<=Norientations;i++)
    if (Rstrings(i).length()>Ncol) Ncol=Rstrings(i).length();
    
 fname=convert_to_charstar("all_results.dat");
 out.clear();
 out.open(fname);
 delete [] fname;
 if (!out){
    cerr << " *** error opening file all_results.dat"<<endl;
    return 1;}
 for (k=1;k<=Nlevels;k++){
    out << endl;
    out << " orientation    |R|        1/|R|             "<<levelnames(k)
        <<"       st_dev"<<endl<<endl;
    for (i=1;i<=Norientations;i++){
       if (ans_stdev(k,i)>0.0){
          out <<"  ";
          out.width(Ncol);
          out << Rstrings(i)<<"   "<<Rmags(i)<<"   "<<Rmagsinv(i)<<"      ";
          out << double_to_string(ans_mean(k,i),6,10)<<"  "
              << double_to_string(ans_stdev(k,i),6,10)<<endl;
          }}
    out << endl;    
    }
 out.close();
 
 fname=convert_to_charstar("all_gap_results.dat");
 out.clear();
 out.open(fname);
 delete [] fname;
 if (!out){
    cerr << " *** error opening file all_gap_results.dat"<<endl;
    return 1;}
 for (k=1;k<=Ngaplevels;k++){
    out << endl;
    out << " orientation    |R|        1/|R|             "<<gaplevelnames(k)
        <<"       st_dev"<<endl<<endl;
    for (i=1;i<=Norientations;i++){
       if (gapans_stdev(k,i)>0.0){
          out <<"  ";
          out.width(Ncol);
          out << Rstrings(i)<<"   "<<Rmags(i)<<"   "<<Rmagsinv(i)<<"      ";
          out << double_to_string(gapans_mean(k,i),6,10)<<"  "
              << double_to_string(gapans_stdev(k,i),6,10)<<endl;
          }}
    out << endl;    
    }
 out.close();
 
 return 0;
}

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

int main(){

 charString directory;
 charString mainoutfile,run_name;
 Svector labels;
 int Nirreps=0,Nspace_dim,orientation;
 double glueball_mass;

 if (read_filename(" directory for html file?", directory)) exit(-1);
 if (read_filename(" name of main html file?", mainoutfile)) exit(-1);
 mainoutfile=mainoutfile+".html";

 read_string(" name of run? ",run_name);
 read_int(" number of spatial dimensions? ",Nspace_dim);
 read_int(" orientation? (0=on-axis, 1=planar diagonal, 2=cubic diagonal) ",orientation);
 read_double(" a_t times glueball mass? ",glueball_mass);

 if (Nspace_dim==3){
    if (orientation==0){  // on-axis
       Nirreps=10;
       labels.reset(Nirreps);
       labels(1) = "A1g";
       labels(2) = "A1u";
       labels(3) = "A2g";
       labels(4) = "A2u";
       labels(5) = "B1g";
       labels(6) = "B1u";
       labels(7) = "B2g";
       labels(8) = "B2u";
       labels(9) = "Eg";
       labels(10)= "Eu";}
    else if (orientation==2){  // cubic-diagonal
       Nirreps=6;
       labels.reset(Nirreps);
       labels(1) = "A1g";
       labels(2) = "A1u";
       labels(3) = "A2g";
       labels(4) = "A2u";
       labels(5) = "Eg";
       labels(6) = "Eu";}
    else if (orientation==1){  // planar-diagonal
       Nirreps=8;
       labels.reset(Nirreps);
       labels(1) = "A1g";
       labels(2) = "A1u";
       labels(3) = "A2g";
       labels(4) = "A2u";
       labels(5) = "B1g";
       labels(6) = "B1u";
       labels(7) = "B2g";
       labels(8) = "B2u";}
    else{
       cerr << "bad orientation code"<<endl;
       exit(1);}
    }
 else if (Nspace_dim==2){
    Nirreps=4;      // 2+1 dimensional
    labels.reset(Nirreps);
    labels(1) = "Sg";
    labels(2) = "Su";
    labels(3) = "Ag";
    labels(4) = "Au";
    if ((orientation!=0)&&(orientation!=1)){
       cerr << "bad orientation code"<<endl;
       exit(1);}
    }

 write_pages(directory,mainoutfile,run_name,Nspace_dim,Nirreps,labels,glueball_mass);

 return 0;
}
// ***************************************************************