Skip to content
Snippets Groups Projects
  • nadolski's avatar
    26f4b74e
    REname function · 26f4b74e
    nadolski authored
    Phase space introduced with tracking arounf 6D COD
    Phase space available in user interface
    26f4b74e
    History
    REname function
    nadolski authored
    Phase space introduced with tracking arounf 6D COD
    Phase space available in user interface
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
max4.cc 15.06 KiB
#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;
  const bool All = TRUE;
  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; 

  
 
  // overview, on energy: 25-15
  //const double  x_max_FMA = 20e-3, y_max_FMA = 10e-3; //const x_max_FMA = 25e-3, y_max_FMA = 15e-3;
  //const int     n_x = 80, n_y = 80, n_tr = 2048;
  // overview, off energy: 10-10
  const double  x_max_FMA = 10e-3, delta_FMA = 10e-2;
  const int     n_x = 80, n_dp = 80, n_tr = 2048;
  //
  // zoom, on energy: 8-2.5
  //const double  x_max_FMA = 8e-3, y_max_FMA = 2.5e-3;
  //const int     n_x = 64, n_y = 15, n_tr = 2048;
  // zoom, off energy: 7-3
  //const double  x_max_FMA = 3e-3, delta_FMA = 7e-2;
  //const int     n_x = 28, n_dp = 56, n_tr = 2048;
  
  double nux=0.0, nuz=0.0, ksix=0.0, ksiz=0.0;
  
  bool chroma;
  double dP = 0.0;
  long lastpos = -1L;
  char str1[S_SIZE];
  
  /************************************************************************
      start read in files and flags
  *************************************************************************/
  read_script(argv[1], true);


 /************************************************************************
    end  read in files and flags
  *************************************************************************/
    
  
  
  
  
  
//  if (true)
//  //  Read_Lattice("/home/nadolski/codes/tracy/maille/soleil/solamor2_tracy3"); 
//    Read_Lattice(argv[1]); //sets some globval params
//  else
//    rdmfile("flat_file.dat"); //instead of reading lattice file, get data from flat file 

  //no_sxt(); //turns off sextupoles
 // Ring_GetTwiss(true, 0e-2); //gettwiss computes one-turn matrix arg=(w or w/o chromat, dp/p) 
  //get_matching_params_scl();
  //get_alphac2();
  //GetEmittance(ElemIndex("cav"), true);
  
//prt_lat("linlat.out", globval.bpm, true);  /* print lattice file for nsrl-ii*/
prtmfile("flat_file.dat"); // writes flat file   /* very important file for debug*/
//prt_chrom_lat(); //writes chromatic functions into chromlat.out
//  printlatt();  /* print out lattice functions */
 

/* print lattice file */
//  prt_lat("linlatBNL.out", globval.bpm, All); // BNL print for all elements
  printlatt();  /* SOLEIL print out lattice functions */
  printglob();
  
  
  
    // Flag factory
//  bool TuneTracFlag = true; 
//  bool ChromTracFlag = true;
  
  
  //*************************************************************
  //=============================================================

    
    // 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("Apertures.dat"); /* read vacuum chamber from a file "Apertures.dat" , soleil version*/
//LoadApers("Apertures.dat", 1.0, 1.0);  /* read vacuum chamber definition for bnl */
  PrintCh();
 
 
  // 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){
    GetChromTrac(2L, 1026L, 1e-5, &ksix, &ksiz);
    fprintf(stdout,"From tracking: ksix= % f ksiz= % f \n",ksix,ksiz);
  }


  if (FitTuneFlag == true){
    fprintf(stdout, "\n Fitting tunes\n");
    FitTune(ElemIndex("qp7"),ElemIndex("qp9"), 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\n");
    FitChrom(ElemIndex("sx9"),ElemIndex("sx10"), targetksix, targetksiz);
    Ring_GetTwiss(chroma=true, 0.0);  /* Compute and get Twiss parameters */
    printglob();                      /* print parameter list */
  }

    //SetKLpar(ElemIndex("QT"), 1, 2L,   0.001026770838382);
  
  // coupling calculation
   if (CouplingFlag == true){
     Ring_GetTwiss(chroma=true, 0.0);  /* Compute and get Twiss parameters */
     printlatt();                      /* dump linear lattice functions into "linlat.dat" */
  //   Coupling_Edwards_Teng();
     printglob();   /* print parameter list */
   }

  // add coupling by random rotating of the quadrupoles
  if (ErrorCouplingFlag == true){
    SetErr();
    Ring_GetTwiss(chroma=true, 0.0);  /* Compute and get Twiss parameters */
    printlatt();                      /* dump linear lattice functions into "linlat.dat" */
//    Coupling_Edwards_Teng();
    printglob();   /* print parameter list */
  }

  // WARNING Fit tunes and chromaticities before applying errors !!!!
  //set multipoles in all magnets
  if (MultipoleFlag == true ){
    if (ThinsextFlag ==true){
      fprintf(stdout, "\n Setting Multipoles for lattice with thin sextupoles \n");
      Multipole_thinsext();  /* 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();  /* for thick sextupoles */
      
      Ring_GetTwiss(chroma=true, 0.0);  /* Compute and get Twiss parameters */
      printglob(); 
    }
  }
                       /* print parameter list */
  

  // PX2 chicane
//  if (PX2Flag ==true){
//  setPX2chicane(); 
//  //get closed orbit    
//  getcod (0.0, &lastpos);
//  printcod();
//  Ring_GetTwiss(chroma=true, 0.0);  /* Compute and get Twiss parameters */
//  printglob();                      /* print parameter list */
//  }
  
 // Computes FMA
  if (FmapFlag == true){
    if (ChamberFlag == true ){
      if (ExperimentFMAFlag == true)
         fmap(40,12,258,-20e-3,5e-3,0.0,true); // for experimental
      if (DetailedFMAFlag == true)
        fmap(100,50,1026,20e-3,5e-3,0.0,true);
      }
      else{
        if (ExperimentFMAFlag == true)
          fmap(40,12,258,-32e-3,5e-3,0.0,true);
        if (DetailedFMAFlag == true)
          fmap(200,100,1026,32e-3,7e-3,0.0,true);
      }
  }
  
  if (CodeComparaisonFlag){
          // SOLEIL
          fmap(100,50,1026,32e-3,7e-3,0.0,true);
          //fmap(200,100,1026,-32e-3,7e-3,0.0,true);
  }
  if (MomentumAccFlag == true){
    //MomentumAcceptance(10L, 28L, 0.01, 0.05, 4L, -0.01, -0.05, 4L);
     MomentumAcceptance(1L, 28L, 0.01, 0.05, 40L, -0.01, -0.05, 40L);
  //  MomentumAcceptance(1L, 108L, 0.01, 0.05, 100L, -0.01, -0.05, 100L);
  }

  // computes Tuneshift with amplitudes
  if (TuneShiftFlag == true){
    if (ChamberFlag == true ){
      TunesShiftWithAmplitude(31L,21L,516L,0.025,0.005,dP);
      //NuDp(31L,516L,0.06);
      //NuDp(31L,1026L,0.06);
      }
      else{
        TunesShiftWithAmplitude(50L,30L,516L,0.035,0.02,dP);
        TunesShiftWithEnergy(31L,1026L,0.06);
      }

  }
 
//  if (SigmaFlag == true){printsigma();
//  }
//  

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

  
  //*********************************************************************************
  //---------------------------------------------------------------------------------------------------------------------------------
  
  // Delicated for max4 lattice. To load alignment error files and do correction
  if (false) {
    //prt_cod("cod.out", globval.bpm, true); //prints a specific closed orbit with corrector strengths
    
    
    globval.bpm = ElemIndex("bpm_m");                   //definition for max4 lattice, bpm
  //  globval.bpm = ElemIndex("bpm");                         
    globval.hcorr = ElemIndex("corr_h"); globval.vcorr = ElemIndex("corr_v");    //definition for max4 lattice, correctors
   // globval.hcorr = ElemIndex("ch"); globval.vcorr = ElemIndex("cv");
    globval.gs = ElemIndex("GS"); globval.ge = ElemIndex("GE");   //definition for max4 lattice, girder maker
    
    
    //compute response matrix (needed for OCO)
    gcmat(globval.bpm, globval.hcorr, 1);  gcmat(globval.bpm, globval.vcorr, 2);
     

    //print response matrix (routine in lsoc.cc)
    //prt_gcmat(globval.bpm, globval.hcorr, 1);  prt_gcmat(globval.bpm, globval.vcorr, 2);
       
    //gets response matrix, does svd, evaluates correction for N seed orbits
    //get_cod_rms_scl(dx, dy, dr, n_seed)
    //get_cod_rms_scl(100e-6, 100e-6, 1.0e-3, 100); //trim coils aren't reset when finished
       

    //for alignments errors check LoadAlignTol (in nsls_ii_lib.cc) and AlignErr.dat
    //CorrectCOD_N("/home/simon/projects/src/lattice/AlignErr.dat", 3, 1, 1);
   
    //for field errors check LoadFieldErr(in nsls_ii_lib.cc) and FieldErr.dat
    //LoadFieldErr("/home/simon/projects/src/lattice/FieldErr.dat", true, 1.0, true);
   
    //for alignments errors check LoadAlignTol (in nsls_ii_lib.cc) and AlignErr.dat
    //LoadAlignTol("/home/simon/projects/src/lattice/AlignErr.dat", true, 1.0, true, 1);
    //LoadAlignTol("/home/simon/projects/out/20091126/AlignErr.dat", true, 1.0, true, 1);
    //prt_cod("cod_err.out", globval.bpm, true); //prints a specific closed orbit with corrector strengths
   
      
    // delicated for max4 lattice
    //load alignment errors and field errors, correct orbit, repeat N times, and get statistics
    get_cod_rms_scl_new(1); //trim coils aren't reset when finished
  
    
    //for aperture limitations use LoadApers (in nsls_ii_lib.cc) and Apertures.dat
    //globval.Aperture_on = true;
    //LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1);
    
  }



//*******************************************************************************
//-------------------------------------------------------------------------------------------------------------------------------------  

  // Call nsls-ii_lib.cc
  // tune shift with amplitude 
  double delta = 0;
  if (false) {
    cout << endl;
    cout << "computing tune shifts" << endl;
    dnu_dA(20e-3, 10e-3, 0.0);
    get_ksi2(delta); // this gets the chromas and writes them into chrom2.out
 // get_ksi2(5.0e-2); // this gets the chromas and writes them into chrom2.out
  }
  
  if (false) {
    //fmap(n_x, n_y, n_tr, x_max_FMA, y_max_FMA, 0.0, true, false);
//    fmapdp(n_x, n_dp, n_tr, x_max_FMA, -delta_FMA, 1e-3, true, false); // always use -delta_FMA (+delta_FMA appears broken)
    fmapdp(n_x, n_dp, n_tr, x_max_FMA, -delta_FMA, 1e-3, true); // always use -delta_FMA (+delta_FMA appears broken)
  } else {
    globval.Cavity_on = true; // this gives longitudinal motion
    globval.radiation = false; // this adds ripple around long. ellipse (needs many turns to resolve damp.)
    //globval.Aperture_on = true;
    //LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1);
    //get_dynap_scl(delta, 512);
  }
  
  

  
  
  
  
  //
  // 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;
  
  if (false) {
    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;
    }
    
    //sigma_delta = 7.70e-04;      //410:7.70e-4,  411:9.57e-4,  412:9.12e-4
    //globval.eps[X_] = 0.326e-9;  //410:3.26e-10, 411:2.63e-10, 412:2.01e-10
    globval.eps[Y_] = 0.008e-9;
    //sigma_s = 9.73e-3;           //410:9.73e-3,  411:12.39e-3, 412:12.50e-3/10.33e-3
    //globval.eps[Z_] = sigma_delta*sigma_s;
    //globval.delta_RF given by cav voltage in lattice file
     //globval.delta_RF = 6.20e-2; //410:6.196e-2, 411:5.285e-2, 412:4.046e-2/5.786e-2
    n_turns = 490;               //410:490(735), 411:503(755), 412:439(659)/529(794)
    

    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
    }
    
    //globval.eps[X_] = 0.362e-9;
    //sigma_delta = 1.04e-3;
    //sigma_s = 14.8e-3;
    Touschek(Qb, globval.delta_RF, globval.eps[X_], globval.eps[Y_],
	     sigma_delta, sigma_s);
    
    
    // IBS
    if (false) {       
      // 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
    if (false) {       
      globval.Aperture_on = true;
      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);
    }
  
  }
}