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

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

extern bool freq_map;

#include "tracy_lib.h"

//***************************************************************************************
//
//  MAIN CODE
//
//****************************************************************************************
int main(int argc, char *argv[]) {
  const long seed = 1121;
  iniranf(seed);
  setrancut(2.0);

  // turn on globval.Cavity_on and globval.radiation to get proper synchr radiation damping
  // IDs accounted too if: wiggler model and symplectic integrator (method = 1)
  globval.H_exact = false;
  /*The following values are set in the Read_lattice( ) in soleilcommon.cc*/
  //  globval.quad_fringe = false;              // quadrupole fringe field
  // globval.radiation   = false;             // synchrotron radiation
  // globval.emittance   = false;             // emittance
  // globval.pathlength  = false;
  //  globval.bpm         = 0;

  //  const double  x_max_FMA = 10e-3,  delta_FMA = 10e-2;
  //  const int     n_x = 801, n_dp = 80, n_tr = 2048;

  double nux = 0, nuy = 0, ksix = 0, ksiy = 0;

  bool chroma;
  double dP = 0.0;
  long lastpos = -1L;
  char str1[S_SIZE];

  // for time handling
  uint32_t start, stop;

  /************************************************************************
   read in files and flags
   *************************************************************************/
  if (argc > 1) {
    read_script(argv[1], true);
  } else {
    fprintf(stdout, "Not enough parameters\nSyntax is program parameter file\n");
    exit_(1);
  }

  /************************************************************************
   writes flat file with all the design values of the lattice, very important file for debug
   *************************************************************************/
  prtmfile("flat_file.dat");

  // print cod
  getcod(dP, lastpos);
  prt_cod("cod.out", globval.bpm, true);

  // Flag factory

  //set RF voltage  
  if (RFvoltageFlag == true) {
    printf("\nSetting RF voltage:\n");
    printf("    Old RF voltage is: %lf [MV]\n", get_RFVoltage(ElemIndex("cav"))
        / 1e6);
    set_RFVoltage(ElemIndex("cav"), RFvolt);
    printf("    New RF voltage is: %lf [MV]\n", get_RFVoltage(ElemIndex("cav"))
        / 1e6);
  }

  // Chamber factory
 // if (ChamberFlag == false)
  //  ChamberOff(); // turn off vacuum chamber setting, use the default +-1 meters
  if (ReadChamberFlag == true)
    ReadCh(chamber_file); /* read vacuum chamber from a file "Apertures.dat" , soleil version*/
  PrintCh(); // print chamber into chamber.out

  //get_matching_params_scl(); // get tunes and beta functions at entrance
  get_alphac2(); // compute up to 3rd order mcf
  //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


  // compute tunes by tracking (should be the same as by tps)
  if (TuneTracFlag == true) {
    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)
  if (ChromTracFlag == true) {
    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) {
    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(qf), 1)].Elem.M->PB[HOMmax + 2];
    qd_bn = Cell[Elem_GetPos(ElemIndex(qd), 1)].Elem.M->PB[HOMmax + 2];

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

    /* integrated field strength after fitting*/
    qf_bn_fit = Cell[Elem_GetPos(ElemIndex(qf), 1)].Elem.M->PB[HOMmax + 2];
    qd_bn_fit = Cell[Elem_GetPos(ElemIndex(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", qf, qf_bn, qd, qd_bn);
    printf("After the fitting, the quadrupole field strengths are: \n");
    printf(" %s = %f,    %s = %f\n", qf, qf_bn_fit, 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) {
    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(qf1), 1)].Elem.M->PB[HOMmax + 2];
    qf2_bn = Cell[Elem_GetPos(ElemIndex(qf2), 1)].Elem.M->PB[HOMmax + 2];
    qd1_bn = Cell[Elem_GetPos(ElemIndex(qd1), 1)].Elem.M->PB[HOMmax + 2];
    qd2_bn = Cell[Elem_GetPos(ElemIndex(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",
        qf1, qf2, qd1, qd2, targetnux, targetnuz);
    FitTune4(ElemIndex(qf1), ElemIndex(qf2), ElemIndex(qd1), ElemIndex(qd2),
        targetnux, targetnuz);

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

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

  /* fit the chromaticities*/
  if (FitChromFlag == true) {

    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(sxm1), 1)].Elem.M->PB[HOMmax + 3];
    sxm2_bn = Cell[Elem_GetPos(ElemIndex(sxm2), 1)].Elem.M->PB[HOMmax + 3];

    /* fit the chromaticities*/
    fprintf(
        stdout,
        "\n Fitting chromaticities: %s %s, targetksix = %f,  targetksiz = %f\n",
        sxm1, sxm2, targetksix, targetksiz);
    FitChrom(ElemIndex(sxm1), ElemIndex(sxm2), targetksix, targetksiz);

    sxm1_bn_fit = Cell[Elem_GetPos(ElemIndex(sxm1), 1)].Elem.M->PB[HOMmax + 3];
    sxm2_bn_fit = Cell[Elem_GetPos(ElemIndex(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", sxm1, sxm1_bn, sxm2, sxm2_bn);
    printf("After the fitting, the sextupole field strengths are: \n");
    printf("    %s = %f,    %s = %f\n", sxm1, sxm1_bn_fit, sxm2, sxm2_bn_fit);

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

  // coupling calculation
  if (CouplingFlag == true) {
    Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
    printlatt(); /* 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) {
    prtmfile("flat_file.dat"); //print the elements without rotation errors
    SetErr(err_seed, 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(); /* 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
  if (ErrorCoupling2Flag == true) {
    prtmfile("flat_file.dat"); //print the elements without rotation errors
    SetErr2(err_seed, 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(); /* 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 */
  }

  // WARNING Fit tunes and chromaticities before applying errors !!!!
  //set multipoles in all magnets
  // read multipole error from a file
  if (ReadMultipoleFlag == true) {
    fprintf(stdout,
        "\n Read Multipoles file for lattice with thick sextupoles \n");
    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();
  }

  //Old version to multipole errors in soleil lattice, is obsoleted
  if (MultipoleFlag == true) {
    if (ThinsextFlag == true) {
      fprintf(stdout,
          "\n Setting Multipoles for lattice with thin sextupoles \n");
      Multipole_thinsext(fic_hcorr, fic_vcorr, fic_skew); /* for thin sextupoles */
      prtmfile("flat_file_errmultipole.dat"); // writes flat file   /* very important file for debug*/
      Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
      printglob();
    } else {
      fprintf(stdout,
          "\n Setting Multipoles for lattice with thick sextupoles \n");
      Multipole_thicksext(fic_hcorr, fic_vcorr, fic_skew); /* for thick sextupoles */
      prtmfile("flat_file_errmultipole.dat"); // writes flat file   /* very important file for debug*/
      Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
      printglob();
    }
  }

  /******************************************************************************************/
  // COMPUTATION PART after setting the model
  /******************************************************************************************/

  // computes TuneShift with amplitudes
  if (AmplitudeTuneShiftFlag == true) {
   // if (ChamberFlag == true) {
      TunesShiftWithAmplitude(_AmplitudeTuneShift_nxpoint, _AmplitudeTuneShift_nypoint,
          _AmplitudeTuneShift_nturn, _AmplitudeTuneShift_xmax,
          _AmplitudeTuneShift_ymax, _AmplitudeTuneShift_delta);
      //NuDx(31L,21L,516L,0.025,0.005,dP);
    } //else { // Utility ?
     // TunesShiftWithAmplitude(50L, 30L, 516L, 0.035, 0.02, dP);
   // }

  //}
  if (EnergyTuneShiftFlag == true) {
   // if (ChamberFlag == true) {
      TunesShiftWithEnergy(_EnergyTuneShift_npoint, _EnergyTuneShift_nturn,
          _EnergyTuneShift_deltamax);
      //NuDp(31L,1026L,0.06);
    }// else { // utility ?
     // TunesShiftWithEnergy(31L, 1026L, 0.06);
    //}

  //}

  // Computes FMA
  if (FmapFlag == true) {
    fmap(_FmapFlag_nxpoint, _FmapFlag_nypoint, _FmapFlag_nturn, _FmapFlag_xmax,
        _FmapFlag_ymax, _FmapFlag_delta, _FmapFlag_diffusion);
  }

  // Compute FMA dp
  if (FmapdpFlag == true) {
    fmapdp(_FmapdpFlag_nxpoint, _FmapdpFlag_nepoint, _FmapdpFlag_nturn,
        _FmapdpFlag_xmax, _FmapdpFlag_emax, _FmapdpFlag_z,
        _FmapdpFlag_diffusion);
  }

  if (CodeComparaisonFlag) {
    fmap(200, 100, 1026, -32e-3, 7e-3, 0.0, true);
  }

  // MOMENTUM ACCEPTANCE
  if (MomentumAccFlag == true) {

    bool cavityflag, radiationflag;
    /* record the initial values*/
    cavityflag = globval.Cavity_on;
    radiationflag = globval.radiation;

    /* set the dimension for the momentum tracking*/
    if (strncmp("6D", TrackDim, 2) == 0) {
      globval.Cavity_on = true;
      globval.radiation = true;
    } else if (strncmp("4D", TrackDim, 2) == 0) {
      globval.Cavity_on = false;
      globval.radiation = false;
    } else {
      printf("MomentumAccFlag: Error!!! DimTrack must be '4D' or '6D'\n");
      exit_(1);
    };

    /* calculate momentum acceptance*/
    MomentumAcceptance(_MomentumAccFlag_istart, _MomentumAccFlag_istop,
        _MomentumAccFlag_deltaminp, _MomentumAccFlag_deltamaxp,
        _MomentumAccFlag_nstepp, _MomentumAccFlag_deltaminn,
        _MomentumAccFlag_deltamaxn, _MomentumAccFlag_nstepn);

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

  // induced amplitude 
  if (InducedAmplitudeFlag == true) {
    InducedAmplitude(193L);
  }

  if (EtaFlag == true) {
    // 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(); /* 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);
    }
  }

  if (PhaseSpaceFlag == true) {
    bool cavityflag, radiationflag;
    /* record the initial values*/
    cavityflag = globval.Cavity_on;
    radiationflag = globval.radiation;

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

    start = stampstart();
    Phase(_Phase_X, _Phase_Px, _Phase_Y, _Phase_Py, _Phase_delta, _Phase_ctau, _Phase_nturn);
    printf("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;
  }
  
 

  //  ????????????? 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

  if (TouschekFlag == true) {
    double sum_delta[globval.Cell_nLoc + 1][2];
    double sum2_delta[globval.Cell_nLoc + 1][2];

    GetEmittance(ElemIndex("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 (IBSFlag == true) {
      // 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 (TousTrackFlag == true) {
      globval.Aperture_on = true;
      ReadCh(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);
    }

  }

  return 0;
}