soltracy.cc

#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; 
  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, nuz=0, ksix=0, ksiz=0;

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

  // for time handling
  uint32_t         start, stop;

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



 /************************************************************************
    end  read in files and flags
  *************************************************************************/
  prtmfile("flat_file.dat"); // writes flat file   /* very important file for debug*/
  printlatt();  /* SOLEIL print out lattice functions */
  printglob();
  
  
  
  // 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 one
  else if (ChamberNoU20Flag == true)
     DefineChNoU20();  // using vacuum chamber setting but without undulator U20
  else 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 DA)
  if (TuneTracFlag == true) {
    GetTuneTrac(1026L, 0.0, &nux, &nuz);
    fprintf(stdout,"From tracking: nux = % f nuz = % f \n",nux,nuz);
  }

  // compute chromaticities by tracking (should be the same as by DA)
  if (ChromTracFlag == true){
	start = stampstart();
	GetChromTrac(2L, 1026L, 1e-5, &ksix, &ksiz);
	stop = stampstop(start);

    fprintf(stdout,"From tracking: ksix= % f ksiz= % f \n",ksix,ksiz);
  }


  if (FitTuneFlag == true){
    fprintf(stdout, "\n Fitting tunes: %s %s, targetnux = %f, targetnuz = %f \n",qm1,qm2,targetnux,targetnuz);
    FitTune(ElemIndex(qm1),ElemIndex(qm2), targetnux, targetnuz);
    Ring_GetTwiss(chroma=true, 0.0);  /* Compute and get Twiss parameters */
    printglob();                      /* print parameter list */
  }

  if (FitChromFlag == true){
    fprintf(stdout, "\n Fitting chromaticities: %s %s, targetksix = %f,  targetksiz = %f\n",sxm1,sxm2,targetksix,targetksiz);
    FitChrom(ElemIndex(sxm1),ElemIndex(sxm2), targetksix, targetksiz);
    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);
     Coupling_Edwards_Teng();
     printglob();   /* print parameter list */
   }

  // add coupling by random rotating of the quadrupoles
  if (ErrorCouplingFlag == true){
    SetErr(err_seed,err_rms);
    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);
    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);
  }
  
  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 */
      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 */
      Ring_GetTwiss(chroma=true, 0.0);  /* Compute and get Twiss parameters */
      printglob(); 
    }
  }

 //first print the full lattice with error as a flat file
  prtmfile("flat_file_error.dat"); // writes flat file   /* very important file for debug*/
  printlatt();  /* SOLEIL print out lattice functions */
  printglob();
  
 /******************************************************************************************/
 // COMPUTATION PART after setting the model
 /******************************************************************************************/

  

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

  }
  if (EnergyTuneShiftFlag == true){
    if (ChamberFlag == true ){
    	NuDp(_EnergyTuneShift_npoint,
    	     _EnergyTuneShift_nturn, _EnergyTuneShift_deltamax);
      //NuDp(31L,1026L,0.06);
      }
      else{ // utility ?
        NuDp(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){
	  MomentumAcceptance(
	  _MomentumAccFlag_istart, _MomentumAccFlag_istop,
	  _MomentumAccFlag_deltaminp, _MomentumAccFlag_deltamaxp, _MomentumAccFlag_nstepp,
	  _MomentumAccFlag_deltaminn, _MomentumAccFlag_deltamaxn, _MomentumAccFlag_nstepn);
  }


  // 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){
		start = stampstart();
		Phase(0.001,0,0.001,0,0.001,0, 1);
	    printf("the simulation time for phase space in tracy 3 is \n");
		stop = stampstop(start);
  }
  return 0;
}