Skip to content
Snippets Groups Projects
max4.cc 15.3 KiB
Newer Older
zhang's avatar
zhang committed
#define ORDER 1          

#include "tracy_lib.h"

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

zhang's avatar
ZJ  
zhang committed
extern bool  freq_map;



//***************************************************************************************
//
//  MAIN CODE
//
//****************************************************************************************
zhang's avatar
zhang committed
 int main(int argc, char *argv[])
{

  const long  seed = 1121;
nadolski's avatar
nadolski committed
  const bool All = TRUE;
zhang's avatar
zhang committed
  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.Cavity_on   = false;                 // include RF cavity 
  globval.radiation   = false;                // synchrotron radiation
  globval.emittance   = false;                 // emittance
  globval.pathlength  = false; 
  globval.bpm         = 0;
  
  // 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;
  
zhang's avatar
ZJ  
zhang committed
  double nux = 0.0 , nuz = 0.0, ksix = 0.0, ksiz =  0.0;
nadolski's avatar
nadolski committed
  bool chroma;
  double dP = 0.0;
  long lastpos = -1L;
  char str1[S_SIZE];
zhang's avatar
zhang committed
  if (true)
nadolski's avatar
nadolski committed
    Read_Lattice("/home/nadolski/codes/tracy/maille/soleil/solamor2_tracy3"); 
  //  Read_Lattice(argv[1]); //sets some globval params
zhang's avatar
zhang committed
  else
    rdmfile("flat_file.dat"); //instead of reading lattice file, get data from flat file
zhang's avatar
ZJ  
zhang committed
    
    
  globval.quad_fringe = true;                // include quadrupole fringe field
nadolski's avatar
nadolski committed
  /*globval.bpm = ElemIndex("bpm");  //definition for max4 lattice, bpm                  
  globval.hcorr = ElemIndex("ch"); 
  globval.vcorr = ElemIndex("cv");*/   
zhang's avatar
zhang committed

zhang's avatar
ZJ  
zhang committed
  
nadolski's avatar
nadolski committed
  //no_sxt(); //turns off sextupoles
zhang's avatar
ZJ  
zhang committed
  
nadolski's avatar
nadolski committed
  Ring_GetTwiss(true, 0.0); 
  printglob(); //gettwiss computes one-turn matrix arg=(w or w/o chromat, dp/p)
  get_matching_params_scl();
  get_alphac2(); // compute up to 3rd order mcf
  GetEmittance(ElemIndex("cav"), true);
zhang's avatar
ZJ  
zhang committed

nadolski's avatar
nadolski committed
/* print lattice file */
  prt_lat("linlatBNL.out", globval.bpm, All); // BNL print for all elements
  printlatt();  /* SOLEIL print out lattice functions */
zhang's avatar
ZJ  
zhang committed
  
nadolski's avatar
nadolski committed
  prtmfile("flat_file.dat"); // writes flat file
  prt_chrom_lat(); //writes chromatic functions into chromlat.out
zhang's avatar
zhang committed

zhang's avatar
ZJ  
zhang committed
  
  // Using tracking to get tunes and chromaticities
nadolski's avatar
nadolski committed
  if (false) {
    GetTuneTrac(1026L, 0.0, &nux, &nuz);
    fprintf(stdout,"From tracking: nux = % f nuz = % f \n",nux,nuz);
    GetChromTrac(2L, 1026L, 1e-5, &ksix, &ksiz);
    fprintf(stdout,"From tracking: ksix= % f ksiz= % f \n",ksix,ksiz);
zhang's avatar
ZJ  
zhang committed
  }
  

nadolski's avatar
nadolski committed
    // Flag factory
  bool TuneTracFlag = true, ChromTracFlag = true;
  bool FmapFlag = false, ExperimentFMAFlag = false, DetailedFMAFlag = false;
  bool TuneShiftFlag = false;
  bool ErrorCouplingFlag = false;  bool CouplingFlag = false;
  bool MomentumAccFlag = false;
  bool MultipoleFlag = false;
  bool FitTuneFlag  = false; double targetnux = 18.202, targetnuz = 10.317;
  bool FitChromFlag = false; double targetksix = 2.0, targetksiz = 2.6;
  bool ChamberFlag = false, U20ChamberFlag = false;
  bool GirderErrorFlag = false;
  bool SigmaFlag = false;
  bool PX2Flag = false;
  bool InducedAmplitudeFlag = false;
  bool CodeComparaisonFlag = false;
  bool EtaFlag = false;

zhang's avatar
ZJ  
zhang committed
  //*************************************************************
  //=============================================================
  
  
nadolski's avatar
nadolski committed
//    // Chamber factory
//   if (ChamberFlag == false)
//     ChamberOff(); // no chamber limitations
//   else if (U20ChamberFlag == true)
//     DefineChNoU20;
  PrintCh();
  LoadApers("Apertures.dat", 1.0, 1.0);
zhang's avatar
ZJ  
zhang committed

nadolski's avatar
nadolski committed
  // 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);
zhang's avatar
ZJ  
zhang committed
  }
nadolski's avatar
nadolski committed

  // compute chromaticities by tracking (should be the same as by DA)
  if (ChromTracFlag == true){
    GetChromTrac(2, 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){
    fprintf(stdout, "\n Setting Multipoles \n");
    //Multipole_new(); /* for thick sextupoles */
    Multipole(); /* for thin 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
zhang's avatar
ZJ  
zhang committed
  if (FmapFlag == true){
    if (ChamberFlag == true ){
      if (ExperimentFMAFlag == true)
nadolski's avatar
nadolski committed
         fmap(40,12,258,-20e-3,5e-3,0.0,true); // for experimental
zhang's avatar
ZJ  
zhang committed
      if (DetailedFMAFlag == true)
nadolski's avatar
nadolski committed
        fmap(100,50,1026,20e-3,5e-3,0.0,true);
zhang's avatar
ZJ  
zhang committed
      }
      else{
        if (ExperimentFMAFlag == true)
nadolski's avatar
nadolski committed
          fmap(40,12,258,-32e-3,5e-3,0.0,true);
zhang's avatar
ZJ  
zhang committed
        if (DetailedFMAFlag == true)
nadolski's avatar
nadolski committed
          fmap(200,100,1026,32e-3,7e-3,0.0,true);
zhang's avatar
ZJ  
zhang committed
      }
  }
  
nadolski's avatar
nadolski committed
  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);
  }

zhang's avatar
ZJ  
zhang committed
  if (MomentumAccFlag == true){
nadolski's avatar
nadolski committed
    MomentumAcceptance(1L, 108L, 0.01, 0.05, 100L, -0.01, -0.05, 100L);
  }

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

zhang's avatar
ZJ  
zhang committed
  }
nadolski's avatar
nadolski committed
 
//  if (SigmaFlag == true){printsigma();
//  }
//  

  // induced amplitude 
  if (InducedAmplitudeFlag == true){
      InducedAmplitude(193L);  
  }
  
 if (EtaFlag == true){
  // compute cod and twiss parameters for different energy offsets    
    for (double 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);
    }
}    
  
  return 0;
zhang's avatar
ZJ  
zhang committed
  
  
  
  //---------------------------------------------------------
  // test region
  
 //  trace = true;
 // Get_Disp_dp();
 
  //InducedAmplitude(20);
  
 // Hcofonction(2, 0.001);
  
// fmapdp(20L, 21L, 1026L, 25e-3, 0.06, 0.3e-3, true);
 // return 1;
 
  //--------------------------------------------------------
  
  
  
  
  //*********************************************************************************
  //---------------------------------------------------------------------------------------------------------------------------------
  
  // 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
    
zhang's avatar
zhang committed
    
    //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
   
zhang's avatar
ZJ  
zhang committed
      
    // delicated for max4 lattice
zhang's avatar
zhang committed
    //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);
    
  }
zhang's avatar
ZJ  
zhang committed



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

  // Call nsls-ii_lib.cc
  // tune shift with amplitude 
nadolski's avatar
nadolski committed
  double delta = 0;
zhang's avatar
zhang committed
  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);
zhang's avatar
ZJ  
zhang committed
//    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)
zhang's avatar
zhang committed
  } 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);
  }
  
  
zhang's avatar
ZJ  
zhang committed

  
  
  
zhang's avatar
zhang committed
  
  //
  // 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);
    }
  
  }
}