Skip to content
Snippets Groups Projects
Commit f5410e99 authored by zhang's avatar zhang
Browse files

ZJ new file

parent d5d78b9b
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
tracy/tracy/src/soleilcommon.cc 0 → 100644
+ 1526
− 0
View file @ f5410e99
/*
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
****************************************************************************/
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

The place to be for your source code