/************************************
 Current revision $Revision$
 On branch $Name$
 Latest change $Date$ by $Author$
*************************************/
#define ORDER 1    
//#define ORDER 4   

int no_tps = ORDER; // arbitrary TPSA order is defined locally

extern bool freq_map;

#if HAVE_CONFIG_H
#include <config.h>
#endif

#if MPI_EXEC
  //library of parallel computation,must be before "stdio.h"
  #include <mpi.h>
#endif

  #include "tracy_lib.h"


//***************************************************************************************
//
//  MAIN CODE
//
//****************************************************************************************
int main(int argc, char *argv[]) {
  
  /* for time handling */
  uint32_t start, stop;
  
  /*set the random value for random generator*/
  const long seed = 1121; //the default random seed number
  iniranf(seed); //initialize the seed
  setrancut(2.0); //default value of the normal cut for the normal distribution
  // turn on globval.Cavity_on and globval.radiation to get proper synchro radiation damping
  // IDs accounted too if: wiggler model and symplectic integrator (method = 1)
  globval.H_exact = false;

  /* parameters to read the user input script .prm */
  long i=0L; //initialize the for loop to read command string
  long j=0L; //initialize the for loop to read the files from parallel computing
  char CommandStr[max_str];
  double nux = 0.0, nuy = 0.0, ksix = 0.0, ksiy = 0.0;
  bool chroma=true;
  double dP = 0.0;
  long lastpos = -1L;
  char str1[S_SIZE];
  
  // paramters to read the command line in user input script
  long CommandNo = -1L;  //the number of commands, since this value is also
                         // the index of array UserCommandFlag[], so this
			 // value is always less than 1 in the really case. 
  const int NCOMMAND = 500;
  UserCommand UserCommandFlag[NCOMMAND];
 
#if MPI_EXEC
  //Initialize parallel computation
    MPI_Init(&argc, &argv);
#endif
  /************************************************************************
   read in files and flags
   *************************************************************************/

   
  if (argc > 1) {
    read_script(argv[1], true,CommandNo, UserCommandFlag);
  } else {
    fprintf(stdout, "Not enough parameters\nSyntax is program parameter file\n");
    exit_(1);
  }
     
  /* get the family index of the special elements, prepare for printglob()*/
  if(strcmp(bpm_name,"")==0)
    globval.bpm = 0;
  else
    globval.bpm = ElemIndex(bpm_name);
  if(strcmp(skew_quad_name,"")==0)
    globval.qt = 0;  
  else
    globval.qt = ElemIndex(skew_quad_name);
  if(strcmp(hcorr_name,"")==0)
    globval.hcorr = 0;
  else
    globval.hcorr = ElemIndex(hcorr_name);
  if(strcmp(vcorr_name,"")==0)
    globval.vcorr = 0;
  else
    globval.vcorr = ElemIndex(vcorr_name);
  if(strcmp(gs_name,"")==0)
    globval.gs = 0;
  else
    globval.gs = ElemIndex(gs_name);
  if(strcmp(ge_name,"")==0)
    globval.ge = 0;
  else
    globval.ge = ElemIndex(ge_name);
  
//  globval.g = ElemIndex("g");  /* get family index of  girder*/
   
  /* print the summary of the element in lattice */  
  printglob();   
  /************************************************************************
    print files, very important file for debug
   *************************************************************************/
  
  //print flat file with all the design values of the lattice,
  prtmfile("flat_file.dat");
  // print location, twiss parameters and close orbit at all elements position to a file
  getcod(dP, lastpos);
  prt_cod("cod.out", globval.bpm, true);

  //get_matching_params_scl(); // get tunes and beta functions at entrance
  // compute up to 3rd order momentum compact factor
  get_alphac2(); 
  //cout << endl << "computing tune shifts" << endl;
  //dnu_dA(10e-3, 5e-3, 0.002);
  //get_ksi2(0.0); // this gets the chromas and writes them into chrom2.out
  
    
  
   
  
/********************************************************************* 
                            Flag factory 
*********************************************************************/
   
  for(i=0; i<=CommandNo; i++){//read user defined command by sequence
   //assign user defined command
    strcpy(CommandStr,UserCommandFlag[i].CommandStr); 
    
       //turn on flag for quadrupole fringe field
   if(strcmp(CommandStr,"QuadFringeOnFlag") == 0) {
      globval.quad_fringe = true;
      cout << "\n";
      cout << "globval.quad_fringe is " << globval.quad_fringe << "\n";
   }
   
         //deactive quadrupole fringe field
   else if(strcmp(CommandStr,"QuadFringeOffFlag") == 0) {
      globval.quad_fringe = false;
      cout << "\n";
      cout << "globval.quad_fringe is " << globval.quad_fringe << "\n";
   }
   
  //set RF voltage  
  else if(strcmp(CommandStr,"RFvoltageFlag") == 0) {
    printf("\nSetting RF voltage:\n");
    printf("    Old RF voltage is: %lf [MV]\n", get_RFVoltage(globval.cav)
        / 1e6);
    set_RFVoltage(globval.cav, UserCommandFlag[i].RFvolt);
    printf("    New RF voltage is: %lf [MV]\n", get_RFVoltage(globval.cav)
        / 1e6);
  }

  // Chamber factory
  else if(strcmp(CommandStr,"ReadChamberFlag") == 0) {
    ReadCh(UserCommandFlag[i].chamber_file); /* read vacuum chamber from a file "Apertures.dat" , soleil version*/
    PrintCh(); // print chamber into chamber.out
  }
  
  // read the misalignment errors to the elements
  //  Based on the function error_and_correction() in nsls-ii_lib_templ.h
  else if (strcmp(CommandStr,"ReadaefileFlag") == 0) {
   // load misalignments
    cout << "Read misalignment errors from file: " << UserCommandFlag[i].ae_file << endl;
    LoadAlignTol(UserCommandFlag[i].ae_file, true, 1.0, true, 0);
    Ring_GetTwiss(true, 0.0);
  }
    //do COD correction using SVD method.
  else if(strcmp(CommandStr,"CODCorrectFlag") == 0) {
    cout << "Begin Closed Orbit correction: " << endl;
    
    if(n_orbit < 1){
    cout << "interation number n_obit should NOT small than 1!!!" << endl;
    exit_(1);
    }
    
    CODCorrect(hcorr_file,vcorr_file, n_orbit,nwh,nwv);
    Ring_GetTwiss(true, 0.0);
    prt_cod("summary_miserr_codcorr.out", globval.bpm, true);
  }
   // set the field error into the lattice 
  // The corresponding field error is replaced by the new value.
  // This feature is generic, works for all lattices
  else if (strcmp(CommandStr,"ReadfefileFlag") == 0) {
    fprintf(stdout,"\n Read multipole field errors from file:      %s \n", UserCommandFlag[i].fe_file);
    LoadFieldErrs(UserCommandFlag[i].fe_file, true, 1.0, true, 1);
	
    Ring_GetTwiss(true, 0.0);
    prtmfile("flat_file_fefile.dat");
  }

  // read multipole errors from a file; specific for soleil lattice
  else if(strcmp(CommandStr,"ReadMultipoleFlag") == 0) {
    fprintf(stdout,"\n Read Multipoles field error for SOLEIL lattice with thick sextupoles, from file:\n");
    fprintf(stdout,"        %s\n",multipole_file);
    ReadFieldErr(multipole_file);
    //first print the full lattice with error as a flat file
    prtmfile("flat_file_errmultipole.dat"); // writes flat file with all element errors  /* very important file for debug*/
    Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
    printglob();
  }
   // read the sources of coupling from a file; for soleil lattice
  else if(strcmp(CommandStr,"ReadVirtualSkewquadFlag") == 0) {
    fprintf(stdout,"\n Read virtual skew quadrupole setting of SOLEIL lattice from file:\n");
    fprintf(stdout,"    %s \n",virtualskewquad_file);
    
    ReadVirtualSkewQuad(virtualskewquad_file);
    //first print the full lattice with sources of coupling as a flat file
    prtmfile("flat_file_skewquad.dat");  /* very important file for debug*/
    //get the coupling
    Coupling_Edwards_Teng(); 
  }

     //print the twiss paramaters in a file defined by the name
   else if(strcmp(CommandStr,"PrintTwissFlag") == 0) {
      cout << "\n";
      cout << "print the twiss parameters to file: "
           << UserCommandFlag[i]._PrintTwiss_twiss_file << "\n";
	   
      printlatt(UserCommandFlag[i]._PrintTwiss_twiss_file);  
   }
   //print the close orbit
   else if(strcmp(CommandStr,"PrintCODFlag") == 0) {
      cout << "\n";
      cout << "print the close orbit to file: "
           << UserCommandFlag[i]._PrintCOD_cod_file << "\n";
      getcod(dP, lastpos);
      prt_cod(UserCommandFlag[i]._PrintCOD_cod_file, globval.bpm, true);
   }
       //print coordinates at lattice elements obtained by tracking
   else if(strcmp(CommandStr,"PrintTrackFlag") == 0) {
      cout << "\n";
      cout << "print the tracked coordinates to file: "
           << UserCommandFlag[i]._PrintTrack_track_file << "\n";
      	
      PrintTrack(UserCommandFlag[i]._PrintTrack_track_file,
                 UserCommandFlag[i]._PrintTrack_x, UserCommandFlag[i]._PrintTrack_px, 
                 UserCommandFlag[i]._PrintTrack_y, UserCommandFlag[i]._PrintTrack_py, 
                 UserCommandFlag[i]._PrintTrack_delta,UserCommandFlag[i]._PrintTrack_ctau,
		 UserCommandFlag[i]._PrintTrack_nmax);  
   }
   
    //print the girder
  // else if(strcmp(CommandStr,"PrintGirderFlag") == 0) {
   //   cout << "\n";
   //   cout << "print the information of girder to file: "<< girder_file << "\n";
      
      
     // getcod(dP, lastpos);
     // prt_cod(cod_file, globval.bpm, true);
  // }
   

  
  // compute tunes by tracking (should be the same as by tps)
  else if (strcmp(CommandStr,"TuneTracFlag") == 0) {
    GetTuneTrac(1026L, 0.0, &nux, &nuy);
    fprintf(stdout, "From tracking: nux = % f nuz = % f \n", nux, nuy);
    }

  // compute chromaticities by tracking (should be the same as by DA)
  else if(strcmp(CommandStr,"ChromTracFlag") == 0) {
    start = stampstart();
    GetChromTrac(2L, 1026L, 1e-5, &ksix, &ksiy);
    stop = stampstop(start);
    fprintf(stdout, "From tracking: ksix= % f ksiz= % f \n", ksix, ksiy);
    }
    

  //generic function, to fit tunes using 1 family of quadrupoles
 // if (FitTuneFlag == true) {
  else if(strcmp(CommandStr,"FitTuneFlag") == 0) {
    double qf_bn = 0.0, qd_bn = 0.0;
    double qf_bn_fit = 0.0, qd_bn_fit = 0.0;

    /* quadrupole field strength before fitting*/
    qf_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf), 1)].Elem.M->PB[HOMmax + 2];
    qd_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd), 1)].Elem.M->PB[HOMmax + 2];

    /* fitting tunes*/
    fprintf(stdout,
        "\n Fitting tunes: %s %s, targetnux = %f, targetnuz = %f \n", UserCommandFlag[i].qf, 
	UserCommandFlag[i].qd,UserCommandFlag[i].targetnux, UserCommandFlag[i].targetnuz);
    FitTune(ElemIndex(UserCommandFlag[i].qf), ElemIndex(UserCommandFlag[i].qd), 
            UserCommandFlag[i].targetnux, UserCommandFlag[i].targetnuz);

    /* integrated field strength after fitting*/
    qf_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf), 1)].Elem.M->PB[HOMmax + 2];
    qd_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd), 1)].Elem.M->PB[HOMmax + 2];
    /* print out the quadrupole strengths before and after the fitting*/
    printf("Before the fitting, the quadrupole field strengths are: \n");
    printf(" %s = %f,    %s = %f\n", UserCommandFlag[i].qf, qf_bn, UserCommandFlag[i].qd, qd_bn);
    printf("After the fitting, the quadrupole field strengths are: \n");
    printf(" %s = %f,    %s = %f\n", UserCommandFlag[i].qf, qf_bn_fit, UserCommandFlag[i].qd, qd_bn_fit);

    /* Compute and get Twiss parameters */
    Ring_GetTwiss(chroma = true, 0.0);
    printglob(); /* print parameter list */
  }

  /* specific for soleil ring in which the quadrupole is cut into 2 parts*/
  //if (FitTune4Flag == true) {
  else if(strcmp(CommandStr,"FitTune4Flag") == 0) {
    double qf1_bn = 0.0, qf2_bn = 0.0, qd1_bn = 0.0, qd2_bn = 0.0;
    double qf1_bn_fit = 0.0, qf2_bn_fit = 0.0, qd1_bn_fit = 0.0, qd2_bn_fit =
        0.0;

    /* quadrupole field strength before fitting*/
    qf1_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf1), 1)].Elem.M->PB[HOMmax + 2];
    qf2_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf2), 1)].Elem.M->PB[HOMmax + 2];
    qd1_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd1), 1)].Elem.M->PB[HOMmax + 2];
    qd2_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd2), 1)].Elem.M->PB[HOMmax + 2];

    /* fitting tunes*/
    fprintf(
        stdout,
        "\n Fitting tunes for Soleil ring: %s %s %s %s, targetnux = %f, targetnuz = %f \n",
        UserCommandFlag[i].qf1, UserCommandFlag[i].qf2, UserCommandFlag[i].qd1, UserCommandFlag[i].qd2, 
	UserCommandFlag[i].targetnux, UserCommandFlag[i].targetnuz);
    FitTune4(ElemIndex(UserCommandFlag[i].qf1), ElemIndex(UserCommandFlag[i].qf2),
             ElemIndex(UserCommandFlag[i].qd1), ElemIndex(UserCommandFlag[i].qd2),
             UserCommandFlag[i].targetnux, UserCommandFlag[i].targetnuz);

    /* integrated field strength after fitting*/
    qf1_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf1), 1)].Elem.M->PB[HOMmax + 2];
    qf2_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf2), 1)].Elem.M->PB[HOMmax + 2];
    qd1_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd1), 1)].Elem.M->PB[HOMmax + 2];
    qd2_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd2), 1)].Elem.M->PB[HOMmax + 2];
    /* print out the quadrupole strengths before and after the fitting*/
    printf("Before the fitting, the quadrupole field strengths are: \n");
    printf("    %s = %f,    %s = %f,    %s = %f,    %s = %f\n", UserCommandFlag[i].qf1, qf1_bn,
        UserCommandFlag[i].qf2, qf2_bn, UserCommandFlag[i].qd1, qd1_bn, UserCommandFlag[i].qd2, qd2_bn);
    printf("After the fitting, the quadrupole field strengths are: \n");
    printf("    %s = %f,    %s = %f,    %s = %f,    %s = %f\n", UserCommandFlag[i].qf1,
        qf1_bn_fit, UserCommandFlag[i].qf2, qf2_bn_fit, UserCommandFlag[i].qd1, qd1_bn_fit, 
	UserCommandFlag[i].qd2, qd2_bn_fit);

    /* Compute and get Twiss parameters */
    Ring_GetTwiss(chroma = true, 0.0);
    printglob(); /* print parameter list */
  }

  /* fit the chromaticities*/
  else if(strcmp(CommandStr,"FitChromFlag") == 0) {
  
    double sxm1_bn = 0.0, sxm2_bn = 0.0;
    double sxm1_bn_fit = 0.0, sxm2_bn_fit = 0.0;
    sxm1_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].sxm1), 1)].Elem.M->PB[HOMmax + 3];
    sxm2_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].sxm2), 1)].Elem.M->PB[HOMmax + 3];

    fprintf(stdout,"\n Fitting chromaticities: %s %s, targetksix = %f,  targetksiz = %f\n",
        UserCommandFlag[i].sxm1, UserCommandFlag[i].sxm2, UserCommandFlag[i].targetksix,
	 UserCommandFlag[i].targetksiz);
	 
    FitChrom(ElemIndex(UserCommandFlag[i].sxm1), ElemIndex(UserCommandFlag[i].sxm2), 
             UserCommandFlag[i].targetksix, UserCommandFlag[i].targetksiz);

    sxm1_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].sxm1), 1)].Elem.M->PB[HOMmax + 3];
    sxm2_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].sxm2), 1)].Elem.M->PB[HOMmax + 3];
    /* print out the sextupole strengths before and after the fitting*/
    printf("Before the fitting, the sextupole field strengths are \n");
    printf("    %s = %f,    %s = %f\n", UserCommandFlag[i].sxm1, sxm1_bn, UserCommandFlag[i].sxm2, sxm2_bn);
    printf("After the fitting, the sextupole field strengths are: \n");
    printf("    %s = %f,    %s = %f\n", UserCommandFlag[i].sxm1, sxm1_bn_fit, UserCommandFlag[i].sxm2, sxm2_bn_fit);

    Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
    printglob(); /* print parameter list */
  }

  // coupling calculation
  else if(strcmp(CommandStr,"CouplingFlag") == 0) {
    Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
    printlatt("linlat_coupling.out"); /* dump linear lattice functions into "linlat.dat" */
   // GetEmittance(ElemIndex("cav"), true);  //only for test
    Coupling_Edwards_Teng();
   // Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
  //  printglob(); /* print parameter list */
  }

  // add coupling by random rotating of the full quadrupole magnets
  //if (ErrorCouplingFlag == true) {
  else if(strcmp(CommandStr,"ErrorCouplingFlag") == 0) {
    prtmfile("flat_file.dat"); //print the elements without rotation errors
    SetErr(UserCommandFlag[i].err_seed, UserCommandFlag[i].err_rms);
    prtmfile("flat_file_errcoupling_full.dat"); //print the elements with rotation errors
    Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
    printlatt("linlat_errcoupling.out"); /* dump linear lattice functions into "linlat.dat" */
    Coupling_Edwards_Teng();
   // GetEmittance(ElemIndex("cav"), true); //only for test
    Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
  //  printlatt();
    printglob(); /* print parameter list */
  }
  
//  add coupling by random rotating of the half quadrupole magnets, delicated for soleil
  else if(strcmp(CommandStr,"ErrorCoupling2Flag") == 0) {
    prtmfile("flat_file.dat"); //print the elements without rotation errors
    SetErr2(UserCommandFlag[i].err_seed, UserCommandFlag[i].err_rms);
    prtmfile("flat_file_errcoupling_half.dat"); //print the elements with rotation errors
    Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
    printlatt("linlat_errcoupling2.out"); /* dump linear lattice functions into "linlat.dat" */
    Coupling_Edwards_Teng();
   // GetEmittance(ElemIndex("cav"), true);
   Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
    printglob(); /* print parameter list */
  }


  
  

  /******************************************************************************************/
  // COMPUTATION PART after setting the model
  /******************************************************************************************/
  // computes TuneShift with amplitudes
  else if(strcmp(CommandStr,"AmplitudeTuneShiftFlag") == 0) {
      TunesShiftWithAmplitude(UserCommandFlag[i]._AmplitudeTuneShift_nudx_file,
		              UserCommandFlag[i]._AmplitudeTuneShift_nudz_file,
                              UserCommandFlag[i]._AmplitudeTuneShift_nxpoint,
                              UserCommandFlag[i]._AmplitudeTuneShift_nypoint,
                              UserCommandFlag[i]._AmplitudeTuneShift_nturn, 
			      UserCommandFlag[i]._AmplitudeTuneShift_xmax,
                              UserCommandFlag[i]._AmplitudeTuneShift_ymax, 
			      UserCommandFlag[i]._AmplitudeTuneShift_delta);
    } 
  // compute tuneshift with energy
 else if(strcmp(CommandStr,"EnergyTuneShiftFlag") == 0) {
     TunesShiftWithEnergy(UserCommandFlag[i]._EnergyTuneShift_nudp_file,
                          UserCommandFlag[i]._EnergyTuneShift_npoint, 
                          UserCommandFlag[i]._EnergyTuneShift_nturn,
                          UserCommandFlag[i]._EnergyTuneShift_deltamax);
   }

  // Computes FMA
  else if(strcmp(CommandStr,"FmapFlag") == 0) {
    printf("\n begin Fmap calculation for on momentum particles: \n");
    
   #if MPI_EXEC

    //initialization for parallel computing
    int myid=0, numprocs=0;
    int namelen=0;
    char processor_name[MPI_MAX_PROCESSOR_NAME]="";
    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
    MPI_Get_processor_name(processor_name, &namelen);
    printf("\n Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name);  
    
    printf("\n begin Fmap calculation for on momentum particles: \n");
    //Each core or process, characterized by myid, will track different region of fmap 
    fmap_p(UserCommandFlag[i]._FmapFlag_fmap_file,
	   UserCommandFlag[i]._FmapFlag_nxpoint, 
	   UserCommandFlag[i]._FmapFlag_nypoint, 
	   UserCommandFlag[i]._FmapFlag_nturn, 
	   UserCommandFlag[i]._FmapFlag_xmax,
	   UserCommandFlag[i]._FmapFlag_ymax, 
	   UserCommandFlag[i]._FmapFlag_delta, 
	   UserCommandFlag[i]._FmapFlag_diffusion,
       UserCommandFlag[i]._FmapFlag_printloss,
       numprocs,myid);

    //Synchronize all cores, till all cores finish fmap of different region,then all cores will proceed.
    MPI_Barrier(MPI_COMM_WORLD);
   
    //Collecting data from all files generated by cores, then write to fmap_p.out
    if(myid==0)
      {
	//  system("cp 0fmap.out fmap_p.out");
	
	FILE *outf;
	if ((outf = fopen(UserCommandFlag[i]._FmapFlag_fmap_file, "w")) == NULL)
	  {
	    fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapFlag_fmap_file);
	    exit_(1);
	  }

	FILE *fp;
	char buffer[81];
	char FmapFile[50];
	for(int j=0; j<numprocs; j++)
	  {
	    FmapFile[0]='\0';
	    sprintf(FmapFile,"%d",j);
	    strcat(FmapFile,UserCommandFlag[i]._FmapFlag_fmap_file);
      
	    if((fp = fopen(FmapFile, "r")) == NULL)  
	      {
		fprintf(stdout, "%s: error while opening file.\n", FmapFile);
		exit_(1);
	      }

	    while(fgets(buffer, 80, fp) != NULL)
	      {
		fputs(buffer, outf);
		buffer[0]='\0';
	      }
	    fclose(fp);
	  }
	fclose(outf);
      }
    MPI_Barrier(MPI_COMM_WORLD); 
    printf("Process %d finished for on momentum fmap.\n", myid);

 #else
      fmap(UserCommandFlag[i]._FmapFlag_fmap_file,
	   UserCommandFlag[i]._FmapFlag_nxpoint, 
	   UserCommandFlag[i]._FmapFlag_nypoint, 
	   UserCommandFlag[i]._FmapFlag_nturn, 
	   UserCommandFlag[i]._FmapFlag_xmax,
	   UserCommandFlag[i]._FmapFlag_ymax, 
	   UserCommandFlag[i]._FmapFlag_delta, 
	   UserCommandFlag[i]._FmapFlag_diffusion,
	   UserCommandFlag[i]._FmapFlag_printloss);
#endif
  }

  // Compute FMA dp
  else if(strcmp(CommandStr,"FmapdpFlag") == 0) {
    printf("\n begin Fmap calculation for off momentum particles: \n");
    
   #if MPI_EXEC

    //initialization for parallel computing
    int myid=0, numprocs=0;
    int namelen=0;
    char processor_name[MPI_MAX_PROCESSOR_NAME]="";
    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
    MPI_Get_processor_name(processor_name, &namelen);
    printf("Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name);  

    printf("\n begin Fmap calculation for off momentum particles: \n");
    fmapdp_p(UserCommandFlag[i]._FmapdpFlag_fmapdp_file,
	     UserCommandFlag[i]._FmapdpFlag_nxpoint, 
	     UserCommandFlag[i]._FmapdpFlag_nepoint, 
	     UserCommandFlag[i]._FmapdpFlag_nturn,
	     UserCommandFlag[i]._FmapdpFlag_xmax, 
	     UserCommandFlag[i]._FmapdpFlag_emax, 
	     UserCommandFlag[i]._FmapdpFlag_z,
	     UserCommandFlag[i]._FmapdpFlag_diffusion,
         UserCommandFlag[i]._FmapFlag_printloss
         numprocs,myid);

    //Synchronize all cores, till all cores finish fmap of different region,then all cores will proceed.
    MPI_Barrier(MPI_COMM_WORLD);
   
    //Collecting data from all files generated by cores, then write to fmap_p.out
    if(myid==0)
     {
       //  system("cp 0fmap.out fmap_p.out");

       FILE *outf;
       if ((outf = fopen(UserCommandFlag[i]._FmapdpFlag_fmapdp_file, "w")) == NULL)
	 {
	   fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapdpFlag_fmapdp_file);
	   exit_(1);
	 }

       FILE *fp;
       char buffer[81];
       char FmapFile[50];
       for(int j=0; j<numprocs; j++)
	 {
	   FmapFile[0]='\0';
	   sprintf(FmapFile,"%d",j);
	   strcat(FmapFile,UserCommandFlag[i]._FmapdpFlag_fmapdp_file);
      
	   if((fp = fopen(FmapFile, "r")) == NULL)  
	     {
	       fprintf(stdout, "%s: error while opening file.\n", FmapFile);
	       exit_(1);
	     }

	   while(fgets(buffer, 80, fp) != NULL)
	     {
	       fputs(buffer, outf);
	       buffer[0]='\0';
	     }
	   fclose(fp);
	 }
       fclose(outf);
     }
    MPI_Barrier(MPI_COMM_WORLD); 
    printf("Process %d finished off momentum fmap.\n", myid);

 #else
       fmapdp(UserCommandFlag[i]._FmapdpFlag_fmapdp_file,
	      UserCommandFlag[i]._FmapdpFlag_nxpoint, 
	      UserCommandFlag[i]._FmapdpFlag_nepoint, 
	      UserCommandFlag[i]._FmapdpFlag_nturn,
	      UserCommandFlag[i]._FmapdpFlag_xmax, 
	      UserCommandFlag[i]._FmapdpFlag_emax, 
	      UserCommandFlag[i]._FmapdpFlag_z,
	      UserCommandFlag[i]._FmapdpFlag_diffusion,
	       UserCommandFlag[i]._FmapdpFlag_printloss);
 #endif
  }

//  // if (CodeComparaisonFlag) {
//   else if(strcmp(CommandStr,"CodeComparaisonFlag") == 0) {
//     fmap(200, 100, 1026, -32e-3, 7e-3, 0.0, true);
//   }

  // MOMENTUM ACCEPTANCE
  else if(strcmp(CommandStr,"MomentumAccFlag") == 0) {
    bool cavityflag, radiationflag;
     
    /* record the initial values*/
    cavityflag = globval.Cavity_on;
    radiationflag = globval.radiation; 
    
    if (strcmp(UserCommandFlag[i]._MomentumAccFlag_TrackDim,"6D") == 0) {
      globval.Cavity_on = true;
      globval.radiation = true;
    }else if (strcmp(UserCommandFlag[i]._MomentumAccFlag_TrackDim,"4D") == 0) {
      globval.Cavity_on = false;
      globval.radiation = false;
    } else {
      printf("MomentumAccFlag: Error!!! DimTrack must be '4D' or '6D'\n");
      exit_(1);
    };

 /* calculate momentum acceptance*/
    printf("\n Calculate momentum acceptance: \n");

 #if MPI_EXEC
    
    /* calculate momentum acceptance*/
    //initialization for parallel computing
    int myid=0, numprocs=0;
    int namelen=0;
    char processor_name[MPI_MAX_PROCESSOR_NAME]="";
    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
    MPI_Get_processor_name(processor_name, &namelen);
    printf("Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name);  

    printf("\n Calculate momentum acceptance: \n");
    MomentumAcceptance_p(UserCommandFlag[i]._MomentumAccFlag_momacc_file,
                       UserCommandFlag[i]._MomentumAccFlag_istart, 
                       UserCommandFlag[i]._MomentumAccFlag_istop,
                       UserCommandFlag[i]._MomentumAccFlag_deltaminp, 
		       UserCommandFlag[i]._MomentumAccFlag_deltamaxp,
                       UserCommandFlag[i]._MomentumAccFlag_nstepp, 
		       UserCommandFlag[i]._MomentumAccFlag_deltaminn,
                       UserCommandFlag[i]._MomentumAccFlag_deltamaxn, 
		       UserCommandFlag[i]._MomentumAccFlag_nstepn,
		       UserCommandFlag[i]._MomentumAccFlag_nturn,
		       UserCommandFlag[i]._MomentumAccFlag_zmax,
                       numprocs,myid);
  //merge the files 
    //Synchronize all cores, till finish cal of different region, 
    //then files from the cores will be processed.
    MPI_Barrier(MPI_COMM_WORLD);

    //Collect data from all files generated by different cores, then write to user defined file
    if(myid==0)
      {
	//merge the phase.out from different cpus
	FILE *outf1;  //phase.out, for debugging
	FILE *fp1;
	char buffer1[81];
	char PhaseFile[50];
	if ((outf1 = fopen("phase_p.out", "w")) == NULL)
	  {
	    fprintf(stdout, "psoltracy: error while opening file %s\n", "phase_p.out");
	    exit_(1);
	  }

	for(int j=0; j<numprocs; j++)
	  {
	    PhaseFile[0]='\0';
	    sprintf(PhaseFile,"%d",j);
	    strcat(PhaseFile,"phase.out");
	    
	    if((fp1 = fopen(PhaseFile, "r")) == NULL){
		fprintf(stdout, "%s: error while opening file.\n", PhaseFile);
		exit_(1);
	      }

	    while(fgets(buffer1, 80, fp1) != NULL)
	      {
		fputs(buffer1, outf1);
		buffer1[0]='\0';
	      }
	    fclose(fp1);
	  }

	//collect data for the momentum acceptance
	FILE *outf2; //momentum acceptance
	FILE *fp2;
	char buffer2[81];
	char MAFile[50];
	if ((outf2 = fopen(UserCommandFlag[i]._MomentumAccFlag_momacc_file, "w")) == NULL)
	  {
	    fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._MomentumAccFlag_momacc_file);
	    exit_(1);
	  }

	//Collect data in Positive region
	for(int j=0; j<numprocs; j++)
	  {
	    MAFile[0]='\0';
	    sprintf(MAFile,"%d",j);
	    strcat(MAFile,UserCommandFlag[i]._MomentumAccFlag_momacc_file);

	    if((fp2 = fopen(MAFile, "r")) == NULL)  
	      {
		fprintf(stdout, "%s: error while opening file.\n", MAFile);
		exit_(1);
	      }
	 
	    while(fgets(buffer2, 80, fp2) != NULL)
	      {
		if(strstr(buffer2,"Negative")!=0) break;

		fputs(buffer2, outf2);
		buffer2[0]='\0'; //reset buffer to NULL
	      }
	    buffer2[0]='\0';
	    fclose(fp2);
	  }

	fprintf(outf2,"\n"); // A void line

	//Collect data in Nagative region
	for(int j=0; j<numprocs; j++)
	  {
	    MAFile[0]='\0';
	    sprintf(MAFile,"%d",j);
	    strcat(MAFile,UserCommandFlag[i]._MomentumAccFlag_momacc_file);
	    
	    if((fp2 = fopen(MAFile, "r")) == NULL)  
	      {
		fprintf(stdout, "%s: error while opening file.\n", MAFile);
		exit_(1);
	      }

	    bool found=false;
	    while(fgets(buffer2, 80, fp2) != NULL)
	      {
		if( (strstr(buffer2,"Negative")==0) && (!found)) {  buffer2[0]='\0'; continue;}
		if( (strstr(buffer2,"Negative")!=0) ) { found=true; buffer2[0]='\0'; continue;}
		
		fputs(buffer2, outf2);
		buffer2[0]='\0';
	      }
	    buffer2[0]='\0';
	    fclose(fp2);
	  }
	fclose(outf2);
	
      }
    MPI_Barrier(MPI_COMM_WORLD); 
    printf("process %d finished momentum acceptance.\n", myid);

#else
     MomentumAcceptance(UserCommandFlag[i]._MomentumAccFlag_momacc_file,
			UserCommandFlag[i]._MomentumAccFlag_istart, 
			UserCommandFlag[i]._MomentumAccFlag_istop,
			UserCommandFlag[i]._MomentumAccFlag_deltaminp, 
			UserCommandFlag[i]._MomentumAccFlag_deltamaxp,
			UserCommandFlag[i]._MomentumAccFlag_nstepp, 
			UserCommandFlag[i]._MomentumAccFlag_deltaminn,
			UserCommandFlag[i]._MomentumAccFlag_deltamaxn, 
			UserCommandFlag[i]._MomentumAccFlag_nstepn,
			UserCommandFlag[i]._MomentumAccFlag_nturn,
			UserCommandFlag[i]._MomentumAccFlag_zmax);
  #endif   

    /* restore the initial values*/
    globval.Cavity_on = cavityflag;
    globval.radiation = radiationflag;
  }

  // induced amplitude 
  else if(strcmp(CommandStr,"InducedAmplitudeFlag") == 0) {
   printf("\n Calculate induced amplitude: \n");
    InducedAmplitude(193L);
  }


  else if(strcmp(CommandStr,"EtaFlag") == 0) {
    // compute cod and Twiss parameters for different energy offsets
    for (int ii = 0; ii <= 40; ii++) {
      dP = -0.02 + 0.001 * ii;
      Ring_GetTwiss(chroma = false, dP); /* Compute and get Twiss parameters */
      printlatt("linlat_eta.out"); /* dump linear lattice functions into "linlat.dat" */
      getcod(dP, lastpos);
      //     printcod();
      prt_cod("cod.out", globval.bpm, true);
      //system("mv linlat.out linlat_ooo.out");
      sprintf(str1, "mv cod.out cod_%02d.out", ii);
      system(str1);
      sprintf(str1, "mv linlat.out linlat_%02d.out", ii);
      system(str1);
    }
  }

// track to get phase space
  else if(strcmp(CommandStr,"PhaseSpaceFlag") == 0) {
    bool cavityflag, radiationflag;
    /* record the initial values*/
    cavityflag = globval.Cavity_on;
    radiationflag = globval.radiation;

    /* set the dimension for the momentum tracking*/
    if (strcmp("6D", UserCommandFlag[i]._Phase_Dim) == 0) {
      globval.Cavity_on = true;
    } else if (strcmp("4D", UserCommandFlag[i]._Phase_Dim) == 0) {
      globval.Cavity_on = false;
    } else {
      printf("MomentumAccFlag: Error!!! _Phase_Dim must be '4D' or '6D'\n");
      exit_(1);
    };
    /* setting damping */
    if ( UserCommandFlag[i]._Phase_Damping == true) {
      globval.radiation = true;
    } else  {
      globval.radiation = false;
    }

    start = stampstart();
    Phase(UserCommandFlag[i]._Phase_phase_file,
          UserCommandFlag[i]._Phase_X, 
          UserCommandFlag[i]._Phase_Px, 
	  UserCommandFlag[i]._Phase_Y, 
	  UserCommandFlag[i]._Phase_Py, 
	  UserCommandFlag[i]._Phase_delta, 
	  UserCommandFlag[i]._Phase_ctau, 
	  UserCommandFlag[i]._Phase_nturn);
    printf("6D phase space tracking: \n the simulation time for phase space in tracy 3 is \n");
    stop = stampstop(start);

    /* restore the initial values*/
    globval.Cavity_on = cavityflag;
    globval.radiation = radiationflag;
  }
  
 
  

 else if (strcmp(CommandStr,"TouschekFlag") == 0 ||strcmp(CommandStr,"IBSFlag") == 0 ||
          strcmp(CommandStr,"TousTrackFlag") == 0 ){ 
  //  ????????????? NSRL version, Check with Soleil version "MomentumAcceptance"
  // IBS & TOUSCHEK
  int k, n_turns;
  double sigma_s, sigma_delta, tau, alpha_z, beta_z, gamma_z, eps[3];
  FILE *outf;
  const double Qb = 5e-9; // e charge in one bunch

  
//  else if(strcmp(CommandStr,"TouschekFlag") == 0) {
    double sum_delta[globval.Cell_nLoc + 1][2];
    double sum2_delta[globval.Cell_nLoc + 1][2];

    GetEmittance(globval.cav, true);

    // initialize momentum aperture arrays
    for (k = 0; k <= globval.Cell_nLoc; k++) {
      sum_delta[k][0] = 0.0;
      sum_delta[k][1] = 0.0;
      sum2_delta[k][0] = 0.0;
      sum2_delta[k][1] = 0.0;
    }

    globval.eps[Y_] = 0.008e-9;
    n_turns = 40;

    // get the twiss parameters
    alpha_z = -globval.Ascr[ct_][ct_] * globval.Ascr[delta_][ct_]
        - globval.Ascr[ct_][delta_] * globval.Ascr[delta_][delta_];
    beta_z = sqr(globval.Ascr[ct_][ct_]) + sqr(globval.Ascr[ct_][delta_]);
    gamma_z = (1 + sqr(alpha_z)) / beta_z;

    sigma_delta = sqrt(gamma_z * globval.eps[Z_]);
    sigma_s = sqrt(beta_z * globval.eps[Z_]);//50e-3
    beta_z = sqr(sigma_s) / globval.eps[Z_];
    alpha_z = sqrt(beta_z * gamma_z - 1);

    // INCLUDE LC (LC changes sigma_s and eps_z, but has no influence on sigma_delta)
    if (false) {
      double newLength, bunchLengthening;
      newLength = 50e-3;
      bunchLengthening = newLength / sigma_s;
      sigma_s = newLength;
      globval.eps[Z_] = globval.eps[Z_] * bunchLengthening;
      beta_z = beta_z * bunchLengthening;
      gamma_z = gamma_z / bunchLengthening;
      alpha_z = sqrt(beta_z * gamma_z - 1); // this doesn't change
    }

    Touschek(Qb, globval.delta_RF, globval.eps[X_], globval.eps[Y_],
        sigma_delta, sigma_s);

 // Intra Beam Scattering(IBS)
    if (strcmp(CommandStr,"IBSFlag") == 0) {
      // initialize eps_IBS with eps_SR
      for (k = 0; k < 3; k++)
        eps[k] = globval.eps[k];
      for (k = 0; k < 20; k++) //prototype (looping because IBS routine doesn't check convergence)
        IBS(Qb, globval.eps, eps, alpha_z, beta_z);
    }

    // TOUSCHEK TRACKING
    // Calculate Touschek lifetime
    // with the momentum acceptance which is determined by 
    // the RF acceptance delta_RF and the momentum aperture 
    // at each element location which is tracked over n turns, 
    //the vacuum chamber is read from the file "chamber_file"  
    // finally, write the momentum acceptance to file "mom_aper.out".
    if (strcmp(CommandStr,"TousTrackFlag") == 0) {
        // globval.Aperture_on = true;
        // ReadCh(UserCommandFlag[i].chamber_file);
      //      LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1);
      tau = Touschek(Qb, globval.delta_RF, false, globval.eps[X_],
          globval.eps[Y_], sigma_delta, sigma_s, n_turns, true, sum_delta,
          sum2_delta); //the TRUE flag requires apertures loaded

      printf("Touschek lifetime = %10.3e hrs\n", tau / 3600.0);

      outf = file_write("mom_aper.out");
      for (k = 0; k <= globval.Cell_nLoc; k++)
        fprintf(outf, "%4d %7.2f %5.3f %6.3f\n", k, Cell[k].S, 1e2
            * sum_delta[k][0], 1e2 * sum_delta[k][1]);
      fclose(outf);
    }

  }
  
  /*
To do ID correction.
Based on parts of functions get_param( ) & ID_corr(), etc in nsls-ii_lib.cc 
*/
  else if(strcmp(CommandStr,"IDCorrFlag") == 0) {
   
    int k = 0;
  
    cout << endl;
    cout << "************  Begin ID matching:   **********" <<endl;
 
  // read the family index of quadrupoles used for ID correction
    if (N_calls > 0) {
      if (N_Fam > N_Fam_max) {
        printf("get_param: N_Fam > N_Fam_max: %d (%d)\n",
	       N_Fam, N_Fam_max);
        exit_(0);
      }

    for (k = 0; k < N_Fam; k++)
      Q_Fam[k] = ElemIndex(IDCq_name[k]);
    }

    //For debug
    if(trace){
      cout << "scl_dbetax = " << scl_dbetax << "    scl_dbetay = " << scl_dbetay <<endl;
      cout << "scl_dnux = " << scl_dnux     << "    scl_dnuy = " << scl_dnuy << endl;
      cout << "scl_nux = " << scl_nux       << "    scl_nuy = " << scl_nuy << endl;
      cout << "ID_step = " << ID_step <<endl;
      cout << "N_calls = " << N_calls << "    N_steps = " << N_steps <<endl;
      cout << "Number of quadrupole families used for ID correction:   " << N_Fam << endl; 
      cout << "Quadrupoles used for ID correction: " << endl;
      for (int k = 0; k < N_Fam; k++)
        printf("%d\n",Q_Fam[k]);
    } 
    
    //initialization
    ini_ID_corr();
    printlatt("linlat00.out");

    ID_corr0();
    printlatt("linlat01.out");

  // exit(0);

    ini_ID_corr();
    printlatt("linlat0.out");

    ID_corr(N_calls,N_steps);
    printlatt("linlat1.out");

  }
  else{
    printf("Command string %s is invalid!!!!!\n ",CommandStr);
    exit_(1);
   }//give an error message
 }//end of looking for user defined flag

 
 #if MPI_EXEC
  MPI_Finalize();   //finish parallel computation
 #endif
  
  return 0;
}//end of main()