Commit a8cd2ff9 authored by Andreas Adelmann's avatar Andreas Adelmann

some fixes

parent cfd53652
......@@ -65,11 +65,12 @@ class ClosedOrbitFinder
* @param nradial is the number of radial splits (fieldmap variable)
* @param dr is the radial step size, \f$ \left[\Delta r\right] = \si{m} \f$
* @param fieldmap is the location of the file that specifies the magnetic field
* @param guesss value of radius for closed orbit finder
* @param domain is a boolean (default: true). If "true" the closed orbit is computed over a single sector,
* otherwise over 2*pi.
*/
ClosedOrbitFinder(value_type, value_type, value_type, size_type, value_type, size_type, value_type, value_type, size_type,
value_type, size_type, size_type, value_type, const std::string&, bool = true);
value_type, size_type, size_type, value_type, const std::string&, value_type, bool = true);
/// Returns the inverse bending radius (size of container N+1)
container_type& getInverseBendingRadius();
......@@ -243,6 +244,9 @@ class ClosedOrbitFinder
/// ONLY FOR STAND-ALONE PROGRAM
float** bmag_m;
/// a guesss for the clo finder
value_type rguess_m;
/*!
* This quantity is defined in the paper "Transverse-Longitudinal Coupling by Space Charge in Cyclotrons"
......@@ -270,15 +274,16 @@ ClosedOrbitFinder<Value_type, Size_type, Stepper>::ClosedOrbitFinder(value_type
value_type Emin, value_type Emax, size_type nSector,
value_type rmin, size_type ntheta, size_type nradial,
value_type dr, const std::string& fieldmap,
value_type rguess,
bool domain)
: nxcross_m(0), nzcross_m(0), E_m(E), E0_m(E0), wo_m(wo), N_m(N), dtheta_m(Physics::two_pi/value_type(N)),
gamma_m(E/E0+1.0), ravg_m(0), phase_m(0), converged_m(false), Emin_m(Emin), Emax_m(Emax), nSector_m(nSector),
rmin_m(rmin), ntheta_m(ntheta), nradial_m(nradial), dr_m(dr), lastOrbitVal_m(0.0), lastMomentumVal_m(0.0),
vertOscDone_m(false), fieldmap_m(fieldmap), domain_m(domain), stepper_m()
vertOscDone_m(false), fieldmap_m(fieldmap), domain_m(domain), stepper_m(), rguess_m(rguess)
{
if (Emin_m > Emax_m || E_m < Emin_m || E > Emax_m)
throw std::domain_error("Error in ClosedOrbitFinder: Emin <= E <= Emax and Emin < Emax");
// if (Emin_m > Emax_m || E_m < Emin_m || E > Emax_m)
// throw std::domain_error("Error in ClosedOrbitFinder: Emin <= E <= Emax and Emin < Emax");
// velocity: beta = v/c = sqrt(1-1/(gamma*gamma))
if (gamma_m == 0)
......@@ -483,7 +488,6 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
// count number of crossings (excluding starting point --> idx>0)
nxcross_m += (idx > 0) * (y[4] * xold < 0);
xold = y[4];
++idx;
};
......@@ -541,12 +545,17 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
*
*/
// step size of energy
value_type dE;
if (Emin_m == Emax_m)
dE = 0.0;
else
dE = (E_m - Emin_m) / (Emax_m - Emin_m);
// iterate until suggested energy (start with minimum energy)
value_type E = Emin_m;
// step size of energy
value_type dE = (E_m - Emin_m) / (Emax_m - Emin_m);
// energy increase not more than 0.25
dE = (dE > 0.25) ? 0.25 : dE;
......@@ -561,12 +570,17 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
// set initial values for radius and radial momentum for lowest energy Emin
// orbit, [r] = m; Gordon, formula (20)
// radial momentum; Gordon, formula (20)
container_type init = {beta * acon, 0.0};
container_type init;
if (rguess_m < 0)
init = {beta * acon, 0.0};
else
init = {rguess_m, 0.0};
// store initial values for updating values for higher energies
container_type previous_init = {0.0, 0.0};
do {
// do {
// (re-)set inital values for r and pr
r_m[0] = init[0];
......@@ -610,7 +624,6 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
// compute amplitude of the error
error = std::sqrt(delta[0] * delta[0] + delta[1] * delta[1] * invgamma4) / r_m[0];
} while (error > accuracy && niterations++ < maxit);
// reset iteration counter
......@@ -633,7 +646,7 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
invgamma4 = 1.0 / (gamma2 * gamma2);
} while (E != E_m);
// } while (E != E_m);
/* store last entry, since it is needed in computeVerticalOscillations(), because we have to do the same
* number of integrations steps there.
......@@ -848,4 +861,4 @@ void ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeVerticalOscillati
phase_m *= nSector_m;
}
#endif
\ No newline at end of file
#endif
This diff is collapsed.
......@@ -9,6 +9,12 @@
#ifndef MAGNETICFIELD_H
#define MAGNETICFIELD_H
#include "Utilities/GeneralClassicException.h"
#define CHECK_CYC_FSCANF_EOF(arg) if(arg == EOF)\
throw GeneralClassicException("Cyclotron::getFieldFromFile",\
"fscanf returned EOF at " #arg);
#include <iostream>
#include <cmath>
......@@ -261,25 +267,22 @@ int MagneticField::interpolate1(double *bint,double *brint,double *btint,int ith
void MagneticField::ReadHeader(int *nr, int *nth, double *rmin, double *dr, double *dth, int *nsc, std::string fname){
std::ifstream fin(fname.c_str());
if (!fin || !fin.is_open())
throw std::runtime_error("Magnetic Field read error");
double dthtmp;
fin >> *rmin;
*rmin /= 1000.0;
fin >> *dr;
*dr /= 1000.0;
fin >> *nsc;
fin >> dthtmp;
if (dthtmp<0.0)
*dth = -1.0/dthtmp;
FILE *f = fopen(fname.c_str(),"r");
if (f != NULL) {
CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", rmin));
*rmin /= 1000.0;
CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", dr));
*dr /= 1000.0;
CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", nsc));
CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", dth));
if (*dth<0.0)
*dth = -1.0 / *dth;
CHECK_CYC_FSCANF_EOF(fscanf(f, "%d", nth));
CHECK_CYC_FSCANF_EOF(fscanf(f, "%d", nr));
fclose(f);
}
else
*dth = dthtmp;
fin >> *nth;
fin >> *nr;
fin.close();
throw GeneralClassicException("MagneticField::ReadHeader", "can not open " + fname);
}
// Sum=0: Read Field Map with nr,nth entries from file "fname" into array "b"
......@@ -290,19 +293,23 @@ int MagneticField::ReadSectorMap(float** b,int nr,int nth,int nsc, std::string f
double tmp;
int i,j,k,err=0;
FILE *f;
if ((nr <= 0) || (nth <= 0))
throw GeneralClassicException("MagneticField::ReadSectorMap", "nr or nth <= 0 ");
f = fopen(fname.c_str(),"r");
if (f != NULL) {
// overread the 6 header lines
for (int n=0; n<5; n++)
for (int n=0; n<6; n++)
err=(fscanf(f,"%lf",&tmp)==EOF);
int dread = 0;
j=0;
while ((j<nr)&&(err==0)) { // j: distance (radius)
for (i=0;i<nth;i++) { // i: angles
err=(fscanf(f,"%lf",&tmp)==EOF);
err=(fscanf(f,"%16lE",&tmp)==EOF);
if (!err) {
dread++;
if (sum) {
for (k=0;k<nsc;k++)
b[i+nth*k][j]+=sum*tmp; // b is a (1440 x 8) x 141 matrix
......@@ -316,8 +323,12 @@ int MagneticField::ReadSectorMap(float** b,int nr,int nth,int nsc, std::string f
if (!err) j++;
}
fclose(f);
return j;
} else std::cerr << "Error" << std::endl; //ERRORMSG("Can not open " << fname << endl);
if (dread != (nr*nth))
throw GeneralClassicException("MagneticField::ReadSectorMap", "Expected number of datapoint not matched");
else
return j;
} else
throw GeneralClassicException("MagneticField::ReadSectorMap", "Can not open" + fname);
return 0;
}
......
......@@ -109,9 +109,10 @@ class SigmaGenerator
* @param maxit is the maximum number of iterations (the program stops if within this value no stationary distribution was found)
* @param maxitOrbit is the maximum number of iterations for finding closed orbit
* @param angle defines the start of the sector (one can choose any angle between 0° and 360°)
* @param guess value of radius for closed orbit finder
* @param harmonic is a boolean. If "true" the harmonics are used instead of the closed orbit finder.
*/
bool match(value_type, size_type, size_type, value_type, bool);
bool match(value_type, size_type, size_type, value_type, value_type, bool);
/// Block diagonalizes the symplex part of the one turn transfer matrix
/** It computes the transformation matrix <b>R</b> and its inverse <b>invR</b>. It returns a vector containing the four eigenvalues
......@@ -340,7 +341,7 @@ SigmaGenerator<Value_type, Size_type>::SigmaGenerator(value_type I, value_type e
}
template<typename Value_type, typename Size_type>
bool SigmaGenerator<Value_type, Size_type>::match(value_type accuracy, size_type maxit, size_type maxitOrbit, value_type angle, bool harmonic) {
bool SigmaGenerator<Value_type, Size_type>::match(value_type accuracy, size_type maxit, size_type maxitOrbit, value_type angle, value_type rguess, bool harmonic) {
/* compute the equilibrium orbit for energy E_
* and get the the following properties:
* - inverse bending radius h
......@@ -364,7 +365,7 @@ bool SigmaGenerator<Value_type, Size_type>::match(value_type accuracy, size_type
ClosedOrbitFinder<value_type, size_type, boost::numeric::odeint::runge_kutta4<container_type> > cof(E_m, m_m, wo_m, N_m, accuracy,
maxitOrbit, Emin_m, Emax_m,
nSector_m, rmin_m, ntheta_m,
nradial_m, dr_m, fieldmap_m,
nradial_m, dr_m, fieldmap_m, rguess,
false);
// properties of one turn
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment