From f5410e990f3a337c5ae12ebd4ab8566603285995 Mon Sep 17 00:00:00 2001
From: zhang <zhang@9a6e40ed-f3a0-4838-9b4a-bf418f78e88d>
Date: Mon, 5 Jul 2010 15:00:46 +0000
Subject: [PATCH] ZJ new file
---
tracy/tracy/src/soleilcommon.cc | 1526 +++++++++++++++++++++++++++++++
1 file changed, 1526 insertions(+)
create mode 100644 tracy/tracy/src/soleilcommon.cc
diff --git a/tracy/tracy/src/soleilcommon.cc b/tracy/tracy/src/soleilcommon.cc
new file mode 100644
index 0000000..c80561b
--- /dev/null
+++ b/tracy/tracy/src/soleilcommon.cc
@@ -0,0 +1,1526 @@
+/*
+
+ Transferred from Tracy 2.7 soleilcommon.c and
+ modified based on Read_Lattice() in Tracy III physlib.cc
+*/
+
+
+/***************************************************************************
+ soleilcommon.c - description
+ -------------------
+ begin : Thu Oct 30 2003
+ copyright : (C) 2003 by nadolski
+ email : nadolski@synchrotron-soleil.fr
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ ***************************************************************************/
+
+//#include "tracy.h"
+
+//#include "soleilcommon.h"
+//#include "soleillib.h"
+//#include "datatyp.h"
+//#include "physlib.h"
+
+
+
+/****************************************************************************/
+/* void Read_Lattice(char *fic)
+
+ Purpose:
+ Read lattice file (fic.lat) and print statistics
+ Generate debugging info if problem in fic.lax
+ Initialize Tracy
+ Initialize the RING
+ Define the vacuum chamber
+ Compute Twiss parameters and chromaticities for on momentum particle
+ and save them in the file named linlat.lat
+
+ Input:
+ fic lattice file w/o its mandatory extension .lat
+
+ Output:
+ none
+
+ Return:
+ none
+
+ Global variables:
+ globval
+ S_SIZE
+
+ Specific functions:
+
+ t2init, Lattice_Read
+ Cell_Init()
+ DefineCh
+ Ring_GetTwiss
+ printglob
+
+ Comments:
+ 27/04/03 energy RF acceptance added set to 6%
+ 29/04/03 eps added for energy RF acceptance
+ 28/10/03 modified for transfer lines, filename added in output
+ 02/06/08 energy RF accpetance set to 1 just to avoid averflow during tracking
+ 22/06/10 add new globval. flag and modify new name of variables
+ from Tracy III Read_Lattice(), which is defined in physlib.cc
+
+****************************************************************************/
+void Read_Lattice(char *fic)
+{
+ char fic_maille[S_SIZE + 4] = "";
+ char fic_erreur[S_SIZE + 4] = "";
+ bool status;
+ bool chroma = true;
+ double dP = 0.0;
+ const double RFacceptance = 1.0; // maximum exusrtion during tracking
+ Vector2 beta, alpha, eta, etap;
+ Vector codvect;
+// double beta[2], alpha[2], eta[2], etap[2], codvect[6];
+ int i;
+
+ strcpy(fic_maille, fic);
+ strcpy(fic_erreur, fic);
+
+ /* generation automatique du nom du fichier maille et erreur */
+ strcat(fic_maille, ".lat");
+ strcat(fic_erreur, ".lax");
+
+ /* Initialisation de Tracy */
+
+ t2init();
+
+ /* open the lattice Input file */
+
+ if ((fi = fopen(fic_maille, "r")) == NULL)
+ {
+ fprintf(stdout,
+ "ReadLattice: Error while opening file %s \n",
+ fic_maille);
+ fprintf(stdout, "The lattice file has to end by .lat \n");
+ exit(1);
+ }
+
+ /* opens the lattice Output file */
+ if ((fo = fopen(fic_erreur, "w")) == NULL)
+ {
+ fprintf(stdout,
+ "ReadLattice: Error while opening file %s \n Access issue",
+ fic_erreur);
+ exit(1);
+ }
+
+ /* Reads lattice and set principle parameters
+ * Energy CODeps and energy offset
+ * print statistics
+ */
+ status = Lattice_Read(&fi, &fo);
+
+ if (status == false)
+ {
+ fprintf(stdout,
+ "Lattice_Read function has returned false\n");
+ fprintf(stdout, "See file %s \n", fic_erreur);
+ exit(1);
+ }
+ fprintf(stdout, "Lattice file: %s\n", fic_maille);
+
+ /* initializes cell structure: construction of the RING */
+ /* Creator of all the matrices for each element */
+ Cell_Init();
+
+
+ if (globval.RingType == 1)
+ { // for a ring
+ /* define x/y physical aperture */
+ //ChamberOff();
+ ChamberOn();
+
+ /* Defines global variables for Tracy code */
+
+ globval.H_exact = false; // Small Ring Hamiltonian
+ globval.quad_fringe = false; // quadrupole fringe fields on/off
+ globval.EPU = false; // Elliptically Polarizing Undulator
+ globval.IBS = false; /* diffusion on/off */
+
+ globval.MatMeth = false; /* matrix method */
+ globval.Cavity_on = false; /* Cavity on/off */
+ globval.radiation = false; /* radiation on/off */
+ globval.emittance = false; /* emittance on/off */
+ globval.pathlength = false; /* Path lengthening computation */
+ globval.CODimax = 40L; /* maximum number of iterations for COD algo */
+ globval.dPcommon = 1e-10; /* Common energy step for energy differentiation */
+ globval.delta_RF = RFacceptance;/* energy acceptance for SOLEIL */
+ globval.bpm = get_bpm_number(); /* number of bpm into the latice */
+ globval.hcorr = get_hcorr_number(); /* number of horizontal corrector into the lattice */
+ globval.vcorr = get_vcorr_number(); /* number of vertical corrector into the lattice */
+ globval.qt = get_qt_number(); /* number of skew quad into the lattice */
+
+
+/* int lastpos;
+ getcod(0.0, &lastpos);
+// findcod(0.0);
+ printcod();
+ exit(1);*/
+ /* Compute and get Twiss parameters */
+ Ring_GetTwiss(chroma = true, dP = 0.0);
+// GetChromTrac(2, 1026L, 1e-6, &ksix, &ksiz);
+// fprintf(stdout,"\n From tracking: ksix= % f ksiz= % f \n",ksix,ksiz);
+
+ Cell_SetdP(dP); /* added for correcting BUG if non convergence: compute on momentum linear matrices */
+ }
+ else
+ {
+ // for transfer lines
+ /* Initial settings : */
+ beta[0] = 8.1;
+ alpha[0] = 0.0;
+ beta[1] = 8.1;
+ alpha[1] = 0.0;
+ eta[0] = 0.0;
+ etap[0] = 0.0;
+ eta[1] = 0.0;
+ etap[1] = 0.0;
+
+ // for (i = 0; i < matdim; i++)
+ for (i = 0; i < ss_dim; i++) {
+ {
+ codvect[i] = 0.0;
+ globval.CODvect[i] = codvect[i];
+ }
+ dP = codvect[4];
+
+ /* Defines global variables for Tracy code */
+ globval.MatMeth = false; /* matrix method */
+ globval.Cavity_on = false; /* Cavity on/off */
+ globval.radiation = false; /* radiation on/off */
+ globval.emittance = false; /* emittance on/off */
+ globval.pathlength = false; /* Path lengthening computation */
+ globval.CODimax = 10L; /* maximum number of iterations for COD algo */
+ globval.dPcommon = 1e-10; /* Common energy step for energy differentiation */
+ globval.delta_RF = 0.10; /* 6% + epsilon energy acceptance for SOLEIL */
+ globval.bpm = get_bpm_number(); /* number of bpm into the latice */
+ globval.hcorr = get_hcorr_number(); /* number of horizontal corrector into the lattice */
+ globval.vcorr = get_vcorr_number(); /* number of vertical corrector into the lattice */
+ globval.qt = get_qt_number(); /* number of skew quad into the lattice */
+ globval.dPparticle = dP;
+
+ ChamberOff();
+ TransTwiss(alpha, beta, eta, etap, codvect);
+ }
+ }
+ printlatt(); /* print out lattice functions */
+ printglob();
+
+ bool a;
+ a =false;
+}
+
+
+
+
+/****************************************************************************/
+/* void ChamberOn(void)
+
+ Purpose:
+ Switch on the vacuum chamber
+
+ Input:
+ none
+
+ Output:
+ none
+
+ Return:
+ none
+
+ Global variables:
+ none
+
+ specific functions:
+ DefineCh
+
+ Comments:
+ none
+
+****************************************************************************/
+// void ChamberOn(void)
+// {
+// DefineCh();
+// status.chambre = true;
+// }
+
+/****************************************************************************/
+/* void ChamberOff(void)
+
+ Purpose:
+ Switch off the vacuum chamber
+ Still a check at 1 meter for avoiding numerical divergences
+
+ Input:
+ none
+
+ Output:
+ none
+
+ Return:
+ none
+
+ Global variables:
+ globval
+
+ specific functions:
+ none
+
+ Comments:
+ 16/06/03 Introduction of an asymmetrical vaccuum vessel
+
+// ****************************************************************************/
+// void ChamberOff(void)
+// {
+// int i;
+
+// for (i = 0; i <= globval.Cell_nLoc; i++)
+// {
+// globval.maxamplH[i][0] = -1.0;
+// globval.maxamplH[i][1] = 1.0;
+// globval.maxamplV[i][0] = -1.0;
+// globval.maxamplV[i][1] = 1.0;
+// }
+// status.chambre = false;
+// }
+
+/****************************************************************************/
+/* void PrintCh(void)
+
+ Purpose:
+ Print the vacuum chamber in the file chambre.out
+
+ Input:
+ none
+
+ Output:
+ none
+
+ Return:
+ none
+
+ Global variables:
+ globval
+
+ specific functions:
+ getglobv_, getelem
+
+ Comments:
+ 16/06/03 Introduction of an asymmetrical vaccuum vessel
+
+****************************************************************************/
+// void PrintCh(void)
+// {
+// FILE *f;
+// long i = 0, j = 0;
+// const char *fic = "chambre.out";
+// const long Taille = 5;
+// CellType Cell;
+// struct tm *newtime;
+
+// /* Get time and date */
+// newtime = GetTime();
+
+
+// /* open the lattice Input file */
+// if ((f = fopen(fic, "w")) == NULL)
+// {
+// fprintf(stderr, "PrintCh: Error while opening file %s \n",
+// fic);
+// exit(1);
+// }
+// fprintf(f, "# TRACY II v. SYNCHROTRON SOLEIL -- %s -- %s \n", fic,
+// asctime2(&newtime));
+// fprintf(f, "# i name s -xch(mm) +xch(mm) zch (mm)\n#\n");
+
+// for (i = 1; i <= globval.Cell_nLoc; i++)
+// {
+// getelem(i, &Cell);
+// fprintf(f, "%4ld ", i);
+// for (j = 0; j < Taille; j++)
+// {
+// if (Cell.Elem.PName[j] != ' ')
+// putc(Cell.Elem.PName[j], f);
+// else
+// putc(' ', f);
+// }
+// fprintf(f, " %6.2f %7.3f %7.3f %7.3f\n", Cell.S,
+// globval.maxamplH[i][0] * 1E3,
+// globval.maxamplH[i][1] * 1E3,
+// globval.maxamplV[i][1] * 1E3);
+// }
+
+// fclose(f);
+// }
+
+/****************************************************************************/
+/* void GetChromTrac(long Nb, long Nbtour, double emax, double *xix, double *xiz)
+
+ Purpose:
+ Computes chromaticities by tracking
+
+ Input:
+ Nb point number
+ Nbtour turn number
+ emax energy step
+
+ Output:
+ xix horizontal chromaticity
+ xiz vertical chromaticity
+
+ Return:
+ none
+
+ Global variables:
+ trace
+
+ Specific functions:
+ Trac_Simple, Get_NAFF
+
+ Comments:
+ 27/04/03 chromaticities are now output arguments
+
+****************************************************************************/
+// #define nterm 2
+// void GetChromTrac(long Nb, long Nbtour, double emax, double *xix, double *xiz)
+// {
+
+// double Tab[6][NTURN], fx[nterm], fz[nterm];
+// int nb_freq[2] = { 0, 0 }; /* frequency number to look for */
+// int i = 0;
+// boolean status = true;
+
+// double x = 1e-6, xp = 0.0, z = 1e-6, zp = 0.0;
+// double x0 = 1e-6, xp0 = 0.0, z0 = 1e-6, zp0 = 0.0;
+// double nux1, nux2, nuz1, nuz2;
+
+// /* initializations */
+// for (i = 0; i < nterm; i++)
+// {
+// fx[i] = 0.0;
+// fz[i] = 0.0;
+// }
+// /* end init */
+
+// /* Tracking for delta = emax and computing tunes */
+// x = x0;
+// xp = xp0;
+// z = z0;
+// zp = zp0;
+
+// Trac_Simple(x, xp, z, zp, emax, 0.0, Nbtour, Tab, &status);
+// Get_NAFF(nterm, Nbtour, Tab, fx, fz, nb_freq);
+
+// nux1 = (fabs (fx[0]) > 1e-8 ? fx[0] : fx[1]);
+// nuz1 = fz[0];
+
+// if (trace)
+// fprintf(stdout,
+// "\n Entering routine for chroma using tracking\n");
+// if (trace)
+// fprintf(stdout, "emax= % 10.6e nux1=% 10.6e nuz1= % 10.6e\n",
+// emax, nux1, nuz1);
+
+// /* Tracking for delta = -emax and computing tunes */
+// x = x0;
+// xp = xp0;
+// z = z0;
+// zp = zp0;
+
+// Trac_Simple(x, xp, z, zp, -emax, 0.0, Nbtour, Tab, &status);
+// Get_NAFF(nterm, Nbtour, Tab, fx, fz, nb_freq);
+
+// if (trace)
+// fprintf(stdout,
+// "nturn=%6ld x=% 10.5g xp=% 10.5g z=% 10.5g zp=% 10.5g delta=% 10.5g ctau=% 10.5g \n",
+// Nbtour, Tab[0][Nbtour - 1], Tab[1][Nbtour - 1],
+// Tab[2][Nbtour - 1], Tab[3][Nbtour - 1],
+// Tab[4][Nbtour - 1], Tab[5][Nbtour - 1]);
+
+// nux2 = (fabs(fx[0]) > 1e-8 ? fx[0] : fx[1]);
+// nuz2 = fz[0];
+
+// if (trace)
+// fprintf(stdout,
+// "emax= % 10.6e nux2= % 10.6e nuz2= % 10.6e\n", -emax,
+// nux2, nuz2);
+
+// /* Computing chromaticities */
+// *xix = (nux2 - nux1) * 0.5 / emax;
+// *xiz = (nuz2 - nuz1) * 0.5 / emax;
+
+// if (trace)
+// fprintf (stdout,
+// " Exiting routine for chroma using tracking\n\n");
+// }
+// #undef nterm
+
+// /****************************************************************************/
+// /* void GetTuneTrac(long Nbtour, double emax, double *nux, double *nuz)
+
+// Purpose:
+// Computes chromaticities by tracking
+
+// Input:
+// Nb point number
+// Nbtour turn number
+// emax energy step
+
+// Output:
+// none
+
+// Return:
+// none
+
+// Global variables:
+// trace
+
+// Specific functions:
+// Trac_Simple, Get_NAFF
+
+// Comments:
+// none
+
+// ****************************************************************************/
+// #define nterm 2
+// void GetTuneTrac(long Nbtour, double emax, double *nux, double *nuz)
+// {
+// double Tab[6][NTURN], fx[nterm], fz[nterm];
+// int nb_freq[2];
+// boolean status;
+
+// double x = 1e-6, xp = 0.0, z = 1e-6, zp = 0.0;
+
+// Trac_Simple(x, xp, z, zp, emax, 0.0, Nbtour, Tab, &status);
+// Get_NAFF(nterm, Nbtour, Tab, fx, fz, nb_freq);
+
+// *nux = (fabs (fx[0]) > 1e-8 ? fx[0] : fx[1]);
+// *nuz = fz[0];
+// }
+// #undef nterm
+
+/****************************************************************************/
+/* long get_bpm_number(void)
+
+ Purpose: called by Read_Lattice
+ Compute number of bpm in the lattice
+ The bpm name has to begin with 'bpm'
+
+ Input:
+ none
+
+ Output:
+ bpm bpm number found into the lattice
+
+ Return:
+ none
+
+ Global variables:
+ none
+
+ specific functions:
+ none
+
+ Comments:
+ none
+
+****************************************************************************/
+long get_bpm_number(void)
+{
+ int i;
+ long nbpm = 0L;
+ CellType Cell;
+
+ for (i = 0; i <= globval.Cell_nLoc; i++)
+ {
+ getelem(i, &Cell); /* get element */
+ if (strncmp(Cell.Elem.PName, "bpm", 3) == 0)
+ {
+ globval.bpm_list[nbpm] = i;
+ nbpm++;
+ }
+ }
+ return nbpm;
+}
+
+/****************************************************************************/
+/* long get_bpm_number(void)
+
+ Purpose: called by Read_Lattice
+ Compute number of corrector into the lattice
+ The corrector name has to begin with 'ch'
+
+ Input:
+ none
+
+ Output:
+ corrector number found in the lattice
+
+ Return:
+ none
+
+ Global variables:
+ none
+
+ specific functions:
+ none
+
+ Comments:
+ none
+
+****************************************************************************/
+long get_hcorr_number(void)
+{
+ int i;
+ long nb = 0L;
+ CellType Cell;
+
+ for (i = 0; i <= globval.Cell_nLoc; i++)
+ {
+ getelem(i, &Cell); /* get element */
+ if (strncmp(Cell.Elem.PName, "ch", 2) == 0)
+ {
+ globval.hcorr_list[nb] = i;
+ nb++;
+ }
+ }
+ return nb;
+}
+
+/****************************************************************************/
+/* long get_vcorr_number(void)
+
+ Purpose: called by Read_Lattice
+ Compute number of vertical corrector in the lattice
+ The corrector name has to begin with 'cv'
+
+ Input:
+ none
+
+ Output:
+ corrector number found into the lattice
+
+ Return:
+ none
+
+ Global variables:
+ none
+
+ specific functions:
+ none
+
+ Comments:
+ none
+
+****************************************************************************/
+long get_vcorr_number(void)
+{
+ int i;
+ long nb = 0L;
+ CellType Cell;
+
+ for (i = 0; i <= globval.Cell_nLoc; i++)
+ {
+ getelem(i, &Cell); /* get element */
+ if (strncmp(Cell.Elem.PName, "cv", 2) == 0)
+ {
+ globval.vcorr_list[nb] = i;
+ nb++;
+ }
+ }
+ return nb;
+}
+
+/****************************************************************************/
+/* long get_qt_number(void)
+
+ Purpose: called by Read_Lattice
+ Compute number of skew quadrupoles in the lattice
+ The skew quad name has to begin with 'qt'
+
+ Input:
+ none
+
+ Output:
+ skew quadrupole number found into the lattice
+
+ Return:
+ none
+
+ Global variables:
+ none
+
+ specific functions:
+ none
+
+ Comments:
+ none
+
+****************************************************************************/
+long get_qt_number(void)
+{
+ int i;
+ long nb = 0L;
+ CellType Cell;
+
+ for (i = 0; i <= globval.Cell_nLoc; i++)
+ {
+ getelem(i, &Cell); /* get element */
+ if (strncmp(Cell.Elem.PName, "qt", 2) == 0)
+ {
+ globval.qt_list[nb] = i;
+ Cell.Elem.M->Porder = 2L;
+ nb++;
+ }
+ }
+ return nb;
+}
+
+/****************************************************************************/
+/* void TransTwiss(double *alpha, double *beta, double *eta, double *etap, double *codvect)
+
+ Purpose: high level application
+ Calculate Twiss functions for a transport line
+
+ Input:
+ alpha alpha fonctions at the line entrance
+ beta beta fonctions at the line entrance
+ eta disperion fonctions at the line entrance
+ etap dipersion derivatives fonctions at the line entrance
+ codvect closed orbit fonctions at the line entrance
+
+ Output:
+ none
+
+ Return:
+ none
+
+ Global variables:
+
+
+ Specific functions:
+ TransTrace
+
+ Comments:
+ none
+
+****************************************************************************/
+// void TransTwiss(double *alpha, double *beta, double *eta, double *etap,
+// double *codvect)
+// {
+// TransTrace(0L, globval.Cell_nLoc, alpha, beta, eta, etap, codvect);
+// }
+
+
+
+// /****************************************************************************/
+// /* void findcodS(double dP)
+
+// Purpose:
+// Search for the closed orbit using a numerical method
+// Algo: Newton_Raphson method
+// Quadratic convergence
+// May need a guess starting point
+// Simple precision algorithm
+
+// Input:
+// dP energy offset
+
+// Output:
+// none
+
+// Return:
+// none
+
+// Global variables:
+// none
+
+// specific functions:
+// Newton_Raphson
+
+// Comments:
+// Method introduced because of bad convergence of DA for ID using RADIA maps
+
+// ****************************************************************************/
+// void findcodS(double dP)
+// {
+// float *vcod, *vector();
+// double x0[6];
+// const int ntrial = 40; // maximum number of trials for closed orbit
+// const float tolx = 1e-8; // numerical precision
+// int k;
+// int dim; // 4D or 6D tracking
+// long lastpos;
+// void free_vector();
+// CellType Cell;
+
+// getelem(1,&Cell);
+
+// vcod = vector(1, 6);
+
+// // starting point
+// for (k = 1; k <= 6; k++)
+// vcod[k] = (float) 0.0;
+
+// vcod[5] = (float) dP; // energy offset
+
+// if (globval.Cavity_on){
+// dim = 6; /* 6D tracking*/
+// fprintf(stderr,"Error looking for cod in 6D\n");
+// exit(1);
+// }
+// else{
+// dim = 4; /* 4D tracking */
+// vcod[1] = (float) Cell.Eta[0]*dP;
+// vcod[2] = (float) Cell.Etap[0]*dP;
+// vcod[3] = (float) Cell.Eta[1]*dP;
+// vcod[4] = (float) Cell.Etap[1]*dP;
+// }
+
+// Newton_RaphsonS(ntrial, vcod, dim, (float) tolx);
+
+// if (status.codflag == false)
+// fprintf(stdout, "Error No COD found\n");
+// if (trace)
+// {
+// for (k = 1; k <= 6; k++)
+// x0[k - 1] = (double) vcod[k];
+// fprintf(stdout,
+// "Before cod % .5e % .5e % .5e % .5e % .5e % .5e \n",
+// x0[0], x0[1], x0[2], x0[3], x0[4], x0[5]);
+// Cell_Pass(1L, globval.Cell_nLoc, x0, &lastpos);
+// fprintf(stdout,
+// "After cod % .5e % .5e % .5e % .5e % .5e % .5e \n",
+// x0[0], x0[1], x0[2], x0[3], x0[4], x0[5]);
+// Cell_Pass(1L, globval.Cell_nLoc, x0, &lastpos);
+// }
+// free_vector(vcod,1,6);
+// }
+
+// /****************************************************************************/
+// /* void findcod(double dP)
+
+// Purpose:
+// Search for the closed orbit using a numerical method
+// Algo: Newton_Raphson method
+// Quadratic convergence
+// May need a guess starting point
+// Simple precision algorithm
+// 4D
+// Starting point: linear closed orbit
+
+// 6D
+// Starting point: zero
+// if radiation on : x[5] is the synchroneous phase (equilibrium RF phase)
+// off: x[5] is zero
+
+// Input:
+// dP energy offset
+
+// Output:
+// none
+
+// Return:
+// none
+
+// Global variables:
+// none
+
+// specific functions:
+// Newton_Raphson
+
+// Comments:
+// Method introduced because of bad convergence of DA for ID using RADIA maps
+
+// ****************************************************************************/
+// void findcod(double dP)
+// {
+// double vcod[6];
+// const int ntrial = 40; // maximum number of trials for closed orbit
+// const double tolx = 1e-10; // numerical precision
+// int k, dim = 0;
+// long lastpos;
+// CellType Cell;
+
+// getelem(1,&Cell);
+
+// // initializations
+// for (k = 0; k <= 5; k++)
+// vcod[k] = 0.0;
+
+// if (globval.Cavity_on){
+// fprintf(stderr,"warning looking for cod in 6D\n");
+// dim = 6;
+// }
+// else{ // starting point linear closed orbit
+// dim = 4;
+// vcod[0] = Cell.Eta[0]*dP;
+// vcod[1] = Cell.Etap[0]*dP;
+// vcod[2] = Cell.Eta[1]*dP;
+// vcod[3] = Cell.Etap[1]*dP;
+// vcod[4] = dP; // energy offset
+// }
+
+// Newton_Raphson(dim, vcod, ntrial, tolx);
+
+// if (status.codflag == false)
+// fprintf(stdout, "Error No COD found\n");
+
+// CopyVec(6L, vcod, globval.CODvect); // save closed orbit at the ring entrance
+
+// if (trace)
+// {
+// fprintf(stdout,
+// "Before cod2 % .5e % .5e % .5e % .5e % .5e % .5e \n",
+// vcod[0], vcod[1], vcod[2], vcod[3], vcod[4], vcod[5]);
+// Cell_Pass(1L, globval.Cell_nLoc, vcod, &lastpos);
+// fprintf(stdout,
+// "After cod2 % .5e % .5e % .5e % .5e % .5e % .5e \n",
+// vcod[0], vcod[1], vcod[2], vcod[3], vcod[4], vcod[5]);
+// }
+// }
+// /****************************************************************************/
+// /* void computeFandJS(float *x, int n, float **fjac, float *fvect)
+
+// Purpose:
+// Simple precision algo
+// Tracks x over one turn. And computes the Jacobian matrix of the
+// transformation by numerical differentiation.
+// using forward difference formula : faster but less accurate
+// using symmetric difference formula
+
+// Input:
+// x vector for evaluation
+// n dimension 4 or 6
+
+// Output:
+// fvect transport of x over one turn
+// fjac Associated jacobian matrix
+
+// Return:
+// none
+
+// Global variables:
+// none
+
+// specific functions:
+// none
+
+// Comments:
+// none
+
+// ****************************************************************************/
+
+// void computeFandJS(float *x, int n, float **fjac, float *fvect)
+// {
+// int i, k;
+// long lastpos = 0L;
+// double x0[6], fx[6], fx1[6], fx2[6];
+// const double deps = 1e-8; //stepsize for numerical differentiation
+
+// for (i = 1; i <= 6; i++)
+// x0[i - 1] = (double) x[i];
+
+// Cell_Pass(1L, globval.Cell_nLoc, x0, &lastpos);
+
+// for (i = 1; i <= n; i++)
+// {
+// fvect[i] = (float) x0[i - 1];
+// fx[i - 1] = x0[i - 1];
+// }
+
+// // compute Jacobian matrix by numerical differentiation
+// for (k = 0; k < n; k++)
+// {
+// for (i = 1; i <= 6; i++)
+// x0[i - 1] = (double) x[i];
+// x0[k] += deps; // differential step in coordinate k
+
+// Cell_Pass(1L, globval.Cell_nLoc, x0, &lastpos); // tracking along the ring
+// for (i = 1; i <= 6; i++)
+// fx1[i - 1] = x0[i - 1];
+
+// for (i = 1; i <= 6; i++)
+// x0[i - 1] = (double) x[i];
+// x0[5] = 0.0;
+// x0[k] -= deps; // differential step in coordinate k
+
+// Cell_Pass(1L, globval.Cell_nLoc, x0, &lastpos); // tracking along the ring
+// for (i = 1; i <= 6; i++)
+// fx2[i - 1] = x0[i - 1];
+
+// for (i = 1; i <= n; i++) // symmetric difference formula
+// fjac[i][k + 1] =
+// (float) (0.5 * (fx1[i - 1] - fx2[i - 1]) / deps);
+// //~ for (i = 1; i <= n; i++) // forward difference formula
+// //~ fjac[i][k + 1] = (float) ((x0[i - 1] - fx[i - 1]) / deps);
+// }
+// }
+
+// /****************************************************************************/
+// /* void computeFand(int n, float *x, float **fjac, float *fvect)
+
+// Purpose:
+// Tracks x over one turn. And computes the Jacobian matrix of the
+// transformation by numerical differentiation.
+// using symmetric difference formula
+// double precision algorithm
+
+// Input:
+// x vector for evaluation
+
+// Output:
+// fvect transport of x over one turn
+// fjac Associated jacobian matrix
+
+// Return:
+// none
+
+// Global variables:
+// none
+
+// specific functions:
+// none
+
+// Comments:
+// none
+
+// ****************************************************************************/
+// void computeFandJ(int n, double *x, vector *fjac, double *fvect)
+// {
+// int i, k;
+// long lastpos = 0L;
+// double x0[6], fx1[6], fx2[6];
+// const double deps = 1e-8; //stepsize for numerical differentiation
+
+// CopyVec(6L, x, x0);
+
+// Cell_Pass(1L, globval.Cell_nLoc, x0, &lastpos);
+// CopyVec((long) n, x0, fvect);
+
+// // compute Jacobian matrix by numerical differentiation
+// for (k = 0; k < n; k++)
+// {
+// CopyVec(6L, x, x0);
+// x0[k] += deps; // differential step in coordinate k
+
+// Cell_Pass(1L, globval.Cell_nLoc, x0, &lastpos); // tracking along the ring
+// CopyVec(6L, x0, fx1);
+
+// CopyVec(6L, x, x0);
+// x0[k] -= deps; // differential step in coordinate k
+
+// Cell_Pass(1L, globval.Cell_nLoc, x0, &lastpos); // tracking along the ring
+// CopyVec(6L, x0, fx2);
+
+// for (i = 0; i < n; i++) // symmetric difference formula
+// fjac[i][k] = 0.5 * (fx1[i] - fx2[i]) / deps;
+// }
+// }
+
+// /****************************************************************************/
+// /* void Newton_RaphsonS(int ntrial,float x[],int n,float tolx, float tolf)
+
+// Purpose:
+// Newton_Rapson algorithm from Numerical Recipes
+// single precision algorithm
+// Robustess: quadratic convergence
+// Hint: for n-dimensional problem, the algo can be stuck on local minimum
+// In this case, it should be enough to provide a resonable starting
+// point.
+
+// Method:
+// look for closed orbit solution of f(x) = x
+// This problems is equivalent to finding the zero of g(x)= f(x) - x
+// g(x+h) ~= f(x) - x + (Jacobian(f) -Id) h + O(h*h)
+// Then at first order we solve h:
+// h = - inverse(Jacobian(f) -Id) * (f(x)-x)
+// the new guess is then xnew = x + h
+// By iteration, this converges quadratically.
+
+// The algo is stopped whenever |x -xnew| < tolx
+
+// f(x) is computes by tracking over one turn
+// Jacobian(f) is computed numerically by numerical differentiation
+// These two operations are provided by the function computeFandJ
+
+// Input:
+// ntrial number of iterations for closed zero search
+// n number of dimension 4 or 6
+// x intial guess for the closed orbit
+// tolx tolerance over the solution x
+// tolf tolerance over the evalution f(x)
+
+// Output:
+// x closed orbit
+
+// Return:
+// none
+
+// Global variables:
+// status
+
+// specific functions:
+// computeFandJS
+// ludcmp,lubksb
+
+// Comments:
+// none
+
+// ****************************************************************************/
+
+// #define FREERETURN {free_matrix(alpha,1,n,1,n);free_vector(bet,1,n); free_vector(fvect,1,n); free_ivector(indx,1,n);return;}
+
+// void Newton_RaphsonS(int ntrial, float x[], int n, float tolx)
+// {
+// int k, i, *indx, *ivector();
+// float errx, d, *bet, *fvect, **alpha, *vector(), **matrix();
+// void usrfun(), ludcmp(), lubksb(), free_ivector(), free_vector(),
+// free_matrix();
+
+// errx = (float) 0.0;
+// // NR arrays start from 1 and not 0 !!!
+// indx = ivector(1, n);
+// bet = vector(1, n);
+// fvect = vector(1, n);
+// alpha = matrix(1, n, 1, n);
+
+// for (k = 1; k <= ntrial; k++)
+// { // loop over number of iterations
+// computeFandJS(x, n, alpha, fvect); // supply function values at x in fvect and Jacobian matrix in fjac
+
+// // Jacobian -Id
+// for (i = 1; i <= n; i++)
+// alpha[i][i] -= (float) 1.0;
+// for (i = 1; i <= n; i++)
+// bet[i] = x[i] - fvect[i]; // right side of linear equation
+// // solve linear equations using LU decomposition using NR routines
+// ludcmp(alpha, n, indx, &d);
+// lubksb(alpha, n, indx, bet);
+// errx = (float) 0.0; // check root convergence
+// for (i = 1; i <= n; i++)
+// { // update solution
+// errx += fabs(bet[i]);
+// x[i] += bet[i];
+// }
+
+// if (trace)
+// fprintf(stdout,
+// "%02d: cod % .5e % .5e % .5e % .5e % .5e % .5e errx =% .5e\n",
+// k, x[1], x[2], x[3], x[4], x[5], x[6], errx);
+// if (errx <= tolx)
+// {
+// status.codflag = true;
+// FREERETURN}
+// }
+// // check whever closed orbit found out
+// if ((k >= ntrial) && (errx >= tolx * 100))
+// {
+// status.codflag = false;
+// FREERETURN}
+// }
+
+// #undef FREERETURN
+
+// /****************************************************************************/
+// /* int Newton_Raphson(int n, double x[], int ntrial, double tolx)
+
+// Purpose:
+// Newton_Rapson algorithm from Numerical Recipes
+// double precision algorithm
+// Robustess: quadratic convergence
+// Hint: for n-dimensional problem, the algo can be stuck on local minimum
+// In this case, it should be enough to provide a resonable starting
+// point.
+
+// Method:
+// look for closed orbit solution of f(x) = x
+// This problems is equivalent to finding the zero of g(x)= f(x) - x
+// g(x+h) ~= f(x) - x + (Jacobian(f) -Id) h + O(h*h)
+// Then at first order we solve h:
+// h = - inverse(Jacobian(f) -Id) * (f(x)-x)
+// the new guess is then xnew = x + h
+// By iteration, this converges quadratically.
+
+// The algo is stopped whenever |x -xnew| < tolx
+
+// f(x) is computes by tracking over one turn
+// Jacobian(f) is computed numerically by numerical differentiation
+// These two operations are provided by the function computeFandJ
+
+// Input:
+// ntrial number of iterations for closed zero search
+// x intial guess for the closed orbit
+// tolx tolerance over the solution x
+// tolf tolerance over the evalution f(x)
+
+// Output:
+// x closed orbit
+
+// Return:
+// none
+
+// Global variables:
+// status
+
+// specific functions:
+// computeFandJ
+// InvMat, LinTrans
+
+// Comments:
+// none
+
+// ****************************************************************************/
+// int Newton_Raphson (int n, double x[], int ntrial, double tolx)
+// {
+// int k, i;
+// double errx, bet[6], fvect[6];
+// vector alpha[6];
+
+// errx = 0.0;
+
+// for (k = 1; k <= ntrial; k++)
+// { // loop over number of iterations
+// computeFandJ(n,x, alpha, fvect); // supply function values at x in fvect and Jacobian matrix in fjac
+
+// // Jacobian - Id
+// for (i = 0; i < n; i++)
+// alpha[i][i] -= 1.0;
+// for (i = 0; i < n; i++)
+// bet[i] = x[i] - fvect[i]; // right side of linear equation
+// // inverse matrix using gauss jordan method from Tracy (from NR)
+// if (!InvMat((long) n,alpha)) fprintf(stdout,"Matrix non inversible ...\n");
+// LinTrans((long) n, alpha, bet); // bet = alpha*bet
+// errx = 0.0; // check root convergence
+// for (i = 0; i < n; i++)
+// { // update solution
+// errx += fabs(bet[i]);
+// x[i] += bet[i];
+// }
+
+// if (trace)
+// fprintf(stdout,
+// "%02d: cod2 % .5e % .5e % .5e % .5e % .5e % .5e errx =% .5e\n",
+// k, x[0], x[1], x[2], x[3], x[4], x[5], errx);
+// if (errx <= tolx)
+// {
+// status.codflag = true;
+// return 1;
+// }
+// }
+// // check whever closed orbit found out
+// if ((k >= ntrial) && (errx >= tolx))
+// {
+// status.codflag = false;
+// return 1;
+// }
+// return 0;
+// }
+
+
+// /****************************************************************************/
+// /* void Get_NAFF(int nterm, long ndata, double T[DIM][NTURN],
+// double *fx, double *fz, int nb_freq[2])
+
+// Purpose:
+// Compute quasiperiodic approximation of a phase space trajectory
+// using NAFF Algorithm ((c) Laskar, IMCCE)
+
+// Input:
+// nterm number of frequencies to look for
+// if not multiple of 6, truncated to lower value
+// ndata size of the data to analyse
+// T 6D vector to analyse
+
+// Output:
+// fx frequencies found in the H-plane
+// fz frequencies found in the V-plane
+// nb_freq number of frequencies found out in each plane
+
+// Return:
+// none
+
+// Global variables:
+// g_NAFVariable see modnaff.c
+// M_2_PI defined in math.h
+// trace ON or TRUE for debugging
+
+// Specific functions:
+// naf_initnaf, naf_cleannaf
+// naf_mftnaf
+
+// Comments:
+// none
+
+// ****************************************************************************/
+// /* Frequency Map Analysis */
+// #include "modnaff.h"
+// #include "complexeheader.h"
+// /* Analyse en Frequence */
+// double pi = M_PI;
+// void Get_NAFF(int nterm, long ndata, double T[DIM][NTURN],
+// double *fx, double *fz, int nb_freq[2])
+// {
+// /* Test whether ndata is divisible by 6 -- for NAFF -- */
+// /* Otherwise truncate ndata to lower value */
+// long r = 0; /* remainder of the euclidian division of ndata by 6 */
+// int i;
+
+// if ((r = ndata % 6) != 0) {
+// printf("Get_NAFF: Warning ndata = %ld, \n", ndata);
+// ndata -= r;
+// printf("New value for NAFF ndata = %ld \n", ndata);
+// }
+
+// g_NAFVariable.DTOUR = M_2_PI; /* size of one "cadran" */
+// g_NAFVariable.XH = M_2_PI; /* step */
+// g_NAFVariable.T0 = 0.0; /* time t0 */
+// g_NAFVariable.NTERM = nterm; /* max term to find */
+// g_NAFVariable.KTABS = ndata; /* number of data: must be a multiple of 6 */
+// g_NAFVariable.m_pListFen = NULL; /* no window */
+// g_NAFVariable.TFS = NULL; /* will contain frequency */
+// g_NAFVariable.ZAMP = NULL; /* will contain amplitude */
+// g_NAFVariable.ZTABS = NULL; /* will contain data to analyze */
+
+// /****************************************************/
+// /* internal use in naf */
+// g_NAFVariable.NERROR = 0;
+// g_NAFVariable.ICPLX = 1;
+// g_NAFVariable.IPRT = -1; /* 1 for diagnostics */
+// g_NAFVariable.NFPRT = stdout; /* NULL */
+// g_NAFVariable.NFS = 0;
+// g_NAFVariable.IW = 1;
+// g_NAFVariable.ISEC = 1;
+// g_NAFVariable.EPSM = 0;
+// g_NAFVariable.UNIANG = 0;
+// g_NAFVariable.FREFON = 0;
+// g_NAFVariable.ZALP = NULL;
+// g_NAFVariable.m_iNbLineToIgnore = 1; /* unused */
+// g_NAFVariable.m_dneps = 1.e10;
+// g_NAFVariable.m_bFSTAB = FALSE; /* unused */
+// /* end of interl use in naf */
+// /****************************************************/
+
+// /* NAFF initialization */
+// naf_initnaf();
+
+// /**********************/
+// /* Analyse in H-plane */
+// /**********************/
+
+// /* fills up complexe vector for NAFF analysis */
+// for(i = 0; i < ndata; i++) {
+// g_NAFVariable.ZTABS[i].reel = T[0][i]; /* x */
+// g_NAFVariable.ZTABS[i].imag = T[1][i]; /* xp */
+// }
+
+// /* Get out the mean value */
+// naf_smoy(g_NAFVariable.ZTABS);
+
+// naf_prtabs(g_NAFVariable.KTABS,g_NAFVariable.ZTABS, 20);
+// // trace=on;
+// naf_mftnaf(nterm,fabs(g_NAFVariable.FREFON)/g_NAFVariable.m_dneps);
+
+// /* fill up H-frequency vector */
+// for (i = 1; i <= g_NAFVariable.NFS; i++) {
+// fx[i-1] = g_NAFVariable.TFS[i];
+// }
+
+// nb_freq[0] = g_NAFVariable.NFS; /* nb of frequencies found out by NAFF */
+
+// if (trace) /* print out results */
+// {
+// printf("(x,x') phase space: NFS=%d\n",g_NAFVariable.NFS);
+// for (i = 1; i <= g_NAFVariable.NFS; i++) {
+// printf("AMPL=%15.8E+i*%15.8E abs(AMPL)=%15.8E arg(AMPL)=%15.8E FREQ=%15.8E\n",
+// g_NAFVariable.ZAMP[i].reel,g_NAFVariable.ZAMP[i].imag,
+// i_compl_module(g_NAFVariable.ZAMP[i]),
+// i_compl_angle(g_NAFVariable.ZAMP[i]),
+// g_NAFVariable.TFS[i]);
+// }
+// }
+
+// /**********************/
+// /* Analyse in V-plane */
+// /**********************/
+
+// /* fill up complexe vector for NAFF analysis */
+// for (i = 0; i < ndata; i++) {
+// g_NAFVariable.ZTABS[i].reel = T[2][i]; /* z */
+// g_NAFVariable.ZTABS[i].imag = T[3][i]; /*zp */
+// }
+
+// naf_mftnaf(nterm,fabs(g_NAFVariable.FREFON)/g_NAFVariable.m_dneps);
+
+// /* fills up V-frequency vector */
+// for (i = 1; i <= g_NAFVariable.NFS; i++) {
+// fz[i-1] = g_NAFVariable.TFS[i];
+// }
+
+// nb_freq[1] = g_NAFVariable.NFS; /* nb of frequencies found out by NAFF */
+
+// if (trace) /* print out results */
+// {
+// printf("(z,z') phase space: NFS=%d\n",g_NAFVariable.NFS);
+// for (i = 1; i <= g_NAFVariable.NFS; i++) {
+// printf("AMPL=%15.8E+i*%15.8E abs(AMPL)=%15.8E arg(AMPL)=%15.8E FREQ=%15.8E\n",
+// g_NAFVariable.ZAMP[i].reel,g_NAFVariable.ZAMP[i].imag,
+// i_compl_module(g_NAFVariable.ZAMP[i]),
+// i_compl_angle(g_NAFVariable.ZAMP[i]),
+// g_NAFVariable.TFS[i]);
+// }
+// }
+
+// /* free out memory used by NAFF */
+// naf_cleannaf();
+// }
+
+// /****************************************************************************/
+// /* void Get_Tabshift(double Tab[DIM][NTURN], double Tab0[DIM][NTURN],
+// long nbturn, long nshift)
+
+// Purpose: used by fmap
+// Store in Tab0 values of Tab shifted by nshift turn
+
+// Input:
+// Tab tracking tabular
+// nshift nb of turns to shift
+// nbturn nb of turns
+
+// Output:
+// Tab0 tracking tabular
+
+// Return:
+// none
+
+// Global variables:
+// none
+
+// Specific functions:
+// none
+
+// Comments:
+
+
+// ****************************************************************************/
+// void Get_Tabshift(double Tab[DIM][NTURN], double Tab0[DIM][NTURN], long nbturn, long nshift)
+// {
+// long k = 0L, k1 = 0L;
+
+// for (k = 0L; k < nbturn ; k++){
+// k1 = k + nshift;
+// Tab0[0][k] = Tab[0][k1];
+// Tab0[1][k] = Tab[1][k1];
+// Tab0[2][k] = Tab[2][k1];
+// Tab0[3][k] = Tab[3][k1];
+// }
+
+// }
+
+// /****************************************************************************/
+// /* void Get_freq(double *fx, double *fz, double *nux, double *nuz)
+
+// Purpose: used by fmap
+// Looks for tunes from frequency vectors output from NAFF
+
+// Input:
+// fx vector of horizontal frequencies
+// fz vector of verticcal frequencies
+
+// Output:
+// nux H tune
+// nuz V tune
+
+// Return:
+// none
+
+// Global variables:
+// none
+
+// Specific functions:
+// none
+
+// Comments:
+
+
+// ****************************************************************************/
+// void Get_freq(double *fx, double *fz, double *nux, double *nuz)
+// {
+// const double eps0 = 1e-4;
+// const double eps1 = 1e-6;
+
+// // case of nux
+// if (fabs(fx[0]) < eps0){
+// *nux = fx[1];
+// }
+// else {
+// *nux = fx[0];
+// }
+
+// // case of nuz
+// if (fabs(fz[0]) < eps0) {
+// if (fabs(fabs(fz[1]) - fabs(*nux)) < eps1) {
+// if (fabs(fabs(fz[2]) - fabs(*nux)) < eps1) {
+// *nuz = fz[3];
+// }
+// else *nuz = fz[2];
+// }
+// else *nuz = fz[1];
+// }
+// else{
+// if (fabs(fabs(fz[0]) - fabs(*nux)) < eps1) {
+// if (fabs(fabs(fz[1]) - fabs(*nux)) < eps1) {
+// *nuz = fz[2];
+// }
+// else *nuz = fz[1];
+// }
+// else *nuz = fz[0];
+// }
+// *nuz = fabs(*nuz); *nux = fabs(*nux);
+// }
+
+
+// /* DO NOT REMOVE FOR f2c ...*/
+// /* To be fixed later .. */
+// void MAIN__(void)
+// {
+
+// }
+
+/****************************************************************************/
+/*
+
+ Purpose:
+
+ Input:
+ none
+
+ Output:
+
+
+ Return:
+ none
+
+ Global variables:
+ none
+
+ specific functions:
+ none
+
+ Comments:
+ none
+
+****************************************************************************/
--
GitLab