Skip to content

GitLab

Commit 38cd9dc9 authored by adelmann's avatar adelmann 🎗
Browse files

Add first version of the matched distribution generator. 

Not yet in operation and some of the headerfiles need to be
cleaned up. Need to remove physics.h.

At the moment the matched distribution generator is not enabled!
parent 445f8cec
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 3816 additions and 15 deletions
.gitattributes
View file @ 38cd9dc9
......@@ -695,12 +695,22 @@ src/BasicActions/What.h -text
src/CMakeLists.txt -text
src/Copyright.readme -text
src/Distribution/CMakeLists.txt -text
src/Distribution/ClosedOrbitFinder.h -text
src/Distribution/Distribution.cpp -text
src/Distribution/Distribution.h -text
src/Distribution/LaserProfile.cpp -text
src/Distribution/LaserProfile.h -text
src/Distribution/MagneticField.h -text
src/Distribution/MapGenerator.h -text
src/Distribution/SigmaGenerator.h -text
src/Distribution/error.h -text
src/Distribution/halton1d_sequence.cpp -text
src/Distribution/halton1d_sequence.hh -text
src/Distribution/math.h -text
src/Distribution/matrix_vector_operation.h -text
src/Distribution/physical_error.h -text
src/Distribution/physics.h -text
src/Distribution/rdm.h -text
src/Editor/CMakeLists.txt -text
src/Editor/Edit.cpp -text
src/Editor/Edit.h -text
......
src/Distribution/ClosedOrbitFinder.h 0 → 100644
View file @ 38cd9dc9
This diff is collapsed.
src/Distribution/Distribution.cpp
View file @ 38cd9dc9
......@@ -12,7 +12,7 @@
// ------------------------------------------------------------------------
#include "Distribution/Distribution.h"
#include "Distribution/SigmaGenerator.h"
#include <cmath>
#include <cfloat>
#include <iomanip>
......@@ -35,6 +35,7 @@
#include "Algorithms/PartBinsCyc.h"
#include "BasicActions/Option.h"
#include "Distribution/LaserProfile.h"
#include "Elements/OpalBeamline.h"
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_randist.h>
......@@ -145,6 +146,13 @@ namespace AttributesT
VW,
SURFMATERIAL, // Add material type, currently 0 for copper
// and 1 for stainless steel.
EX, // below is for the matched distribution
EY,
ET,
FMAPFN,
HN,
FMSYM,
INTSTEPS,
SIZE
};
}
......@@ -222,8 +230,10 @@ Distribution::Distribution():
paraFNVYZe_m(0.0),
secondaryFlag_m(0),
ppVw_m(0.0),
vVThermal_m(0.0) {
vVThermal_m(0.0),
I_m(0.0),
E_m(0.0)
{
SetAttributes();
Distribution *defaultDistribution = clone("UNNAMED_Distribution");
......@@ -314,7 +324,10 @@ Distribution::Distribution(const std::string &name, Distribution *parent):
tFall_m(parent->tFall_m),
sigmaRise_m(parent->sigmaRise_m),
sigmaFall_m(parent->sigmaFall_m),
cutoff_m(parent->cutoff_m){
cutoff_m(parent->cutoff_m),
I_m(parent->I_m),
E_m(parent->E_m)
{
}
Distribution::~Distribution() {
......@@ -343,6 +356,16 @@ Distribution::~Distribution() {
}
}
/*
void Distribution::printSigma(SigmaGenerator<double,unsigned int>::matrix_type& M, Inform& out) {
for(int i=0; i<M.size1(); ++i) {
for(int j=0; j<M.size2(); ++j) {
*gmsg << M(i,j) << " ";
}
*gmsg << endl;
}
}
*/
/**
* At the moment only write the header into the file dist.dat
......@@ -391,6 +414,9 @@ void Distribution::Create(size_t &numberOfParticles, double massIneV) {
switch (distrTypeT_m) {
case DistrTypeT::MATCHEDGAUSS:
CreateMatchedGaussDistribution(numberOfParticles, massIneV);
break;
case DistrTypeT::FROMFILE:
CreateDistributionFromFile(numberOfParticles, massIneV);
break;
......@@ -1668,6 +1694,56 @@ void Distribution::CreateDistributionFromFile(size_t numberOfParticles, double m
inputFile.close();
}
void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, double massIneV) {
/// ADA
*gmsg << "----------------------------------------------------" << endl;
*gmsg << "About to find closed orbit and matched distribution " << endl;
*gmsg << "I= " << I_m*1E3 << " (mA) E= " << E_m*1E-6 << " (MeV)" << endl;
*gmsg << "EX= " << Attributes::getReal(itsAttr[AttributesT::EX])
<< " EY= " << Attributes::getReal(itsAttr[AttributesT::EY])
<< " ET= " << Attributes::getReal(itsAttr[AttributesT::ET])
<< " FMAPFN " << Attributes::getString(itsAttr[AttributesT::FMAPFN])
<< " FMSYM= " << Attributes::getReal(itsAttr[AttributesT::FMSYM]) << endl
<< " HN = " << Attributes::getReal(itsAttr[AttributesT::HN])
<< " INTSTEPS = " << Attributes::getReal(itsAttr[AttributesT::INTSTEPS]) << endl;
*gmsg << "----------------------------------------------------" << endl;
double wo = 8.44167E6*2.0*Physics::pi;
double m = 1;
SigmaGenerator<double,unsigned int> siggen(I_m*1E3,
Attributes::getReal(itsAttr[AttributesT::EX])*1E6,
Attributes::getReal(itsAttr[AttributesT::EY])*1E6,
Attributes::getReal(itsAttr[AttributesT::EX])*1E6,
wo,
E_m*1E6,
Attributes::getReal(itsAttr[AttributesT::HN]),m,72,590,
Attributes::getReal(itsAttr[AttributesT::FMSYM]),
Attributes::getReal(itsAttr[AttributesT::INTSTEPS]));
/*
if(siggen.match(1e-8,100,1000,7)) {
*gmsg << "Converged: Sigma-Matrix for " << E_m*1E6 << " MeV" << endl;
for(int i=0; i<siggen.getSigma().size1(); ++i) {
for(int j=0; j<siggen.getSigma().size2(); ++j) {
*gmsg << siggen.getSigma()(i,j) << " ";
}
*gmsg << endl;
}
} else {
*gmsg << "Not converged." << endl;
}
/// call Gauss
CreateDistributionGauss(numberOfParticles, massIneV);
*/
}
void Distribution::CreateDistributionGauss(size_t numberOfParticles, double massIneV) {
SetDistParametersGauss(massIneV);
......@@ -1723,18 +1799,31 @@ void Distribution::CreateBoundaryGeometry(PartBunch &beam, BoundaryGeometry &bg
void Distribution::CreateOpalCycl(PartBunch &beam,
size_t numberOfParticles,
double current, const Beamline &bl,
bool scan) {
/*
* When scan mode is true, we need to destroy particles except
* for the first pass.
*/
/*
* setup data for matched distribution generation
*/
E_m = (beam.getGamma()-1.0)*beam.getM();
I_m = current;
/*
Fixme:
// Beamline myLine = dynamic_cast<Beamline *>(bl.clone());
// FieldList cycl = bl.getElementByType("Cyclotron");
avrgpz_m = beam.getP()/beam.getM();
*/
/*
* When scan mode is true, we need to destroy particles except
* for the first pass.
*/
/*
Fixme:
avrgpz_m = beam.getP()/beam.getM();
*/
size_t numberOfPartToCreate = numberOfParticles;
if (beam.getTotalNum() != 0) {
scan_m = scan;
......@@ -3338,6 +3427,9 @@ void Distribution::PrintDist(Inform &os, size_t numberOfParticles) const {
case DistrTypeT::ASTRAFLATTOPTH:
PrintDistFlattop(os);
break;
case DistrTypeT::MATCHEDGAUSS:
PrintDistMatchedGauss(os);
break;
default:
INFOMSG("Distribution unknown." << endl;);
break;
......@@ -3445,6 +3537,12 @@ void Distribution::PrintDistFromFile(Inform &os) const {
<< Attributes::getString(itsAttr[AttributesT::FNAME]) << endl;
}
void Distribution::PrintDistMatchedGauss(Inform &os) const {
os << "* Distribution type: MATCHEDGAUSS" << endl;
}
void Distribution::PrintDistGauss(Inform &os) const {
os << "* Distribution type: GAUSS" << endl;
os << endl;
......@@ -3692,8 +3790,26 @@ void Distribution::SetAttributes() {
"SURFACEEMISSION, "
"SURFACERANDCREATE, "
"GUNGAUSSFLATTOPTH, "
"ASTRAFLATTOPTH",
"GAUSS");
"ASTRAFLATTOPTH, "
"GAUSS, "
"GAUSSMATCHED");
itsAttr[AttributesT::FMAPFN]
= Attributes::makeString("FMAPFN", "File for reading fieldmap used to create matched distibution ", "");
itsAttr[AttributesT::EX]
= Attributes::makeReal("EX", "Projected normalized emittance EX (m-rad), used to create matched distibution ", 1E-6);
itsAttr[AttributesT::EY]
= 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::HN]
= Attributes::makeReal("HN", "Harmonic number used to create matched distibution ", 1);
itsAttr[AttributesT::FMSYM]
= Attributes::makeReal("FMSYM", "Field map symmetry used to create matched distibution ", 1);
itsAttr[AttributesT::INTSTEPS]
= Attributes::makeReal("INTSTEPS", "Integration steps used to create matched distibution ", 1440);
itsAttr[AttributesT::FNAME]
= Attributes::makeString("FNAME", "File for reading in 6D particle "
......@@ -4056,6 +4172,8 @@ void Distribution::SetDistType() {
distrTypeT_m = DistrTypeT::SURFACEEMISSION;
else if(distT_m == "SURFACERANDCREATE")
distrTypeT_m = DistrTypeT::SURFACERANDCREATE;
else if(distT_m == "GAUSSMATCHED")
distrTypeT_m = DistrTypeT::MATCHEDGAUSS;
}
......
src/Distribution/Distribution.h
View file @ 38cd9dc9
......@@ -27,6 +27,7 @@
#include "Algorithms/PartData.h"
#include "Algorithms/Vektor.h"
#include "Beamlines/Beamline.h"
#include "Ippl.h"
......@@ -36,6 +37,7 @@
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_qrng.h>
class Beam;
class PartBunch;
class PartBins;
......@@ -53,7 +55,8 @@ namespace DistrTypeT
SURFACEEMISSION,
SURFACERANDCREATE,
GUNGAUSSFLATTOPTH,
ASTRAFLATTOPTH
ASTRAFLATTOPTH,
MATCHEDGAUSS
};
}
......@@ -93,6 +96,7 @@ public:
void CreateBoundaryGeometry(PartBunch &p, BoundaryGeometry &bg);
void CreateOpalCycl(PartBunch &beam,
size_t numberOfParticles,
double current, const Beamline &bl,
bool scan);
void CreateOpalE(Beam *beam,
std::vector<Distribution *> addedDistributions,
......@@ -197,6 +201,9 @@ private:
Distribution(const Distribution &);
void operator=(const Distribution &);
// void printSigma(SigmaGenerator<double,unsigned int>::matrix_type& M, Inform& out);
void AddDistributions();
void ApplyEmissionModel(double eZ, double &px, double &py, double &pz);
void ApplyEmissModelAstra(double &px, double &py, double &pz);
......@@ -215,6 +222,7 @@ private:
void CreateDistributionFlattop(size_t numberOfParticles, double massIneV);
void CreateDistributionFromFile(size_t numberOfParticles, double massIneV);
void CreateDistributionGauss(size_t numberOfParticles, double massIneV);
void CreateMatchedGaussDistribution(size_t numberOfParticles, double massIneV);
void DestroyBeam(PartBunch &beam);
void FillEBinHistogram();
void FillParticleBins();
......@@ -234,6 +242,7 @@ private:
void PrintDistFlattop(Inform &os) const;
void PrintDistFromFile(Inform &os) const;
void PrintDistGauss(Inform &os) const;
void PrintDistMatchedGauss(Inform &os) const;
void PrintDistSurfEmission(Inform &os) const;
void PrintDistSurfAndCreate(Inform &os) const;
void PrintEmissionModel(Inform &os) const;
......@@ -386,6 +395,20 @@ private:
/// of secondaries in Vaughan's model.
// AAA This is for the matched distribution
double ex_m;
double ey_m;
double et_m;
double I_m;
double E_m;
std::string bfieldfn_m; /// only temporarly
// Some legacy members that need to be cleaned up.
/// time binned distribution with thermal energy
......
src/Distribution/MagneticField.h 0 → 100644
View file @ 38cd9dc9
#ifndef MAGNETICFIELD_H
#define MAGNETICFIELD_H
#include <iostream>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <string>
#include <vector>
#define MALLOC(n,a) ((a*)malloc((n)*sizeof(a)))
/// Temporary class for comparison with reference program of Dr. Christian Baumgarten.
class MagneticField
{
public:
/// Returns the interpolated values of the magnetic field: B(r,theta), dB/dr, dB/dtheta
/*!
* @param *bint stores the interpolated magnetic field
* @param *brint stores the interpolated radial derivative of the magnetic field
* @param *btint stores the interpolated "theta" derivative of the magnetic field
* @param th is the angle theta where to interpolate
* @param nr is the number of radial steps
* @param nth is the number of angle steps per sector
* @param r is radial position where to interpolate
* @param r0 is the minimal extent of eh cyclotron
* @param dr is the radial step size
* @param **b stores the data of the magnetic field (read in by ReadSectorMap())
*/
static int interpolate(double*, double*, double*, double, int, int, double, double, double, float**);
/// Returns the interpolated values of the magnetic field: B(r,theta), dB/dr, dB/dtheta
/*!
* @param *bint stores the interpolated magnetic field
* @param *brint stores the interpolated radial derivative of the magnetic field
* @param *btint stores the interpolated "theta" derivative of the magnetic field
* @param ith is the angle theta where to interpolate
* @param nr is the number of radial steps
* @param nth is the number of angle steps per sector
* @param r is radial position where to interpolate
* @param r0 is the minimal extent of eh cyclotron
* @param dr is the radial step size
* @param **b stores the data of the magnetic field (read in by ReadSectorMap())
* @param btflag true if dB/dtheta should be computed, else false
*/
static int interpolate1(double*,double*,double*,int,int,int,double,double, double, float**,int);
/// Reads the magnetic field data from a file
/*!
* @param **b stores the read data
* @param nr is the number of radial steps
* @param nth is the number of angle steps per sector
* @param nsc is a symmetry parameter
* @param fname is the filename
* @param sum is 0
*/
static int ReadSectorMap(float**, int, int, int, std::string, double);
static void MakeNFoldSymmetric(float**, int, int, int, int);
static float** malloc2df(int, int);
private:
// ipolmethod==0 => lineare Interpolation
static int ipolmethod;
static void intcf(double, std::vector<double>&);
static void _intcf(double, std::vector<double>&, int);
static double norm360(double);
};
int MagneticField::interpolate(double* bint,double* brint,double* btint,double th,int nr,int nth,
double r,double r0, double dr, float** b) {
double dth,frb,fth;
std::vector<double> cofr(8);
std::vector<double> cofth(8);
int k0,k1,nrb,i,i1,j,ith,j1,xth,ctrl=0;
*bint=0.0;
*brint=0.0;
*btint=0.0;
// Map "th" into interval [0..360[:
th=norm360(th);
// Map radius "r" onto radial grid:
if (r<r0) return 0;
frb=(r-r0)/dr;
// Split radius in fractional and integer value for interpolation:
nrb=(int)frb;
frb-=(double)nrb;
dth=360./(double)nth;
fth=th/dth;
ith=(int)fth;
fth-=(double)ith;
// r-interpolation for b-field;
ctrl=0;
k0=0;
k1=4;
if (nrb+1>nr) return 2;
if (nrb+1==nr) {
/* Extrapolation: */
ctrl=2;
k0=0;
k1=3;
nrb--;
} else
if (nrb+2==nr) {
ctrl=1;
k0=0;
k1=3;
}
// Calculate coefficients for Laplace-Interpolation:
_intcf(frb,cofr,ctrl);
intcf(fth,cofth);
// Sum 'em up:
for (i=k0;i<k1;i++) {
xth=0;
i1=std::abs(nrb+i-1);
if (i1<0) {
i1=std::abs(i1);
xth=nth/2;
}
for (j=0;j<4;j++) {
j1=(ith+j-1+xth) % nth;
if (j1<0) j1+=nth;
*bint+=cofr[i]*cofth[j]*b[j1][i1]; // problem with b, since size not 1440
*brint+=cofr[i+4]*cofth[j]*b[j1][i1];
*btint+=cofr[i]*cofth[j+4]*b[j1][i1];
}
}
*btint*=((double)nth)/(2.0*M_PI);
*brint/=dr;
return 0;
}
int MagneticField::interpolate1(double *bint,double *brint,double *btint,int ith,int nr,int nth,
double r,double r0, double dr, float **b,int btflag) {
double frb,dth;
std::vector<double> cof(8);
int nrb,i,i0,i1,ith1,ith2,ctrl;
// Map "ith" into interval [0..nth-1]:
while (ith>=nth) ith-=nth;
while (ith<0) ith+=nth;
// Map radius "r" onto radial grid:
frb=(r-r0)/dr-1.0;
// Split radius in fractional and integer value for interpolation:
nrb=(int)frb;
frb-=(double)nrb;
ctrl=0;
i0=0;
i1=4;
if (nrb<0) {
ctrl=-1;
i0=1;
i1=4;
}
if (nrb+3==nr) {
ctrl=1;
i0=0;
i1=3;
}
if (nrb+2==nr) {
ctrl=2;
i0=0;
i1=3;
nrb--;
}
if (nrb+2>nr) {
return 2;
}
// r-interpolation for b-field;
// Calculate coefficients for Laplace-Interpolation:
_intcf(frb,cof,ctrl);
*bint=0.0;
*brint=0.0;
// Sum 'em up:
for (i=i0;i<i1;i++) {
*bint+=cof[i]*b[ith][nrb+i];
*brint+=cof[i+4]*b[ith][nrb+i];
}
(*brint)/=dr;
// Theta derivative?
if (btflag) {
*btint=0.0;
dth=2.0*M_PI/(double)nth;
ith1=ith-1;
if (ith1<0) ith1+=nth;
ith2=ith+1;
if (ith2>=nth) ith2-=nth;
for (i=i0;i<i1;i++) {
*btint+=cof[i]*(b[ith2][nrb+i]-b[ith1][nrb+i]);
}
(*btint)/=(2.0*dth);
}
return 0;
}
// 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];
// Return value is BOOL(success).
int MagneticField::ReadSectorMap(float** b,int nr,int nth,int nsc, std::string fname, double sum=0) {
double tmp;
int i,j,k,err=0;
FILE *f;
if (f=fopen(fname.c_str(),"r")) {
j=0;
while ((j<nr)&&(err==0)) { // j: distance (radius)
for (i=0;i<nth;i++) { // i: angles
err=(fscanf(f,"%lf",&tmp)==EOF);
if (!err) {
if (sum) {
for (k=0;k<nsc;k++)
b[i+nth*k][j]+=sum*tmp; // b is a (1440 x 8) x 141 matrix
} else {
for (k=0;k<nsc;k++)
b[i+nth*k][j]=tmp;
}
} else fprintf(stderr,"Error reading %s at (%d,%d)\n",fname.c_str(),j,i);
}
if (!err) j++;
}
fclose(f);
if (!err) fprintf(stderr,"Sector map %s read in.\n",fname.c_str());
return j;
} else fprintf(stderr,"Error opening %s.\n",fname.c_str());
return 0;
}
float** MagneticField::malloc2df(int n,int m) {
float **a;
int i;
a=MALLOC(n,float*);
for (i=0;i<n;i++) {
a[i]=MALLOC(m,float);
std::memset(a[i],0,m*sizeof(float));
}
return a;
}
int MagneticField::ipolmethod=1;
void MagneticField::intcf(double x, std::vector<double>& cof) {
// double 3-pt interpolation for general use;
// dimension cof(8);
double x2,x3;
if (ipolmethod) {
x2=x*x;
x3=x2*x;
cof[0]=x2-.5*(x3+x);
cof[1]=1.5*x3-2.5*x2+1.;
cof[2]=2.*x2-1.5*x3+.5*x;
cof[3]=.5*(x3-x2);
cof[4]=2.*x-1.5*x2-.5;
cof[5]=4.5*x2-5.*x;
cof[6]=4.*x-4.5*x2+.5;
cof[7]=1.5*x2-x;
} else {
cof[0]=0.0;
cof[1]=1.0-x;
cof[2]=x;
cof[3]=0.0;
cof[4]=0.0;
cof[5]=-1.0;
cof[6]=1.0;
cof[7]=0.0;
}