Commit 4ee77c59 authored by adelmann's avatar adelmann 🎗

New: Read field map information from header of the CARBONCYCL type, hence for the

matched distribution the CARBONCYCL type has to be used. For example the PSI Ring 
header would look like: 

1800.0
20.0
8
-4.0
1440
141
......

Important: the third line describes the number of sectors AND is not 
specified in the original definition of the CARBONCYCL type. THIS needs to be 
fixed as soon as the matched distribution is working properly.

At the moment the correct sigma matrix will be produced only with the ring590_bfld.dat:

 3.859  -0.6331 0     0	     0.4538 1.324

-0.6331  0.636  0     0      0.7276 -0.2605

 0	 0      6.256 1.197  0	     0

 0	 0      1.197 0.3845 0	     0

 0.4538  0.7276 0     0	     2.396   0.172

 1.324  -0.2605 0     0	     0.172   0.8693

Check at svn+ssh://savannah02.psi.ch/repos/opal/tests/trunk/RegressionTests/RingCyclotronMatched
parent 5fa7eac8
......@@ -25,7 +25,7 @@
#include "physics.h"
#include "physical_error.h"
#include "MagneticField.h" // ONLY FOR STAND-ALONE PROGRAM
#include "MagneticField.h"
#include <fstream>
......@@ -175,6 +175,10 @@ class ClosedOrbitFinder
/// Location of magnetic field
std::string fieldmap_m;
/// some proerties of the magnetic field map
int nr_m, nth_m, nsc_m;
double rmin_m, dr_m, dth_m;
/// Defines the stepper for integration of the ODE's
Stepper stepper_m;
......@@ -291,11 +295,15 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
* a' = a = acon
*/
// READ IN MAGNETIC FIELD: ONLY FOR STAND-ALONE PROGRAM
int nsc = 8, nr = 141, Nth = 1440, nth = 1440 / 8; value_type r0 = 1.8, dr = 0.02;
bmag_m = MagneticField::malloc2df(Nth,nr);
MagneticField::ReadSectorMap(bmag_m,nr,Nth,1,fieldmap_m,0.0);
MagneticField::MakeNFoldSymmetric(bmag_m,Nth,nr,nth,nsc);
// ada Injector 2 int nsc = 4, nr = 179, Nth = 144, nth = 144 / 4; value_type r0 = 0.24, dr = 0.02;
MagneticField::ReadHeader(&nr_m, &nth_m, &rmin_m, &dr_m, &dth_m, &nsc_m, fieldmap_m);
INFOMSG("Magnetic fiel map properties nr= " << nr_m << " nth= " << nth_m << " rmin= " << rmin_m << " dr= " << dr_m << " dth= " << dth_m << " nsc= " << nsc_m << endl);
bmag_m = MagneticField::malloc2df(nth_m,nr_m);
MagneticField::ReadSectorMap(bmag_m,nr_m,nth_m,1,fieldmap_m,0.0);
MagneticField::MakeNFoldSymmetric(bmag_m,nth_m,nr_m,nth_m,nsc_m);
value_type bint, brint, btint;
// velocity: beta = v/c = sqrt(1-1/(gamma*gamma))
......@@ -343,7 +351,7 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
invptheta = 1.0 / ptheta;
// intepolate values of magnetic field
MagneticField::interpolate(&bint,&brint,&btint,theta * 180 / M_PI,nr,Nth,y[0],r0,dr,bmag_m);
MagneticField::interpolate(&bint,&brint,&btint,theta * 180 / M_PI,nr_m,nth_m,y[0],rmin_m,dr_m,bmag_m);
bint *= invbcon;
brint *= invbcon;
......@@ -535,8 +543,6 @@ void ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeOrbitProperties()
* p. 6
*/
// READ IN MAGNETIC FIELD: ONLY FOR STAND-ALONE PROGRAM
int nsc = 8, nr = 141/*, Nth = 1440*/, nth = 1440 / 8; value_type r0 = 1.8, dr = 0.02;
value_type bint, brint, btint; // B, dB/dr, dB/dtheta
value_type invbcon = 1.0 / physics::bcon(wo_m);
......@@ -553,7 +559,7 @@ void ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeOrbitProperties()
for (size_type i = 0; i < N_m; ++i) {
// interpolate magnetic field
MagneticField::interpolate(&bint,&brint,&btint,theta * 180.0 / M_PI,nr,nth*nsc,r_m[i],r0,dr,bmag_m);
MagneticField::interpolate(&bint,&brint,&btint,theta * 180.0 / M_PI,nr_m,nth_m,r_m[i],rmin_m,dr_m,bmag_m);
bint *= invbcon;
brint *= invbcon;
btint *= invbcon;
......@@ -582,7 +588,7 @@ void ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeVerticalOscillati
vertOscDone_m = true;
// READ IN MAGNETIC FIELD: ONLY FOR STAND-ALONE PROGRAM
int /*nsc = 8,*/ nr = 141, Nth = 1440/*, nth = 1440/8*/; value_type r0 = 1.8, dr = 0.02;
value_type bint, brint, btint; // B, dB/dr, dB/dtheta
value_type en = E_m / physics::E0; // en = E/E0 = E/(mc^2) with kinetic energy E0
......@@ -607,7 +613,7 @@ void ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeVerticalOscillati
invptheta = 1.0 / ptheta;
// intepolate values of magnetic field
MagneticField::interpolate(&bint,&brint,&btint,theta * 180 / M_PI,nr,Nth,y[0],r0,dr,bmag_m);
MagneticField::interpolate(&bint,&brint,&btint,theta * 180 / M_PI,nr_m,nth_m,y[0],rmin_m,dr_m,bmag_m);
bint *= invbcon;
brint *= invbcon;
btint *= invbcon;
......
......@@ -46,6 +46,8 @@
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include "MagneticField.h"
extern Inform *gmsg;
#define DISTDBG1
......@@ -154,7 +156,6 @@ namespace AttributesT
ET,
LINE,
FMAPFN,
INTSTEPS,
RESIDUUM,
MAXSTEPSCO,
MAXSTEPSSI,
......@@ -1740,8 +1741,7 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
<< " ET= " << Attributes::getReal(itsAttr[AttributesT::ET])
<< " FMAPFN " << CyclotronElement->getFieldMapFN()
<< " FMSYM= " << CyclotronElement->getSymmetry()
<< " HN = " << CyclotronElement->getCyclHarm()
<< " INTSTEPS = " << Attributes::getReal(itsAttr[AttributesT::INTSTEPS]) << endl;
<< " HN = " << CyclotronElement->getCyclHarm() << endl;
*gmsg << "----------------------------------------------------" << endl;
const double wo = CyclotronElement->getRfFrequ()*1E6/CyclotronElement->getCyclHarm()*2.0*Physics::pi;
......@@ -1753,7 +1753,13 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
ERRORMSG("FMHIGHE of FMLOWE not set propperly" << endl);
exit(1);
}
int nr, nth, nsc;
double rmin, dr, dth;
MagneticField::ReadHeader(&nr, &nth, &rmin, &dr, &dth, &nsc, Attributes::getString(itsAttr[AttributesT::FMAPFN]));
SigmaGenerator<double,unsigned int> siggen(I_m,
Attributes::getReal(itsAttr[AttributesT::EX])*1E6,
Attributes::getReal(itsAttr[AttributesT::EY])*1E6,
......@@ -1764,8 +1770,8 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
massIneV*1E-6,
fmLowE,
fmHighE,
CyclotronElement->getSymmetry(),
Attributes::getReal(itsAttr[AttributesT::INTSTEPS]),
nsc,
nth,
Attributes::getString(itsAttr[AttributesT::FMAPFN]));
if(siggen.match(Attributes::getReal(itsAttr[AttributesT::RESIDUUM]),
......@@ -3897,8 +3903,6 @@ void Distribution::SetAttributes() {
= Attributes::makeReal("EY", "Projected normalized emittance EY (m-rad) used to create matched distibution ", 1E-6);
itsAttr[AttributesT::ET]
= Attributes::makeReal("ET", "Projected normalized emittance EY (m-rad) used to create matched distibution ", 1E-6);
itsAttr[AttributesT::INTSTEPS]
= Attributes::makeReal("INTSTEPS", "Integration steps used to create matched distibution ", 1440);
itsAttr[AttributesT::E2]
= Attributes::makeReal("E2", "If E2<Eb, we compute the tunes from the beams energy Eb to E2 with dE=0.25 MeV ", 0.0);
itsAttr[AttributesT::RESIDUUM]
......
......@@ -9,6 +9,11 @@
#include <cstring>
#include <string>
#include <vector>
#include <algorithm>
#include <iostream>
#include <fstream>
#include "Utilities/LogicalError.h"
#include <Ippl.h>
......@@ -59,6 +64,18 @@ public:
* @param sum is 0
*/
static int ReadSectorMap(float**, int, int, int, std::string, double);
/// Reads the header of magnetic field from a file
/*!
* @param nr is the number of radial steps
* @param nth is the number of angle steps per sector
* @param rmin minimal radial position (m)
* @param dr radial step (m)
* @parame dth angle per theta step (rad)
* @param fname is the filename
*/
static void ReadHeader(int*, int*, double*, double*, double*, int*, std::string);
static void MakeNFoldSymmetric(float**, int, int, int, int);
static float** malloc2df(int, int);
......@@ -200,6 +217,32 @@ int MagneticField::interpolate1(double *bint,double *brint,double *btint,int ith
return 0;
}
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(LogicalError(
"SectorMagneticFieldMap::IO::ReadLines",
"Failed to open file "+fname
));
}
double dthtmp;
fin >> *rmin;
*rmin /= 1000.0;
fin >> *dr;
*dr /= 1000.0;
fin >> *nsc;
fin >> dthtmp;
if (dthtmp<0.0)
*dth = -1.0/dthtmp;
else
*dth = dthtmp;
fin >> *nth;
fin >> *nr;
fin.close();
}
// Sum=0: Read Field Map with nr,nth entries from file "fname" into array "b"
// Sum=1: Add Field Map with nr,nth entries from file "fname" to array "b"
// b=b[nth*nsc][nr];
......@@ -211,7 +254,12 @@ int MagneticField::ReadSectorMap(float** b,int nr,int nth,int nsc, std::string f
f = fopen(fname.c_str(),"r");
if (f != NULL) {
j=0;
// overread the 6 header lines
for (int n=0; n<5; n++)
err=(fscanf(f,"%lf",&tmp)==EOF);
j=0;
while ((j<nr)&&(err==0)) { // j: distance (radius)
for (i=0;i<nth;i++) { // i: angles
err=(fscanf(f,"%lf",&tmp)==EOF);
......@@ -322,12 +370,14 @@ double MagneticField::norm360(double phi) {
void MagneticField::MakeNFoldSymmetric(float** bmag, int N, int nr, int nth, int nsc) {
double bav;
int nthpersec = nth/nsc;
for (int i=0;i<nr;i++) {
for (int j=0;j<nth;j++) {
for (int j=0;j<nthpersec;j++) {
bav=0.0;
for (int k=0;k<nsc;k++) bav+=bmag[(j+k*nth) % N][i];
for (int k=0;k<nsc;k++) bav+=bmag[(j+k*nthpersec) % N][i];
bav/=(double)nsc;
for (int k=0;k<nsc;k++) bmag[(j+k*nth) % N][i]=bav;
for (int k=0;k<nsc;k++) bmag[(j+k*nthpersec) % N][i]=bav;
}
}
}
......
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