soltracy.cc

/*
 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[]) {
  
  /*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 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;

 
  long i=0L; //initilize the for loop to read command string
  char CommandStr[max_str];
  double nux = 0, nuy = 0, ksix = 0, ksiy = 0;
  bool chroma;
  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. 
  UserCommand UserCommandFlag[500];
 

  // for time handling
  uint32_t start, stop;

 
  /************************************************************************
   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);
  }

  /************************************************************************
    print files, very important file for debug
   *************************************************************************/
  //print flat file with all the design values of the lattice,
  prtmfile("flat_file.dat");
  // print close orbit 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); 
    
    //print the twiss paramters in a file defined by the name
    if(strcmp(CommandStr,"PrintTwissFlag") == 0) {
      cout << "\n";
      cout << "print the twiss parameters to file: "<< twiss_file << "\n";
      printlatt(twiss_file);  
   }
   //print the close orbit
   else if(strcmp(CommandStr,"PrintCODFlag") == 0) {
      cout << "\n";
      cout << "print the close orbit to file: "<< cod_file << "\n";
      getcod(dP, lastpos);
      prt_cod(cod_file, globval.bpm, true);
   }
    
    
    //turn on flag for quadrupole fringe field
   else if(strcmp(CommandStr,"QuadFringeOnFlag") == 0) {
      globval.quad_fringe = true;
      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(ElemIndex("cav"))
        / 1e6);
    set_RFVoltage(ElemIndex("cav"), UserCommandFlag[i].RFvolt);
    printf("    New RF voltage is: %lf [MV]\n", get_RFVoltage(ElemIndex("cav"))
        / 1e6);
  }

  // Chamber factory
  else if(strcmp(CommandStr,"ReadChamberFlag") == 0) {
    ReadCh(chamber_file); /* read vacuum chamber from a file "Apertures.dat" , soleil version*/
  PrintCh(); // print chamber into chamber.out
  }

 
  // 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 */
  }

  //set multipoles in all magnets
  // read multipole error from a file
  else if(strcmp(CommandStr,"ReadMultipoleFlag") == 0) {
    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();
  }

 

  /******************************************************************************************/
  // COMPUTATION PART after setting the model
  /******************************************************************************************/
  // computes TuneShift with amplitudes
  else if(strcmp(CommandStr,"AmplitudeTuneShiftFlag") == 0) {
      TunesShiftWithAmplitude(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_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");
    fmap(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);
  }

  // Compute FMA dp
  else if(strcmp(CommandStr,"FmapdpFlag") == 0) {
    printf("\n begin Fmap calculation for off momentum particles: \n");
    fmapdp(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);
  }

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

    /* set the dimension for the momentum tracking*/
    if (strncmp("6D", UserCommandFlag[i].TrackDim, 2) == 0) {
      globval.Cavity_on = true;
      globval.radiation = true;
    } else if (strncmp("4D", UserCommandFlag[i].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*/
    printf("\n Calculate momentum acceptance: \n");
    MomentumAcceptance(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);

    /* 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 (strncmp("6D", UserCommandFlag[i]._Phase_Dim, 2) == 0) {
      globval.Cavity_on = true;
    } else if (strncmp("4D", UserCommandFlag[i]._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 ( UserCommandFlag[i]._Phase_Damping == true) {
      globval.radiation = true;
    } else  {
      globval.radiation = false;
    }

    start = stampstart();
    Phase(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("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(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 (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(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);
    }

  }
  else
    printf("Wrong!!!!!");
 }//end of looking for user defined flag

  
  return 0;
}//end of main()