#include <stdio.h>
#include <string.h>
#include <time.h>
#include <stdlib.h>
#include <math.h>
#include "spectrum.h"
#include "jessica.h"
#include "extrapolate.h"
#if defined(USE_MPI)
#include "usempi.h"
#endif

#if HINDEX==2
#define PRINTIT 1 /*  Print timing and convergence information  */

/*
              Glueballs in 3+1 dimensions
              on the transverse lattice

    This program searches coupling constant space for a best fit
    to Poincare symmetries that are broken by the transverse lattice
    regulator.  

    We test for rotational invariance in the transverse plane.  
    In practice, we will choose the x^1 direction and the x^1=x^2 
    directions in order to take advandage of residual reflection 
    symmetries.

    This version also demands degeneracy in the lowest parity doublet.
    It also uses the string tension to determine the lattice spacing.


    We measure:                                                         multi

    p-truncation convergence for n=(1,0) winding mode                   15
    longitudinal potential for n=(0,0) and various L                     3
    longitudinal eigenvalue for n=(1,0) and one L                       15
    longitudinal eigenvalue for n=(1,1) and one L                       17
    winding mode potential for various n=(c,0)                          15
    winding mode potential for one n=(c,c)                              17
    spectrum for P_perp = (0,0) and (c,0)                               9, 8
                                                                        7, 10
    spectrum for P_perp = (c,c)                                         11, 12
                                                                        14, 13


************************************+++++++**********************************/

#define FTEPER 0  /* Number of Teper's lattice states to fit. */
                 /* This can be 0 or 1  */
#define FCSQUARED NPERPVALS /* Number of c^2 values in each sector and
                    direction to use in fit. This can be up to NPERPVALS */
#define FPARITY 3/* Number of parity doublets to fit. */
                 /* Suggested: 0, 2, 3, or 5. */
#define FSCALE 1 /*Whether to include consistantcy of the */
                 /*longitudinal potential with the overall scale.
		   The potential itself is calculated anyway.
		   This can be 0 or 1. */
#define FLONG MAX_LONG /* Number of longitudinal points to fit */
                /* This can be 0 to MAX_LONG. */
#define FWIND 1 /* number of winding spectra to fit
		   This can be 0 or 1. */
#define FEXTRA 0 /* Fit criteria that would not be normally used. */  
                   /*For instance, demanding 
                     a specific value for the lattice spacing.  */
#define NPERP 4 /* Number of J^{P_1 C} sectors for nonzero perp momentum 
		   this must be 4 to cover all multiplets */
#define NPERPVALS 4 /* number of states calculated in each sector */
#define NP1 1 /* Number of nonzero values of P_perp in the (c,0) direction*/
#define NP2 1 /* Number of nonzero values of P_perp in the (c,c) direction*/
#define NSTILL 10 /*Number of J^{P_1 C} sectors for zero transverse momentum spectra
                    This must be 10 to cover all multiplets. */ 
#define NSTILLVALS 4  /*Number of states calculated in each sector */
                      /*This must be large enough to use for nonzero momenta */
#define MAXTH 10 /* max number of terms in fit functions */
#define NLONG 3  /*Data points used to estimate longitudinal potential */
#define NWIND 3 /* Data points used to estimate winding potential */
#define MAX_LONG 3 /* max number of criteria for roundness of long. pot.*/


#if 0 /* smaller test basis */
/* This takes about 120 seconds to calculate everything */ 
#define NS 4
const static int kt[NS]={10,8,10,12}, np[NS]={4,6,6,4};
#define WNS 3
static int wkt0[NWIND][WNS]={{14,14,22},{11,11,21},{12,12,22}}, 
  wnp0[NWIND][WNS]={{6,4,4},{7,5,5},{8,6,6}},
  wht0[NWIND][HINDEX]={{2,0},{3,0},{4,0}};
static int wkt[1][WNS]={{14,14,22}},wnp[1][WNS]={{6,4,4}},
           wht[1][HINDEX]={{1,1}};
#define LNS 4
static int lkt0[LNS]={12,12,12,16}, lnp0[LNS]={2,4,2,2},
	     lkt[MAX_LONG][LNS]={{7,7,7,11},{7,7,7,11},
				{8,8,10,12}}, 
   lnp[MAX_LONG][LNS]={{3,5,3,3},{3,5,3,3},{2,4,2,2}};
doublereal lkmax[MAX_LONG][LNS]={{3.0,3.0,2.0,3.0},{3.0,3.0,2.0,3.0},
	{3.0,3.0,2.0,3.0}},lkmax0[LNS]={3.0,3.0,2.0,3.0},
  llong0[NLONG]={2.5,3.0,4.0},llong[MAX_LONG]={2.5,0.5,0.5},
  llong0m[NLONG]={0.1,0.1,0.1},llongm[MAX_LONG]={0.1,0.1,0.1};
#define PNS 4 /* number of sectors in p extrapolation */
static int pnp1[PNS]={1,3,5,7}, pkt1[PNS]={11,11,11,11};

#elif 0 /* small basis */
/* this takes about 300 seconds to calculate everything*/
#define NS 4
const static int kt[NS]={14,12,14,16}, np[NS]={4,6,6,4};
#define WNS 3
static int wkt0[NWIND][WNS]={{12,14,16},{9,11,13},{10,12,14}}, 
  wnp0[NWIND][WNS]={{6,4,4},{7,5,5},{8,6,6}},
  wht0[NWIND][HINDEX]={{2,0},{3,0},{4,0}};
static int wkt[1][WNS]={{12,14,16}},wnp[1][WNS]={{6,4,4}},
    wht[1][HINDEX]={{2,2}};
#define LNS 4
static int lkt0[LNS]={12,12,12,16}, lnp0[LNS]={2,4,2,2},
  lkt[MAX_LONG][LNS]={{7,7,7,11},{7,7,7,11},
    {8,8,10,12}}, 
    lnp[MAX_LONG][LNS]={{3,5,3,3},{3,5,3,3},{2,4,2,2}};
doublereal lkmax0[LNS]={3.0,3.0,2.0,3.0},
      lkmax[MAX_LONG][LNS]={{3.0,3.0,2.0,3.0},
        {3.0,3.0,2.0,3.0},{3.0,3.0,2.0,3.0}},
	llong0[NLONG]={3.0,4.0,5.0},llong[MAX_LONG]={3.0,1.0,1.0},
	llong0m[NLONG]={0.1,0.1,0.1},llongm[MAX_LONG]={0.1,0.1,0.1};
#define PNS 6 /* number of sectors in p extrapolation */
static int pnp1[PNS]={3,3,5,5,7,7},pkt1[PNS]={7,9,7,9,7,9};

#elif 1 /* larger basis */
/* this takes about 25 minutes to calculate everything*/
#define NS 4
const static int kt[NS]={12,12,14,18}, np[NS]={8,6,6,6};
#define WNS 3
#if GLUE_PERIODIC
static int wkt0[NWIND][WNS]={{16,16,26},{18,18,26},{16,16,26}};
#else 
static int wkt0[NWIND][WNS]={{16,16,26},{17,17,25},{16,16,26}}; 
#endif
static int wnp0[NWIND][WNS]={{6,4,4},{7,5,5},{8,6,6}},
  wht0[NWIND][HINDEX]={{2,0},{3,0},{4,0}};
static int wkt[1][WNS]={{14,14,22}},wnp[1][WNS]={{8,6,6}},
    wht[1][HINDEX]={{2,2}};
#define LNS 7
static int lkt0[LNS]={26,26,26,42,42,70,70}, lnp0[LNS]={2,4,2,2,2,2,2};
#if GLUE_PERIODIC
static int lkt[MAX_LONG][LNS]={{14,14,14,24,24,36,36},
		     {14,14,14,24,24,36,36},
			       {22,22,22,44,44,70,70}};
#else
static int lkt[MAX_LONG][LNS]={{13,13,13,23,23,35,35},
		     {13,13,13,23,23,35,35},
			       {22,22,22,44,44,70,70}}; 
#endif
static int lnp[MAX_LONG][LNS]={{3,5,3,3,3,3,3},{3,5,3,3,3,3,3},
		       {2,4,2,2,2,2,2}};
doublereal lkmax0[LNS]={3.5,3.5,3.0,3.5,4.25,4.0,4.5},
      lkmax[MAX_LONG][LNS]={{3.5,3.5,3.0,4.0,4.5,5.0,6.0},
       {3.5,3.5,3.0,4.0,4.5,5.0,6.0},{5.0,5.0,4.5,5.0,6.0,6.0,7.0}},
      llong0[NLONG]={3.577,5.058,7.154},llong0m[NLONG]={-2.726,-3.856,-5.452},
      llong[MAX_LONG]={1.265,5.058,0.0},llongm[MAX_LONG]={-0.964,-3.856,0.0};
#define PNS 7 /* number of sectors in p extrapolation */
static int pnp1[PNS]={3,3,3,5,5,7,7};
#if GLUE_PERIODIC
static int pkt1[PNS]={14,22,46,14,28,14,16};
#else
static int pkt1[PNS]={13,21,45,13,27,13,15};
#endif

#elif 1 /* BIG basis */
/* this takes about 6.654 hours to calculate everything on PC's,
   4.14 CPU-hours on the CRAY T3E.  
   max_els=3101540, max_base=2255950 (but this could grow...)  */
#define NS 4
const static int kt[NS]={16,16,20,26}, np[NS]={8,6,6,6};
#define NWIND 3
#define WNS 4
#if GLUE_PERIODIC
static int wkt0[NWIND][WNS]={{20,20,24,32},{20,20,26,34},{18,18,24,32}}; 
#else
static int wkt0[NWIND][WNS]={{20,20,24,32},{19,19,25,33},{18,18,24,32}}; 
#endif
static int wnp0[NWIND][WNS]={{4,6,4,4},{5,7,5,5},{6,8,6,6}},
  wht0[NWIND][HINDEX]={{2,0},{3,0},{4,0}},wkt[1][WNS]={{16,16,24,28}},
    wnp[1][WNS]={{6,8,6,6}},wht[1][HINDEX]={{2,2}};
#define LNS 7
static int lkt0[LNS]={40,40,40,50,50,70,70},
  lnp0[LNS]={2,4,2,2,2,2,2};
#if GLUE_PERIODIC
static int lkt[MAX_LONG][LNS]={{20,20,20,34,34,70,70},
			 {20,20,20,34,34,70,70},
			       {34,34,34,48,48,70,70}};
#else 
static int lkt[MAX_LONG][LNS]={{19,19,19,33,33,69,69},
			 {19,19,19,33,33,69,69},
			       {34,34,34,48,48,70,70}};
#endif
int lnp[MAX_LONG][LNS]={{3,5,3,3,3,3,3},{3,5,3,3,3,3,3},{2,4,2,2,2,2,2}};
doublereal lkmax0[LNS]={4.0,4.0,3.5,3.75,4.25,4.0,4.5},
      lkmax[MAX_LONG][LNS]={{4.0,4.0,3.5,4.0,5.0,5.0,6.0},
       {4.0,4.0,3.5,4.0,5.0,5.0,6.0},{5.5,5.5,5.0,5.5,6.5,6.0,7.0}},
      llong0[NLONG]={3.577,5.058,7.154},llong0m[NLONG]={-2.726,-3.856,-5.452},
      llong[MAX_LONG]={1.265,5.058,0.0},llongm[MAX_LONG]={-0.964,-3.856,0.0};
#define PNS 7 /* number of sectors in p extrapolation */
static int pnp1[PNS]={3,3,3,5,5,7,7};
#if GLUE_PERIODIC
static int pkt1[PNS]={14,36,46,14,30,14,18};
#else
static int pkt1[PNS]={13,35,45,13,29,13,17};
#endif

#endif  /** Choose different possible bases **/

/* Total number of different bases */
#define N_ALL_BASES ((NSTILL+(NP1+NP2)*NPERP)*NS+(NWIND+FWIND)*WNS+\
                    (FLONG+NLONG)*LNS+PNS)
static specify_calculation all_bases[N_ALL_BASES];
static specify_calculation *s[NP1][NPERP],*t[NP2][NPERP],*ws0[NWIND],*ws[FWIND],
  *ls[FLONG],*ls0[NLONG],*ps1,*stillbases[NSTILL];

element par[NPARAMS];

             /* This is the lattice value for the 0^++ in units of the 
                string tension extrapolated to large N.  I assume that
                the large N extrapolation is the same as for the 2+1 glueball.
                The following is the mass formula:
                 1648./440(1-1/4.11* 2.379/3^2)
                The problem here is that the errors are too large to
                really constrain the value for anything other than
		the lowest state.  */
#if FTEPER<2
doublereal teper[1]={12.67},tepererr[1]={1.28};
#endif

static int wpar[NPARAMS],scalemethod;
static doublereal p1a[NP1][HINDEX]={{0.25,0.0}},
  p2a[NP2][HINDEX]={{0.25,0.25}};
static int pht1[HINDEX]={1,0};  /* winding number for p-extrapolation */

/* The error for: 
     c^2 for each state (either direction), 
     roundness of winding mode potential (two data points), and
     roundness of heavy quark potential,
     parity doublets  (fractional error).           */

#if FCSQUARED <= 4  /* check that array is big enough */
static const element c2err[NPERP][4]={{0.25,2.0,2.0,2.0}, {0.25,2.0,2.0,2.0},
		    {0.25,2.0,2.0,2.0}, {0.25,2.0,2.0,2.0}},
#endif

#if 1 && MAX_LONG==3 /*Normal values for longitudinal and winding errors */
                    /* these errors match what was found in 2+1 ..  */
  lerr0=0.1, lerr[MAX_LONG]={0.1,0.1,0.2}, werr[1]={1.0};
#elif 0 && MAX_LONG==3/*Normal values for longitudinal and winding errors */
  lerr0=0.2, lerr[MAX_LONG]={0.1,0.1,0.2}, werr[1]={1.0};
#elif 0 && MAX_LONG==3 /*larger errors for longitudnal and winding criteria. */
  lerr0=1.0, lerr[MAX_LONG]={1.0,1.0,1.0}, werr[1]={1.0};
#endif

/* These are the parity doublets that are to be compared.
   parityerror is the fractional error assumed for each  
   doublet.
   each element of whichparity contains {charge,multiplet,which state}
   for each parity doublet member. */
#if FPARITY<=5  /* This is just the size of the arrays. */
  static element parityerr[5]={0.25,0.25,0.25,0.25,0.25};
  static struct {int charge; int multi; int state;}
  whichparity[5][2]={{{1,5,0},{1,6,0}},{{-1,5,0},{-1,6,0}},
		      {{-1,7,0},{-1,4,0}},
		      {{1,7,0},{1,5,0}},{{1,7,0},{1,6,0}}};
#endif

/* Error for extra fit criteria. */

#if FEXTRA<2
doublereal extraerr[1]={0.1};
#endif

/* Winding numbers for longitudinal string tension. lht0 must be
   zero since the lattice spacing is unknown.  Also the 
   winding numbers and lattice symmetry multiplet for determining
   roundness of the heavy sourc potential.  */
#if MAX_LONG==3
static int lht[MAX_LONG][HINDEX]={{1,0},{1,0},{1,1}},
  lmulti[MAX_LONG]={15,15,17},wmulti[1]={17};
#endif
static integer wth,lth;
/* These are all the quantities returned by fcn */
static struct {doublereal pextrap; element lvalue[FLONG]; 
  element wvalue[FWIND]; 
  element calclvalue[FLONG]; element calcwvalue[FWIND];
  element wfitting[MAXTH]; element lfitting[MAXTH]; element rescale;
  element spectrumvals1[NP1+1][NPERP][NPERPVALS]; 
  element spectrumvals2[NP2+1][NPERP][NPERPVALS];
  element stillvals[NSTILL][NSTILLVALS]; element c2[2][NPERP][NPERPVALS];
} best;
   
#define CHI_DIM ((NP1+NP2)*NPERP*FCSQUARED+FLONG+FWIND+FPARITY+FTEPER+FSCALE+FEXTRA) 
                    /*Number of dimensions of chi^2 fit*/
#define NOUT 100 /* Output file name length */

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

int MAIN(int argc, char **argv){
  char out[NOUT],*pout;
  /* maxfev is the maximum number of iterations in minimizing chi^2. */
  integer mm=CHI_DIM,npar=1,maxfev,mode=1,ifail=0,info,nfev,ipvt[NPARAMS];
  /*  epsfcn should be a little less than the Lanczos accuracy. (say 10 times)
      ftol should be a little more than epsfcn
      xtol is the desired precision of the coupling constants  
      factor controls the maximum initial step size.  If the
        starting parameters are near the minimum, this must
        be decreased to maybe 1. 
        (recommended range:  0.1 to 100, default 100). */
  /*  epsfcn is the variable used to determine the step size 
      for the finite difference approximation.  Since the 
      extrapolations magnify any error in the eigenvalues,
      this should be comfortably larger than the errors in
      the lanczos routine.  However, it should also be smaller
      than xtol.                                             */
  doublereal xtol=0.01,gtol=0.0,ftol=5.0e-3,epsfcn=5.0e-5,factor=1.0;
  doublereal vpar[NPARAMS],fvec[CHI_DIM],diag[NPARAMS],fjac[CHI_DIM*NPARAMS];
  doublereal qtf[NPARAMS],wa1[NPARAMS],wa2[NPARAMS],wa3[NPARAMS],wa4[CHI_DIM];
  int count_bases=0;
  element angle[HINDEX+2];
  element chi2;
#if NSTILL==6
  int multistill[NSTILL]={3, 5, 5, 6, 6, 7},
    /*    ostill[NSTILL]={1, 1,-1, 1,-1,-1};*/
    charge_still[NSTILL]={1, 1,-1, 1,-1, 1};
#elif NSTILL==10
  int multistill[NSTILL]={3, 3, 4, 4, 5, 5, 6, 6, 7, 7},
    /* ostill[NSTILL]=   {1,-1, 1,-1, 1,-1, 1,-1, 1,-1};*/
    charge_still[NSTILL]={1,-1, 1,-1, 1,-1, 1,-1,-1, 1};
#elif NSTILL==12
  int multistill[NSTILL]={3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8},
    /*    ostill[NSTILL]={1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1};*/
    charge_still[NSTILL]={1,-1, 1,-1, 1,-1, 1,-1,-1, 1,-1, 1};
#endif
  int ht_zero[HINDEX]={0,0}; /* zero winding for ordinary states */
  int i,j,k,l,
    multis[NPERP]={15,15,16,16},
    chargest[NPERP]={1,-1,-1,1},
      multit[NPERP]={17,17,18,18};
  FILE *fp;

#ifdef USE_MPI
  int myrank;
#endif

#ifdef BABBAGE  /* Routine to initialize Fortran on Babbage */
  if(hf_fint((char*)NULL)) exit(1);
#endif

  /* initialize parallelization via MPI when applicable.
     in this case we use a simple master/slave setup to
     distribute the construction of the Hamiltonian and the
     diagonalization of the various bases.      */

#ifdef USE_MPI
  MPI_Init(&argc, &argv);     /* initialize MPI */
  initmpisectors();
  MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
  if(myrank!=WORLD_MASTER){
    slave();
    MPI_Finalize();   /* cleanup MPI */
    return 0;
  }
#endif

  /* Choose scale method and open output file */

  printf("Output file name:  ");
  fgets(out,NOUT,stdin);
  if((pout=strchr(out,'\n'))!=NULL)*pout='\0';
  printf("%s\n",out);
  if((fp = fopen(out,"w")) == NULL){
    fprintf(stderr,__FILE__":  Can't open file %s, exiting\n",out);
    goto cleanup;
  }

  printf("Input maximum number of iterations:  ");
  scanf("%li",&maxfev);

  printf("Choose method for determining the overall scale:\n"
         "    0 use longitudinal string tension\n"
         "    1 fit lowest state in spectrum to lattice value\n"
         "    2 demand minimum chi^2 for c^2 measurements\n"
	 "    3 find scale by demanding overall minimum chi^2.\n");
  scanf("%i",&scalemethod);  
  
  /*  Print out the results in the output file **/
  fprintf(fp,"********* Eigenvalues for the 3+1 transverse lattice *********\n");
  fprintf(fp,"Couplings:\n    ");
  for(i=0; i<NPARAMS; i++){
    fprintf(fp,"  %i:%s",i,coupling_string(i));
    if(i%5==4 && i!=NPARAMS-1)fprintf(fp,"\n    ");
  }    
  fprintf(fp,"\nUse chi^2 fit with %i criteria and tolerance %g.\n",
	  CHI_DIM,xtol);
  switch(scalemethod){
  case 0:
    fprintf(fp,"Overall scale from longitudinal string tension.\n");
    break;
  case 1:
    fprintf(fp,"Overall scale from fitting lowest state to lattice value.\n");
    break;
  case 2:
    fprintf(fp,"Overall scale determined from best fit for c^2.\n");
    break;
  case 3:
    fprintf(fp,"Overall scale from minimizing chi^2.\n");
    break;
  default:
     fprintf(stderr,__FILE__":  Incorrect scale method = %i, exiting.\n",
	     scalemethod);
    goto cleanup;
  }
#ifdef KSQUARED
  fprintf(fp,"1/K^2 term used in all extrapolations.\n");
#endif
#if FEXTRA>0
  fprintf(fp,"Includes %i nonstandard fit criteria (such as a \n"
	  " specific lattice spacing).\n",FEXTRA);
#endif
#if FTEPER>0
  fprintf(fp,"Includes a fit to %i eigenvalues from Teper.\n",FTEPER);
#endif
#if FPARITY>0
  fprintf(fp,"Parity doublets (charge,multi,state) with error for difference\n"
	  "  over average.\n");
  for(i=0; i<FPARITY; i++)
    fprintf(fp,"  (%2i,%2i,%i) and (%2i,%2i,%i), error %g\n",
	    whichparity[i][0].charge,whichparity[i][0].multi,
	    whichparity[i][0].state,
	    whichparity[i][1].charge,whichparity[i][1].multi,
	    whichparity[i][1].state,parityerr[i]);
#else
  fprintf(fp,"No parity doublets.\n");
#endif
  fprintf(fp,"Spectrum for P_perp a = (0,0),");
  for(i=0; i<NP1; i++){
    fprintf(fp," (");
    for(j=0; j<HINDEX; j++)fprintf(fp," %g",p1a[i][j]);
    fprintf(fp,")");
  }
  fprintf(fp,"\n  using %i states (multiplet, charge, c^2 error for each) =",
	  NPERPVALS);
  for(i=0; i<NPERP; i++){
    fprintf(fp,"\n  (%2i, %2i,",multis[i],chargest[i]);
    for(j=0; j<FCSQUARED; j++)
      fprintf(fp," %g",c2err[i][j]);
    fprintf(fp,")");
  }
  fprintf(fp,"\nSpectrum for P_perp a = (0,0),");
  for(i=0; i<NP2; i++){
    fprintf(fp," (");
    for(j=0; j<HINDEX; j++)fprintf(fp," %g",p2a[i][j]);
    fprintf(fp,")");
  }
  fprintf(fp,"\n  using %i states (multiplet, charge, c^2 error for each) =",
	  NPERPVALS);
  for(i=0; i<NPERP; i++){
    fprintf(fp,"\n  (%2i, %2i,",multit[i],chargest[i]);
    for(j=0; j<FCSQUARED; j++)
      fprintf(fp," %g",c2err[i][j]);
    fprintf(fp,")");
  }
#if NSTILL>0
  fprintf(fp,".\nOrdinary spectra multiplets, %i states, (charge,multiplet) =",
	  NSTILLVALS);
  for(i=0, k=5; i<NSTILL; k++, i++){
    if(k%6==0)fprintf(fp,"\n ");
    fprintf(fp," (%2i, %2i) ",charge_still[i],multistill[i]);
  }
#endif
  fprintf(fp,".\nAll spectra extrapolated using (K,p) =");
  for(i=0, k=3; i<NS; k++, i++){
    if(k%5==0)fprintf(fp,"\n ");
    fprintf(fp," (%i/2,%i) ",kt[i],np[i]);
  }
  fprintf(fp,".\nWinding potential using (n,K,p) =");
  for(i=0, k=2; i<NWIND; i++)
    for(j=0; j<WNS; k++, j++){
      if(k%4==0)fprintf(fp,"\n ");
      fprintf(fp," (");
      for(l=0; l<HINDEX; l++)fprintf(fp," %i",wht0[i][l]);
      fprintf(fp,",%i/2,%i) ",wkt0[i][j],wnp0[i][j]);
    }
  fprintf(fp,".\nRoundness of winding using (n,K,p) =");
  for(i=0, k=3; i<FWIND; k=0, i++){
    for(j=0; j<WNS; k++, j++){
      if(k%4==0)fprintf(fp,"\n ");
      fprintf(fp," (");
      for(l=0; l<HINDEX; l++)fprintf(fp," %i",wht[i][l]);
      fprintf(fp,",%i/2,%i) ",wkt[i][j],wnp[i][j]);
    }
    fprintf(fp,"\n  with error %g;",werr[i]);
  }
  fprintf(fp," in G^2 N units.\n"
	  "Heavy potential determined using (n,K,p,K_max) =");
  for(j=0, k=3; j<LNS; k++, j++){
    if(k%3==0)fprintf(fp,"\n ");
    fprintf(fp," (");
    for(l=0; l<HINDEX; l++)fprintf(fp," %i",ht_zero[l]);
    fprintf(fp,",%i/2,%i,%g) ",lkt0[j],lnp0[j],lkmax0[j]);
  }
  fprintf(fp,",\n L =");
  for(i=0; i<NLONG; i++)fprintf(fp," %g%+g*m",llong0[i],llong0m[i]);
  fprintf(fp,"\n (all in G^2 N units)");
  if(FSCALE)fprintf(fp,"; relative scale error=%g",lerr0);
  fprintf(fp,".\nRoundness determined using (n,K,p,K_max) =");
  for(i=0, k=2; i<FLONG; k=0, i++){
    for(j=0; j<LNS; k++, j++){
      if(k%3==0)fprintf(fp,"\n ");
      fprintf(fp," (");
      for(l=0; l<HINDEX; l++)fprintf(fp," %i",lht[i][l]);
      fprintf(fp,",%i/2,%i,%g) ",lkt[i][j],lnp[i][j],lkmax[i][j]);
    }
    fprintf(fp,"\n L=%g%+g*m and error %g;",
	    llong[i],llongm[i],lerr[i]);
  }
  fprintf(fp," all in G^2 N units.\np-extrapolation using n=(");
  for(l=0; l<HINDEX; l++)fprintf(fp," %i",pht1[l]);
  fprintf(fp,") and (K,p) =");
  for(j=0, k=4; j<PNS; k++, j++){
    if(k%5==0)fprintf(fp,"\n ");
    fprintf(fp," (%i/2,%i) ",pkt1[j],pnp1[j]);
  }
  fprintf(fp,".\nResult format:\n"
	  " Fit info, # steps, chi^2, and p damping, scale G^2 N/sigma.\n"
	  " The %i couplings  (G^2 N units);\n"
	  "   which couplings -- if any -- were fit.\n"
	  " Winding potential and heavy souce potential fits.\n"
	  " Roundness of winding potential, %i values, and roundness of\n"
	  "   heavy source potential, %i values, (G^2 N units)\n"
	  "   showing measured value and derived value for each.\n"
	  " The rescaled spectrum for each P_perp*a and c^2 values.\n"
#if NSTILL>0
	  " Finally come the states for the ordinary spectra.\n"
#endif
	  "\n",NPARAMS,FWIND,FLONG);
  fflush(fp);

  /******************* make basis structures  ****************************/
  for(k=0; k<HINDEX+2; k++)angle[k]=0;

  printf("Make %i spectrum bases,",(NP1+NP2)*NPERP*NS);
  /*  Calculate the extrapolated spectrum in the P_perp = (c,0) direction */
  for(l=0; l<NP1; l++){
    for(k=0; k<HINDEX; k++)angle[k]=p1a[l][k];
    for(j=0; j<NPERP; j++){
      s[l][j]=all_bases+count_bases;
      for(i=0; i<NS; i++)
	initspectrum(all_bases+count_bases++,TYPE_GLUE,ht_zero,np[i],kt[i],
		     chargest[j],0,multis[j],1,angle,NPERPVALS+1,0);
    }
  }
  /*  Calculate the extrapolated spectrum in the P_perp = (c,c) direction */
  for(l=0; l<NP2; l++){
    for(k=0; k<HINDEX; k++)angle[k]=p2a[l][k];
    for(j=0; j<NPERP; j++){
      t[l][j]=all_bases+count_bases;
      for(i=0; i<NS; i++)
	initspectrum(all_bases+count_bases++,TYPE_GLUE,ht_zero,np[i],kt[i],
		     chargest[j],0,multit[j],1,angle,NPERPVALS+1,0);
    }
  }
 
  /*  Calculate the ordinary extrapolated spectra. */
  printf(" %i ordinary spectrum bases,",NSTILL);
  for(j=0; j<NSTILL; j++){
    stillbases[j]=all_bases+count_bases;
    for(i=0; i<NS; i++)
      initspectrum(all_bases+count_bases++,TYPE_GLUE,ht_zero,np[i],kt[i],
		   charge_still[j],0,multistill[j],0,NULL,NSTILLVALS,0);
  }

  printf(" %i winding bases,",(NWIND+FWIND)*WNS);
  for(i=0; i<NWIND; i++){
    ws0[i]=all_bases+count_bases;
    for(j=0; j<WNS; j++)
      /* assume o=1 is lowest state */
      initspectrum(all_bases+count_bases++,TYPE_GLUE,wht0[i],wnp0[i][j],
		   wkt0[i][j],0,1,15,0,NULL,1,0);
  }

  for(i=0; i<FWIND; i++){
    ws[i]=all_bases+count_bases;
    for(j=0; j<WNS; j++)
      /* assume o=1 is lowest state */
      initspectrum(all_bases+count_bases++,TYPE_GLUE,wht[i],wnp[i][j],
		   wkt[i][j],0,1,wmulti[i],0,NULL,1,0);
  }

  /*  Calculate the longitudinal string tension  */
  printf("\n %i long bases,",(FLONG+1)*LNS);
  for(i=0; i<NLONG; i++){
    ls0[i]=all_bases+count_bases;
    /* orientation reversal is not a good symmetry for L != 0 */
    for(j=0; j<LNS; j++){
      angle[HINDEX]=lkmax0[j];
      initspectrum(all_bases+count_bases++,TYPE_SOURCE,ht_zero,lnp0[j],lkt0[j],
		   0,0,3,1,angle,1,0);
    }
  }
  
  /*  Calculate potential for n=(c,0), L nonzero  */
  for(i=0; i<FLONG; i++){
    /* orientation reversal is not a good symmetry for L != 0 */
    ls[i]=all_bases+count_bases;
    for(j=0; j<LNS; j++){
      angle[HINDEX]=lkmax[i][j];
      initspectrum(all_bases+count_bases++,TYPE_SOURCE,lht[i],lnp[i][j],
		   lkt[i][j],0,0,lmulti[i],1,angle,1,0);
    }
  }
  
  /* Calculate the convergence in p-truncation */
  printf(" and %i n=%i bases",PNS,*pht1);
  ps1=all_bases+count_bases;
  for(j=0; j<PNS; j++)
    /* assume o=1, multi=15 is lowest state */
    initspectrum(all_bases+count_bases++,TYPE_GLUE,pht1,pnp1[j],pkt1[j],
		 0,1,15,0,NULL,1,0);
  
  printf("\n");
  fflush(stdout);

  if(count_bases!=N_ALL_BASES){
    fprintf(stderr,__FILE__":%i  Number of bases don't match\n",__LINE__);
    goto cleanup;
  }


  /******************* Main loop of program ********************************/

  for(;;){
    printf("How many to fit (0 to %i, -1=stop)\n",NPARAMS-2);
    scanf("%li",&npar);
    if(npar<0)break;
#if 1 /* Input starting parameters */
    printf("Input %i starting parameters\n",NPARAMS);
    for(i=0; i<NPARAMS; i++)scanf("%lf",par+i);
#endif

    if(npar==0){
      info=0; nfev=1;
      fcn(&mm,&npar,vpar,fvec,&info);
    }
    else{
      printf("Input which %li couplings (0=m^2, ...)\n",npar);
      for(i=0;i<npar;i++)scanf("%i",wpar+i);
      for(i=0;i<npar;i++)vpar[i]=par[wpar[i]];

/*    Least squares fitting routine    
      subroutine lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,
 *                  maxfev,epsfcn,diag,mode,factor,nprint,info,nfev,
 *                 fjac,ldfjac,ipvt,qtf,wa1,wa2,wa3,wa4)  */
      LMDIF(FCN_TO_FORTRAN,&mm,&npar,vpar,fvec,&ftol,&xtol,&gtol,
            &maxfev,&epsfcn,diag,&mode,&factor,&ifail,&info,&nfev,
            fjac,&mm,ipvt,qtf,wa1,wa2,wa3,wa4);

      if(info<1||info>4) fprintf(stderr,__FILE__":  Error in lmdif info=%li\n",info);
      if(info==0)goto cleanup;
      /*Map best values into par[] */
      for(i=0; i<npar; i++)par[wpar[i]]=vpar[i];
    }
    if(info<0){
       fprintf(stderr,__FILE__":  Error in fcn info=%li\n",info);
      if(info<=-8)goto cleanup;
    }
    chi2=fdot(fvec,fvec,mm);
    printf("error %li, %li steps chi2=%f with c_p=%f rescale=%f\n",
	   info,nfev,chi2,best.pextrap,best.rescale);
    fprintf(fp,"%li %li " PDIGITS " " PDIGITS " " PDIGITS "\n",
	    info,nfev,chi2,best.pextrap,best.rescale);
    printf("Final couplings: ");
    for(i=0; i<NPARAMS; i++){
      printf(" " PDIGITS,par[i]); fprintf(fp,PDIGITS " ",par[i]);
    }
    if(npar>0){fprintf(fp,"\n");printf("\nwhile fitting");}
    for(i=0; i<npar; i++){
      printf(" %i",wpar[i]); fprintf(fp," %i",wpar[i]);
    }
    printf("\nWinding: "); fprintf(fp,"\n");
    for(i=0; i<wth; i++){
      printf("%f ",best.wfitting[i]); fprintf(fp,"%f ",best.wfitting[i]);
    }
    printf("and long:");
    for(i=0; i<lth; i++){
      printf(" %f",best.lfitting[i]); fprintf(fp," %f",best.lfitting[i]);
    }
    printf("\n Roundness of winding: ");
    fprintf(fp,"\n");
    for(i=0; i<FWIND; i++){
      printf("%f %f  ",best.wvalue[i],best.calcwvalue[i]);
      fprintf(fp,"%f %f  ",best.wvalue[i],best.calcwvalue[i]);
    }
    printf("\n Roundness of potential: ");
    for(i=0; i<FLONG; i++){
      printf("%f %f  ",best.lvalue[i],best.calclvalue[i]);
      fprintf(fp,"%f %f  ",best.lvalue[i],best.calclvalue[i]);
    }
    printf("\n");
    fprintf(fp,"\n");
    for(i=0; i<NPERP; i++)
      for(k=0; k<NPERPVALS; k++){
	for(j=0; j<NP1+1; j++){
	  printf(" %f",best.rescale*best.spectrumvals1[j][i][k]); 
	  fprintf(fp," %f",best.rescale*best.spectrumvals1[j][i][k]);
	}
	printf(" %f\n",best.c2[0][i][k]); fprintf(fp," %f\n",best.c2[0][i][k]); 
      }
    for(i=0; i<NPERP; i++)
      for(k=0; k<NPERPVALS; k++){
	for(j=0; j<NP2+1; j++){
	  printf(" %f",best.rescale*best.spectrumvals2[j][i][k]); 
	  fprintf(fp," %f",best.rescale*best.spectrumvals2[j][i][k]);
	}
	printf(" %f\n",best.c2[1][i][k]); fprintf(fp," %f\n",best.c2[1][i][k]); 
      }
    for(i=0; i<NSTILL; i++){
      printf("charge=% i multi=%i ",charge_still[i],multistill[i]);
      for(j=0; j<NSTILLVALS; j++){
	printf("%f ",best.rescale*best.stillvals[i][j]); 
	fprintf(fp,"%f ",best.rescale*best.stillvals[i][j]);
      }
      fprintf(fp,"\n"); printf("\n");
    }
    fprintf(fp,"\n"); printf("\n");
    fflush(fp);
    fflush(stdout);
  }
  
  
cleanup:
  fclose(fp);

  /* Free up memory used.  This is mainly for debugging */
  free_spectrum();

#if defined(USE_MPI)
  /************* Tell slaves to shut down. **************/
#if ONE_PROCESSOR_SPECTRUM
  MPI_Comm_size(MPI_COMM_WORLD,&myrank);
  for(i=1; i<myrank; i++)
    /* only the message tag=0 is read by the slaves. */
    MPI_Send(&myrank,1,MPI_INT,i,0,MPI_COMM_WORLD);
#else
  MasterShutDown();
#endif
  MPI_Finalize();             /* cleanup MPI */
#else
  printf("Maximum matrix size: n=" SIZE_MEM_PRINT ", nonzero=" SIZE_MEM_PRINT 
	 ", " SIZE_MEM_PRINT " bytes\n",
	 max_els.n,max_els.nonzero,max_els.memory);
  printf("Memory used to construct bases: " SIZE_MEM_PRINT " bytes\n",
	 reduced_basis_memory);
#endif

  return 0;
}


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

  The following routine calculates the fit criteria for
  a given set of fit parameters.  It returns the vector
  fvecc which is then squared to give the chi^2.

  The values that are calculated are reused from previous
  iterations if the relevant couplings are unchanged.
  This includes:

         pextrap, wfitting, lfitting, wvalue1, lvalue1,
	 lvalue2, spectrumvals1, spectrumvals2, stillvals
         
  Information about the spectra are returned in the structure
  best.  The members of the structure best 
  are from the lowest chi^2 iteration of fcn (which is
  not necessarily the last iteration of fcn.)

  On error, fcn sets *iflag to a negative integer.
  *flag<=-8 are "fatal" errors.

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


extern fsub fcn(integer *m, integer *npar, doublereal *xp, 
	       doublereal *fvecc, integer *iflag){
  /* fitc is the index for the fit criterion */
  int i,k,j,l,fitc;
  integer nwind=NWIND,nlong=NLONG,th[1];
  element *vals=NULL; /* this is realloc'ed by the extrapolation routines */
  element sum,temp,rescale,c2[2][NPERP][NPERPVALS],zz;
  doublereal mlong0[NLONG],mlong[MAX_LONG];
  static element oldpar[NPARAMS],bestpar[NPARAMS],bestchi2;
  static doublereal spectrumvals1[NP1+1][NPERP][NPERPVALS],
    spectrumvals2[NP2+1][NPERP][NPERPVALS],
    stillvals[NSTILL][NSTILLVALS];
  static element pextrap,calclvalue[FLONG],calcwvalue[FWIND];
  static element wfitting[MAXTH],lfitting[MAXTH],lvalue[FLONG],
    wvalue[FWIND],wvalue20;
  
  clock_t tif;
#if NWIND>NLONG
  element weights[NWIND*3];
#else
  element weights[NLONG*3];
#endif

  int job;
  element *all_vals[N_ALL_BASES];
#if defined(USE_MPI) && ONE_PROCESSOR_SPECTRUM
#else
  full_basis basis={0};
  eigensystem eigen=INITIAL_EIGENSYSTEM;
#endif
  tif=clock();

  /* bestpar is used later to see if unfitted variables have changed*/
  for(i=0; i<*npar; i++)bestpar[wpar[i]]=par[wpar[i]]=xp[i];

#if PRINTIT
  printf("  Parameters: ");
  for(i=0; i<NPARAMS; i++)printf(" %.10g",par[i]);
  printf("\n");
#endif

  /* Find longitudinal values for a given mass */

  for(i=0;i<nlong;i++){
    mlong0[i]=llong0[i]+sqrt(par[0]/par[1])*llong0m[i];
      for(j=0; j<LNS; j++)
	ls0[i][j].angle[HINDEX+1]=mlong0[i];
  }

  for(i=0;i<FLONG;i++){
    mlong[i]=llong[i]+sqrt(par[0]/par[1])*llongm[i];
    for(j=0; j<LNS; j++)
	ls[i][j].angle[HINDEX+1]=mlong[i];
  }

  /* allocate memory to store eigenvalues */

  for(i=0; i<N_ALL_BASES; i++)
    all_vals[i]=malloc(all_bases[i].nval*sizeof(element));


#if defined(USE_MPI) && ONE_PROCESSOR_SPECTRUM

  /************ now distribute the work to the various slaves **********/

#if 0 /* debugging print */
  printf("starting subroutine master with %i jobs.\n",njob);
#endif
  if(master(mpibase,par,angle,mpinvals,all_vals,mpivals,njob)){
     fprintf(stderr,__FILE__":  master() is returning an error, exiting\n");
    *iflag=-9;
    return SUBRETURN;
  }
#if 0 /* debugging print */
  printf("finished subroutine master \n");
  for(i=0; i<10; i++)
    printf("element %i, val %f, nvals=%i which=%i\n",i,mpivals[i],
	   mpinvals[i],all_vals[i]);
#endif
#else  /****** calculate everything without using MPI **************/

  for(job=0; job<N_ALL_BASES; job++){
    spectrum(all_bases+job,par,&eigen,&basis);
    for(i=0; i<all_bases[job].nval; i++)
      all_vals[job][i]=eigen.values[i];
  }
#endif
  /******** take results and calculate extrapolations etc. **************/

  /* Calculate the convergence in p-truncation */
  pextrap=pextrapolate(PNS,ps1,all_vals+(ps1-all_bases));

  /*  Calculate the transverse string tension.  The wht matrix
      must be reformatted into a form compatable with
      the FORTRAN routince WWEIGHTS.    */
  wweights(weights,&wth,nwind,wht0);
  PSEUDOINVERSE(weights,&wth,&nwind);
  for(j=0; j<wth; j++)wfitting[j]=0.0;
  for(i=0; i<nwind; i++){
    extrapolate(&vals,WNS,ws0[i],all_vals+(ws0[i]-all_bases),th,1,&pextrap);
    /* winding n=(2,0) eigenvalue saved */ 
    if(dot(wht0[i],wht0[i],HINDEX)==4)wvalue20=vals[0];
    for(j=0; j<wth; j++)
      wfitting[j] += vals[0]*weights[i+j*nwind];
  }
  
  /*  Calculate roundness of winding modes */
  for(i=0; i<FWIND; i++){
    extrapolate(&vals,WNS,ws[i],all_vals+(ws[i]-all_bases),th,1,&pextrap);
    wvalue[i]=vals[0];
  }
  
  /*  Calculate the longitudinal string tension  */
  lweights(weights,&lth,nlong,mlong0);
  PSEUDOINVERSE(weights,&lth,&nlong);
  for(j=0; j<lth; j++)lfitting[j]=0.0;
  for(i=0; i<nlong; i++){
    lextrapolate(&vals,LNS,ls0[i],all_vals+(ls0[i]-all_bases),
		  th,1,&pextrap,lkmax0);
    for(j=0; j<lth; j++)
      lfitting[j]+=vals[0]*weights[i+j*nlong];
    }
#if 0 /* debugging print */
  printf("  sigma=%f const=%f 1/L %f\n",
	   lfitting[0],lfitting[1],lfitting[2]);
#endif
  
  /*  Calculate roundness of heavy source potential.  */
  for(i=0; i<FLONG; i++){
    lextrapolate(&vals,LNS,ls[i],all_vals+(ls[i]-all_bases),
		  th,1,&pextrap,lkmax[i]);
    lvalue[i]=vals[0];
  }

  /*  Calculate the ordinary extrapolated spectra. */
  for(i=0; i<NSTILL; i++){
    extrapolate(&vals,NS,stillbases[i],all_vals+(stillbases[i]-all_bases),
		 th,NSTILLVALS,&pextrap);
    for(k=0; k<NSTILLVALS; k++)stillvals[i][k]=vals[k];
  }
  
  /*  Calculate the extrapolated spectrum in the P_perp = (c,0) direction */
  /*  Sort states before adding (and calculate 1 extra value) */
  for(l=0; l<NP1; l++)
    for(i=0; i<NPERP; i++){
      extrapolate(&vals,NS,s[l][i],all_vals+(s[l][i]-all_bases),
		   th,NPERPVALS+1,&pextrap);
      /* sort and add states: we assume NPERPVALS+1 eigenvalues */
      k=0;
      addterms(spectrumvals1[l+1][i],&k,NPERPVALS,vals,NPERPVALS+1);
#if 0 /* debugging print */
      printf("  l=%i, spectrum %i P_perp=(%f,%f)\n",l,i,p1a[l][0],p1a[l][1]);
#endif
    } 
  
  /*  Calculate the extrapolated spectrum in the P_perp = (c,c) direction */
  /*  Sort states before adding (and calculate 1 extra value) */
  for(l=0; l<NP2; l++)
    for(i=0; i<NPERP; i++){
      extrapolate(&vals,NS,t[l][i],all_vals+(t[l][i]-all_bases),
		   th,NPERPVALS+1,&pextrap);
      /* sort and add states: we assume NPERPVALS+1 eigenvalues */
      k=0;
      addterms(spectrumvals2[l+1][i],&k,NPERPVALS,vals,NPERPVALS+1);
#if 0 /* debugging print */
      printf("  l=%i spectrum %i P_perp=(%f,%f)\n",l,i,p2a[l][0],p2a[l][1]);
#endif
    }

  /* Free up memory used.  This is mainly for debugging */
  free_eigensystem(&eigen);
  free_spectrum();
  free_basis(&basis);
  for(i=0; i<N_ALL_BASES; i++)
    free(all_vals[i]);


  /*  Put ordinary spectrum values into nonzero P_perp array*/
  /*  Here we sort the states before adding.  */
  /*  The ordering of the states must be compatable with */
  /*  the ordering obtained when calculating nonzero p_perp */
  for(i=0; i<NPERP; i++){
    for(j=0, k=l=0; j<NSTILL; j++){
      if(stillbases[j]->charge!=s[0][i]->charge)continue;
      /* assume that the two spin 1 multiplets are degenerate */
      if(((stillbases[j]->multi==3 || stillbases[j]->multi==5 ||
	    stillbases[j]->multi==7) && s[0][i]->multi==15) ||
	  ((stillbases[j]->multi==4 || stillbases[j]->multi==6 ||
	    stillbases[j]->multi==7) && s[0][i]->multi==16))
	addterms(spectrumvals1[0][i],&k,NPERPVALS,stillvals[j],NSTILLVALS);
     if(((stillbases[j]->multi==3 || stillbases[j]->multi==6 ||
	    stillbases[j]->multi==7) && t[0][i]->multi==17) ||
	  ((stillbases[j]->multi==4 || stillbases[j]->multi==5 ||
	    stillbases[j]->multi==7) && t[0][i]->multi==18))
	addterms(spectrumvals2[0][i],&l,NPERPVALS,stillvals[j],NSTILLVALS);
    }
    if(k!=NPERPVALS || l!=NPERPVALS){
       fprintf(stderr,__FILE__":  Incorrect multiplet mapping from states "
	       "in ordinary\n"
	      "spectrum to those in nonzero P_perp, exiting\n");
       fprintf(stderr,"          k=%i, l=%i and NPERPVALS=%i\n",
	       k,l,NPERPVALS);
      *iflag=-12;
      return SUBRETURN;
    }	
  }
    
#if 1 /*debugging print of the entire spectrum */
  for(k=0; k<NPERPVALS; k++)
  for(i=0; i<NPERP; i++){
    for(j=0; j<NP1+1; j++)
      printf(" %g",spectrumvals1[j][i][k]);
    for(j=0; j<NP1+1; j++)
      printf(" %g",spectrumvals2[j][i][k]);
    printf("\n");
  } 
#endif
    
  /*  copy couplings into oldparams */
  for(i=0; i<NPARAMS; i++)oldpar[i]=par[i];
  
  /* There are several possible methods for setting
     the overall scale of the Hamiltonian  */
  
  switch(scalemethod){         
  case 0:          /*Find scale by long string tension */
    rescale=1.0/lfitting[0];
    break;
  case 1:    /*Find scale by looking only at lowest state */
    rescale=teper[0]/stillvals[0][0];
    break;
  case 2:    /* Find scale by minimizing chi^2 for c^2 measurement */
    sum=0.0; zz=0.0;
    for(i=0; i<NPERP; i++)
      for(j=0; j<FCSQUARED && j<NPERPVALS; j++){
	rescale=spectrumvals1[1][i][j]-spectrumvals1[0][i][j];
	sum += rescale/pow(c2err[i][j],2);
	zz += pow(rescale/c2err[i][j],2);
      }
    rescale=sum/zz*fdot(p1a[0],p1a[0],HINDEX)/wfitting[0];
    if(rescale>0.0)
      rescale=sqrt(rescale);
    else {
      rescale=0.0;
      *iflag = -1;
    }
    break;
  case 3:    /* Find scale by minimizing overall chi^2 */
    /* This does not include the consistancy of the L nonzero,
       n=1 datum because it is too complicated */
    /* initial guess is from Teper's lowest state */
    for(rescale=teper[0]/stillvals[0][0], j=0; ; j++){
      sum=0.0; zz=0.0;
      for(j=0; j<NPERP; j++)
	for(i=0; i<FCSQUARED && i<NPERPVALS; i++){
	  c2[0][j][i]=wfitting[0]*pow(rescale,2)*
	    (spectrumvals1[1][j][i]-spectrumvals1[0][j][i])/
	    fdot(p1a[0],p1a[0],HINDEX);
	  temp=2.0*c2[0][j][i]/rescale;
	  sum += temp*(c2[0][j][i]-1.0)/pow(c2err[j][i],2);
	  zz +=  pow(temp/c2err[j][i],2);
	  
	  c2[1][j][i]=wfitting[0]*pow(rescale,2)*
	    (spectrumvals2[1][j][i]-spectrumvals2[0][j][i])/
	    fdot(p2a[0],p2a[0],HINDEX);
	  temp=2.0*c2[1][j][i]/rescale;
	  sum += temp*(c2[1][j][i]-1.0)/pow(c2err[j][i],2);
	  zz +=  pow(temp/c2err[j][i],2);
	}
#if FTEPER==1
      sum += stillvals[0][0]*(stillvals[0][0]*rescale-teper[0])/
	pow(tepererr[0],2);
      zz += pow(stillvals[0][0]/tepererr[0],2);
#endif
#if FSCALE
      sum += lfitting[0]*(rescale*lfitting[0]-1.0)/pow(lerr0,2);
      zz += pow(lfitting[0]/lerr0,2);
#endif
      temp=rescale;
      /* Here I include a damping factor which will
	 slow convergence but hopefully stablize it. */
      rescale -= 0.25*sum/zz;
#if 0 /* debugging print */
      if(j<100)printf("rescale:  j=%i rescale=%.16g\n",j,rescale);
#endif
      /* converge as well as we can, assuming round-off errors */
      if(fabs(temp-rescale)<=(2.0+2.0*FCSQUARED)*ELEMENT_EPSILON)break;
      if(j==1000){
	 fprintf(stderr,__FILE__":  Scale determination not converged, "
		 "rescale=%g\n",rescale);
	*iflag=-2;
	break;
      }
    }
    break;
  default:
    fprintf(stderr,__FILE__":  in fcn, scale method=%i is wrong.\n",
	    scalemethod);
    *iflag=-14;
    rescale=teper[0]/stillvals[0][0];
  }
  /* If rescale is negative for some reason, 
     set it to something large and positive.        */
  if(rescale<0.0){
    rescale=50.0;
     fprintf(stderr,__FILE__":  rescale is negative, set recale=%f.\n",
	     rescale);
  }
  
  /*  
      Calculate the desired dispersion using the overall scale 
      determined above.
      Find error of speed of light squared
      with coupling g^2 N and "a" as determined above.
      One can weight the value for our value of sigma.
  */

  fitc=0;
  temp=wfitting[0]*pow(rescale,2)/fdot(p1a[0],p1a[0],HINDEX);
  for(i=0; i<NPERP; i++)
    for(j=0; j<NPERPVALS; j++)
    c2[0][i][j] = (spectrumvals1[1][i][j]-spectrumvals1[0][i][j])*temp;
  temp=wfitting[0]*pow(rescale,2.0)/fdot(p2a[0],p2a[0],HINDEX);
  for(i=0; i<NPERP; i++)
    for(j=0; j<NPERPVALS; j++)
    c2[1][i][j] = (spectrumvals2[1][i][j]-spectrumvals2[0][i][j])*temp;
  for(k=0; k<2; k++)
    for(i=0; i<NPERP; i++)
      for(j=0; j<FCSQUARED && j<NPERPVALS; j++)
      fvecc[fitc++] = (c2[k][i][j]-1.0)/c2err[i][j];
  
  /* If the 0^++ or the n=(2,0) is tachyonic, modify the 
     0^++ c^2 errors to reflect this */

  if(stillvals[0][0]<0.0){
     fprintf(stderr,__FILE__":  Lowest state gone tachyonc.\n");
      fvecc[0]=100.0*fabs(stillvals[0][0]);
  }
#if 0  /*  Add penalty if the winding modes become tachyonic.  */
  if(wvalue20<0.0){
     fprintf(stderr,__FILE__":  Winding n=(2,0) gone tachyonc.\n");
      fvecc[1]=100.0*fabs(wvalue20);
  }
#endif

  /* Add roundness of winding modes to fitting criterion. */
  for(i=0; i<FWIND; i++){
    wweights(weights,&wth,1,&wht[i]);
    calcwvalue[i]=fdot(weights,wfitting,wth);
    fvecc[fitc++]=(wvalue[i]-calcwvalue[i])/werr[i];
  }

  /* Add roundness of heavy source to fitting criterion.  
     Here, sum is the total
     distance in units of G^2 N. */
  for(i=0; i<FLONG; i++){
  calclvalue[i]=pow(mlong[i],2)+wfitting[0]*
    (element) dot(lht[i],lht[i],HINDEX)*pow(rescale,2);
  if(calclvalue[i]>0.0){
    calclvalue[i]=sqrt(calclvalue[i]);
    lweights(weights,&lth,1,&calclvalue[i]);
    calclvalue[i]=fdot(weights,lfitting,lth);
    fvecc[fitc++]=(lvalue[i]-calclvalue[i])/lerr[i];
  } else {
     fprintf(stderr,__FILE__":  value for calclvalue[%i]=%g bad, mlong=%g wfitting[0]=%g\n"
	    "  dot=%g rescale=%g\n",i,calclvalue[i],mlong[i],
	    wfitting[0],(element) dot(lht[i],lht[i],HINDEX),rescale);
    fvecc[fitc++]=100.0*fabs(calclvalue[i])+10.0;
  }
#if 1==0 /* debugging print  */
  printf("fitc=%i fvecc=%f calclvalue[%i]=%f lerr=%f lvalue=%f\n",
	 fitc,fvecc[fitc],calclvalue[i],lerr[i],lvalue[i]);
#endif
  }

  /*  Consistancy of longitudinal string tension with rescale */
#if FSCALE
  fvecc[fitc++]=(rescale*lfitting[0]-1.0)/lerr0;
#endif

  /* Add parity doublets to fitting criterion using fractional
     difference.   */

  for(k=0; k<FPARITY; k++){
    for(i=0; i<NSTILL; i++)
      if(stillbases[i]->multi==whichparity[k][0].multi && 
	 stillbases[i]->charge==whichparity[k][0].charge)break;
    for(j=0; j<NSTILL; j++)
      if(stillbases[j]->multi==whichparity[k][1].multi && 
	 stillbases[j]->charge==whichparity[k][1].charge)break;
    if(i<NSTILL && j<NSTILL)
      fvecc[fitc++]=(stillvals[i][whichparity[k][0].state]-
		     stillvals[j][whichparity[k][1].state])/
	(0.5*(stillvals[i][whichparity[k][0].state]+
	      stillvals[j][whichparity[k][1].state])*parityerr[k]);
    else {
      fprintf(stderr,__FILE__":  Incorrect multiplets in ordinary spectra, "
	      " exiting\n");
      *iflag=-15;
      return SUBRETURN;
    }
  }

  /* Add fit to teper's value for the lowest state
     we assume that the lowest state is in the first sector
     of stillvals.   */
#if FTEPER==1
    fvecc[fitc++]=(rescale*stillvals[0][0]-teper[0])/tepererr[0];
#endif

    /* Demand that the lattice spacing be a specific value. */
#if FEXTRA==1
    fvecc[fitc++]=(1.0+2.0*sqrt(par[0]/par[1])-wfitting[0]*rescale)/
      extraerr[0];
#endif
    
  /* after including all fit criteria, make sure we 
     have actually filled the fvecc vector. */
  if(fitc!=*m){
    fprintf(stderr,__FILE__":  Number of fit criteria wrong.\n");
    *iflag=-75;
    return SUBRETURN;
  }
    
    /* check if this is the best chi^2 so far,
       if so, copy values into the structure best.
       bestpar makes sure that the unfitted variables are
       the same; that is the previous best values are
       from the same fit. */
  sum=fdot(fvecc,fvecc,*m);
  for(i=0; i<NPARAMS && fabs(bestpar[i]-par[i])<=2*ELEMENT_EPSILON; i++);
  if(*iflag>=0 && (i<NPARAMS || sum<bestchi2)){
    bestchi2=sum;
    for(i=0; i<NPARAMS; i++)bestpar[i]=par[i];
    best.pextrap=pextrap; 
    for(i=0; i<FWIND; i++){
      best.wvalue[i]=wvalue[i];
      best.calcwvalue[i]=calcwvalue[i];
    }
    for(i=0; i<FLONG; i++){
      best.lvalue[i]=lvalue[i];
      best.calclvalue[i]=calclvalue[i];
    }
    best.rescale=rescale;
    for(k=0; k<2; k++)
      for(i=0; i<NPERP; i++)
	for(j=0; j<NPERPVALS; j++)
	  best.c2[k][i][j]=c2[k][i][j];
    for(i=0; i<NPERP; i++)
      for(k=0; k<NPERPVALS; k++){
	for(j=0; j<NP1+1; j++)
	  best.spectrumvals1[j][i][k]=spectrumvals1[j][i][k];
	for(j=0; j<NP2+1; j++)
	  best.spectrumvals2[j][i][k]=spectrumvals2[j][i][k];
      }
    for(i=0; i<NSTILL; i++)
      for(j=0; j<NSTILLVALS; j++)
	best.stillvals[i][j]=stillvals[i][j];
    for(i=0; i<wth; i++)best.wfitting[i]=wfitting[i];
    for(i=0; i<lth; i++)best.lfitting[i]=lfitting[i];
  }

#if PRINTIT
  printf("  fvecc values =");
  for(i=0, k=3; i<*m; i++){
    printf(" %g",fvecc[i]);
    if((k++)%6==0 && i+1<*m)printf("\n  ");
  }
  printf("\n  Chi^2=%g in %f CPU seconds.\n\n",
	 fdot(fvecc,fvecc,*m),
	 (float)(clock()-tif)/(float) CLOCKS_PER_SEC);
#endif
#if PRINTIT 
  fflush(stdout);
#endif

  free(vals); /* vals is realloc()ed by extrapolation routines */

  return SUBRETURN; /* normal return from subroutine */

}

#endif /* HINDEX==2 */