Skip to content

GitLab

Commit a21ca6e8 authored by frey_m's avatar frey_m
Browse files

Merge branch '285-matched-distribution-not-matched' into 'master' 

Resolve "Matched Distribution: Not matched?"

See merge request !393
parents cccfc5ee a523040e
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 1587 additions and 1741 deletions
src/Distribution/CMakeLists.txt
View file @ a21ca6e8
set (_SRCS
Distribution.cpp
LaserProfile.cpp
RealDiracMatrix.cpp
SigmaGenerator.cpp
)
include_directories (
......@@ -12,10 +14,10 @@ add_opal_sources(${_SRCS})
set (HDRS
ClosedOrbitFinder.h
Distribution.h
Harmonics.h
LaserProfile.h
MapGenerator.h
matrix_vector_operation.h
RealDiracMatrix.h
SigmaGenerator.h
)
......
src/Distribution/ClosedOrbitFinder.h
View file @ a21ca6e8
......@@ -40,6 +40,7 @@
#include "Utilities/Options.h"
#include "Utilities/Options.h"
#include "Utilities/OpalException.h"
#include "Physics/Physics.h"
#include "AbstractObjects/OpalData.h"
......
src/Distribution/Distribution.cpp
View file @ a21ca6e8
......@@ -7,7 +7,7 @@
// OPAL is licensed under GNU GPL version 3.
#include "Distribution/Distribution.h"
#include "Distribution/SigmaGenerator.h"
#include "Distribution/ClosedOrbitFinder.h"
#include "AbsBeamline/SpecificElementVisitor.h"
#include <cmath>
......@@ -211,16 +211,6 @@ Distribution::~Distribution() {
delete laserProfile_m;
}
/*
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;
}
}
*/
/**
* Calculate the local number of particles evenly and adjust node 0
......@@ -1269,18 +1259,18 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
bool writeMap = true;
typedef SigmaGenerator<double, unsigned int> sGenerator_t;
sGenerator_t *siggen = new sGenerator_t(I_m,
Attributes::getReal(itsAttr[Attrib::Distribution::EX])*1E6,
Attributes::getReal(itsAttr[Attrib::Distribution::EY])*1E6,
Attributes::getReal(itsAttr[Attrib::Distribution::ET])*1E6,
E_m*1E-6,
massIneV*1E-6,
CyclotronElement,
Nint,
Nsectors,
Attributes::getReal(itsAttr[Attrib::Distribution::ORDERMAPS]),
writeMap);
std::unique_ptr<SigmaGenerator> siggen = std::unique_ptr<SigmaGenerator>(
new SigmaGenerator(I_m,
Attributes::getReal(itsAttr[Attrib::Distribution::EX])*1E6,
Attributes::getReal(itsAttr[Attrib::Distribution::EY])*1E6,
Attributes::getReal(itsAttr[Attrib::Distribution::ET])*1E6,
E_m*1E-6,
massIneV*1E-6,
CyclotronElement,
Nint,
Nsectors,
Attributes::getReal(itsAttr[Attrib::Distribution::ORDERMAPS]),
writeMap));
if (siggen->match(accuracy,
Attributes::getReal(itsAttr[Attrib::Distribution::MAXSTEPSSI]),
......@@ -1288,7 +1278,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
CyclotronElement,
denergy,
rguess,
false, full)) {
full)) {
std::array<double,3> Emit = siggen->getEmittances();
......@@ -1302,62 +1292,18 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
for (unsigned int i = 0; i < siggen->getSigma().size1(); ++ i) {
*gmsg << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,0);
for (unsigned int j = 1; j < siggen->getSigma().size2(); ++ j) {
if (siggen->getSigma()(i,j) < 10e-12){
if (std::abs(siggen->getSigma()(i,j)) < 1.0e-15) {
*gmsg << " & " << std::setprecision(4) << std::setw(11) << 0.0;
}
else{
*gmsg << " & " << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,j);
*gmsg << " & " << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,j);
}
}
*gmsg << " \\\\" << endl;
}
/*
Now setup the distribution generator
Units of the Sigma Matrix:
spatial: mm
moment: rad
*/
double gamma = E_m / massIneV + 1.0;
double beta = std::sqrt(1.0 - 1.0 / (gamma * gamma));
auto sigma = siggen->getSigma();
// change units from mm to m
for (unsigned int i = 0; i < 6; ++ i)
for (unsigned int j = 0; j < 6; ++ j) sigma(i, j) *= 1e-6;
for (unsigned int i = 0; i < 3; ++ i) {
if ( sigma(2 * i, 2 * i) < 0 || sigma(2 * i + 1, 2 * i + 1) < 0 )
throw OpalException("Distribution::CreateMatchedGaussDistribution()",
"Negative value on the diagonal of the sigma matrix.");
}
sigmaR_m[0] = std::sqrt(sigma(0, 0));
sigmaP_m[0] = std::sqrt(sigma(1, 1))*beta*gamma;
sigmaR_m[2] = std::sqrt(sigma(2, 2));
sigmaP_m[2] = std::sqrt(sigma(3, 3))*beta*gamma;
sigmaR_m[1] = std::sqrt(sigma(4, 4));
//p_l^2 = [(delta+1)*beta*gamma]^2 - px^2 - pz^2
double pl2 = (std::sqrt(sigma(5,5)) + 1)*(std::sqrt(sigma(5,5)) + 1)*beta*gamma*beta*gamma -
sigmaP_m[0]*sigmaP_m[0] - sigmaP_m[2]*sigmaP_m[2];
double pl = std::sqrt(pl2);
sigmaP_m[1] = gamma*(pl - beta*gamma);
correlationMatrix_m(1, 0) = sigma(0, 1) / (std::sqrt(sigma(0, 0) * sigma(1, 1)));
correlationMatrix_m(3, 2) = sigma(2, 3) / (std::sqrt(sigma(2, 2) * sigma(3, 3)));
correlationMatrix_m(5, 4) = sigma(4, 5) / (std::sqrt(sigma(4, 4) * sigma(5, 5)));
correlationMatrix_m(4, 0) = sigma(0, 4) / (std::sqrt(sigma(0, 0) * sigma(4, 4)));
correlationMatrix_m(4, 1) = sigma(1, 4) / (std::sqrt(sigma(1, 1) * sigma(4, 4)));
correlationMatrix_m(5, 0) = sigma(0, 5) / (std::sqrt(sigma(0, 0) * sigma(5, 5)));
correlationMatrix_m(5, 1) = sigma(1, 5) / (std::sqrt(sigma(1, 1) * sigma(5, 5)));
createDistributionGauss(numberOfParticles, massIneV);
generateMatchedGauss(siggen->getSigma(), numberOfParticles, massIneV);
// update injection radius and radial momentum
CyclotronElement->setRinit(siggen->getInjectionRadius() * 1.0e3);
......@@ -1366,13 +1312,9 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
else {
*gmsg << "* Not converged for " << E_m*1E-6 << " MeV" << endl;
delete siggen;
throw OpalException("Distribution::CreateMatchedGaussDistribution",
"didn't find any matched distribution.");
}
delete siggen;
}
void Distribution::createDistributionGauss(size_t numberOfParticles, double massIneV) {
......@@ -2412,6 +2354,135 @@ void Distribution::generateGaussZ(size_t numberOfParticles) {
gsl_matrix_free(corMat);
}
void Distribution::generateMatchedGauss(const SigmaGenerator::matrix_t& sigma,
size_t numberOfParticles, double massIneV)
{
/* This particle distribution generation is based on a
* symplectic method described in
* https://arxiv.org/abs/1205.3601
*/
/* Units of the Sigma Matrix:
* spatial: m
* moment: rad
*
* Attention: The vertical and longitudinal direction must be interchanged!
*/
for (unsigned int i = 0; i < 3; ++ i) {
if ( sigma(2 * i, 2 * i) < 0 || sigma(2 * i + 1, 2 * i + 1) < 0 )
throw OpalException("Distribution::generateMatchedGauss()",
"Negative value on the diagonal of the sigma matrix.");
}
double bgam = Util::getBetaGamma(E_m, massIneV);
/*
* only used for printing
*/
// horizontal
sigmaR_m[0] = std::sqrt(sigma(0, 0));
sigmaP_m[0] = std::sqrt(sigma(1, 1)) * bgam;
// longitudinal
sigmaR_m[1] = std::sqrt(sigma(4, 4));
sigmaP_m[1] = std::sqrt(sigma(5, 5)) * bgam;
// vertical
sigmaR_m[2] = std::sqrt(sigma(2, 2));
sigmaP_m[2] = std::sqrt(sigma(3, 3)) * bgam;
correlationMatrix_m(1, 0) = sigma(0, 1) / (std::sqrt(sigma(0, 0) * sigma(1, 1)));
correlationMatrix_m(3, 2) = sigma(2, 3) / (std::sqrt(sigma(2, 2) * sigma(3, 3)));
correlationMatrix_m(5, 4) = sigma(4, 5) / (std::sqrt(sigma(4, 4) * sigma(5, 5)));
correlationMatrix_m(4, 0) = sigma(0, 4) / (std::sqrt(sigma(0, 0) * sigma(4, 4)));
correlationMatrix_m(4, 1) = sigma(1, 4) / (std::sqrt(sigma(1, 1) * sigma(4, 4)));
correlationMatrix_m(5, 0) = sigma(0, 5) / (std::sqrt(sigma(0, 0) * sigma(5, 5)));
correlationMatrix_m(5, 1) = sigma(1, 5) / (std::sqrt(sigma(1, 1) * sigma(5, 5)));
inputMoUnits_m = InputMomentumUnitsT::NONE;
/*
* decouple horizontal and longitudinal direction
*/
// extract horizontal and longitudinal directions
RealDiracMatrix::matrix_t A(4, 4);
A(0, 0) = sigma(0, 0);
A(1, 1) = sigma(1, 1);
A(2, 2) = sigma(4, 4);
A(3, 3) = sigma(5, 5);
A(0, 1) = sigma(0, 1);
A(0, 2) = sigma(0, 4);
A(0, 3) = sigma(0, 5);
A(1, 0) = sigma(1, 0);
A(2, 0) = sigma(4, 0);
A(3, 0) = sigma(5, 0);
A(1, 2) = sigma(1, 4);
A(2, 1) = sigma(4, 1);
A(1, 3) = sigma(1, 5);
A(3, 1) = sigma(5, 1);
A(2, 3) = sigma(4, 5);
A(3, 2) = sigma(5, 4);
RealDiracMatrix rdm;
RealDiracMatrix::sparse_matrix_t R1 = rdm.diagonalize(A);
RealDiracMatrix::vector_t variances(8);
for (int i = 0; i < 4; ++i) {
variances(i) = std::sqrt(A(i, i));
}
/*
* decouple vertical direction
*/
A *= 0.0;
A(0, 0) = 1;
A(1, 1) = 1;
A(2, 2) = sigma(2, 2);
A(3, 3) = sigma(3, 3);
A(2, 3) = sigma(2, 3);
A(3, 2) = sigma(3, 2);
RealDiracMatrix::sparse_matrix_t R2 = rdm.diagonalize(A);
for (int i = 0; i < 4; ++i) {
variances(4 + i) = std::sqrt(A(i, i));
}
xDist_m.resize(numberOfParticles);
pxDist_m.resize(numberOfParticles);
yDist_m.resize(numberOfParticles);
pyDist_m.resize(numberOfParticles);
tOrZDist_m.resize(numberOfParticles);
pzDist_m.resize(numberOfParticles);
RealDiracMatrix::vector_t p1(4), p2(4);
for (size_t i = 0; i < numberOfParticles; i++) {
for (int j = 0; j < 4; j++) {
p1(j) = gsl_ran_gaussian(randGen_m, 1.0) * variances(j);
p2(j) = gsl_ran_gaussian(randGen_m, 1.0) * variances(4 + j);
}
p1 = boost::numeric::ublas::prod(R1, p1);
p2 = boost::numeric::ublas::prod(R2, p2);
xDist_m[i] = p1(0);
pxDist_m[i] = p1(1) * bgam;
yDist_m[i] = p1(2);
pyDist_m[i] = p1(3) * bgam;
tOrZDist_m[i] = p2(2);
pzDist_m[i] = p2(3) * bgam;
}
}
void Distribution::generateLongFlattopT(size_t numberOfParticles) {
double flattopTime = tPulseLengthFWHM_m
......@@ -2971,7 +3042,7 @@ void Distribution::printDistMatchedGauss(Inform &os) const {
os << "* SIGMAPX = " << sigmaP_m[0] << " [Beta Gamma]" << endl;
os << "* SIGMAPY = " << sigmaP_m[1] << " [Beta Gamma]" << endl;
os << "* SIGMAPZ = " << sigmaP_m[2] << " [Beta Gamma]" << endl;
os << "* AVRGPZ = " << avrgpz_m << " [Beta Gamma]" << endl;
// os << "* AVRGPZ = " << avrgpz_m << " [Beta Gamma]" << endl;
os << "* CORRX = " << correlationMatrix_m(1, 0) << endl;
os << "* CORRY = " << correlationMatrix_m(3, 2) << endl;
......@@ -2980,12 +3051,12 @@ void Distribution::printDistMatchedGauss(Inform &os) const {
os << "* R62 = " << correlationMatrix_m(5, 1) << endl;
os << "* R51 = " << correlationMatrix_m(4, 0) << endl;
os << "* R52 = " << correlationMatrix_m(4, 1) << endl;
os << "* CUTOFFX = " << cutoffR_m[0] << " [units of SIGMAX]" << endl;
os << "* CUTOFFY = " << cutoffR_m[1] << " [units of SIGMAY]" << endl;
os << "* CUTOFFLONG = " << cutoffR_m[2] << " [units of SIGMAZ]" << endl;
os << "* CUTOFFPX = " << cutoffP_m[0] << " [units of SIGMAPX]" << endl;
os << "* CUTOFFPY = " << cutoffP_m[1] << " [units of SIGMAPY]" << endl;
os << "* CUTOFFPZ = " << cutoffP_m[2] << " [units of SIGMAPY]" << endl;
// os << "* CUTOFFX = " << cutoffR_m[0] << " [units of SIGMAX]" << endl;
// os << "* CUTOFFY = " << cutoffR_m[1] << " [units of SIGMAY]" << endl;
// os << "* CUTOFFLONG = " << cutoffR_m[2] << " [units of SIGMAZ]" << endl;
// os << "* CUTOFFPX = " << cutoffP_m[0] << " [units of SIGMAPX]" << endl;
// os << "* CUTOFFPY = " << cutoffP_m[1] << " [units of SIGMAPY]" << endl;
// os << "* CUTOFFPZ = " << cutoffP_m[2] << " [units of SIGMAPY]" << endl;
}
void Distribution::printDistGauss(Inform &os) const {
......
src/Distribution/Distribution.h
View file @ a21ca6e8
......@@ -19,6 +19,8 @@
#include "AppTypes/SymTenzor.h"
#include "Attributes/Attributes.h"
#include "Distribution/SigmaGenerator.h"
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_qrng.h>
#include <gsl/gsl_rng.h>
......@@ -279,7 +281,6 @@ private:
void eraseTOrZDist();
void eraseBGzDist();
// void printSigma(SigmaGenerator<double,unsigned int>::matrix_type& M, Inform& out);
void addDistributions();
void applyEmissionModel(double lowEnergyLimit, double &px, double &py, double &pz, std::vector<double> &additionalRNs);
void applyEmissModelAstra(double &px, double &py, double &pz, std::vector<double> &additionalRNs);
......@@ -336,6 +337,8 @@ private:
void generateFlattopT(size_t numberOfParticles);
void generateFlattopZ(size_t numberOfParticles);
void generateGaussZ(size_t numberOfParticles);
void generateMatchedGauss(const SigmaGenerator::matrix_t&,
size_t numberOfParticles, double massIneV);
void generateLongFlattopT(size_t numberOfParticles);
void generateTransverseGauss(size_t numberOfParticles);
void initializeBeam(PartBunchBase<double, 3> *beam);
......
src/Distribution/Harmonics.h deleted 100644 → 0
View file @ cccfc5ee
//
// Class Harmonics
// This class computes the cyclotron map based on harmonics.
// All functions are copied and translated to C++ from the original program inj2_ana.c of Dr. C. Baumgarten.
//
// Copyright (c) 2014 - 2015, Christian Baumgarten, Paul Scherrer Institut, Villigen PSI, Switzerland
// Matthias Frey, ETH Zürich
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL 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 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#ifndef HARMONICS_H
#define HARMONICS_H
#include <cmath>
#include <fstream>
#include <string>
#include <utility>
#include <vector>
#include "Physics/Physics.h"
#include <boost/numeric/ublas/matrix.hpp>
#include "matrix_vector_operation.h"
template<typename Value_type, typename Size_type>
class Harmonics
{
public:
/// Type of variable
typedef Value_type value_type;
/// Type of size
typedef Size_type size_type;
/// Type for specifying matrices
typedef boost::numeric::ublas::matrix<value_type> matrix_type;
/// Defines starting point of computation
enum START { SECTOR_ENTRANCE, SECTOR_CENTER, SECTOR_EXIT, VALLEY_CENTER };
/// Constructor
/*!
* @param wo is the nominal orbital frequency in Hz
* @param Emin is the minimal energy in MeV
* @param Emax is the maximal energy in MeV
* @param nth is the number of angle steps
* @param nr is the number of radial steps
* @param nSector is the number of sectors
* @param E is the kinetic energy [MeV]
* @param E0 is the potential energy [MeV]
*/
Harmonics(value_type, value_type, value_type, size_type, size_type, size_type, value_type, value_type);
/// Compute all maps of the cyclotron using harmonics
/*!
* @param filename1 name of file 1
* @param filename2 name of file 2
* @param startmod (in center of valley, start of sector, ...)
*/
std::vector<matrix_type> computeMap(const std::string, const std::string, const int);
/// Returns the radial and vertical tunes
std::pair<value_type, value_type> getTunes();
/// Returns the radius
value_type getRadius();
/// Returns step size
value_type getPathLength();
private:
/// Stores the nominal orbital frequency in Hz
value_type wo_m;
/// Stores the minimum energy in MeV
value_type Emin_m;
/// Stores the maximum energy in MeV
value_type Emax_m;
/// Stores the number of angle splits
size_type nth_m;
/// Stores the number of radial splits
size_type nr_m;
/// Stores the number of sectors
size_type nSector_m;
/// Stores the tunes (radial, vertical)
std::pair<value_type, value_type> tunes_m;
/// Stores the radius
value_type radius_m;
/// Stores the path length
value_type ds_m;
/// Stores the energy for which we perform the computation
value_type E_m;
/// Stores the potential energy [MeV]
value_type E0_m;
/// Compute some matrix (ask Dr. C. Baumgarten)
matrix_type __Mix6(value_type, value_type, value_type);
/// Compute some matrix (ask Dr. C. Baumgarten)
matrix_type __Md6(value_type, value_type);
/// Compute some matrix (ask Dr. C. Baumgarten)
matrix_type __Mb6(value_type, value_type, value_type);
/// Compute some matrix (ask Dr. C. Baumgarten)
matrix_type __Mb6k(value_type, value_type, value_type, value_type);
};
// -----------------------------------------------------------------------------------------------------------------------
// PUBLIC MEMBER FUNCTIONS
// -----------------------------------------------------------------------------------------------------------------------
template<typename Value_type, typename Size_type>
Harmonics<Value_type, Size_type>::Harmonics(value_type wo, value_type Emin, value_type Emax,
size_type nth, size_type nr, size_type nSector, value_type E, value_type E0)
: wo_m(wo), Emin_m(Emin), Emax_m(Emax), nth_m(nth), nr_m(nr), nSector_m(nSector), E_m(E), E0_m(E0)
{}
template<typename Value_type, typename Size_type>
std::vector<typename Harmonics<Value_type, Size_type>::matrix_type> Harmonics<Value_type, Size_type>::computeMap(const std::string filename1, const std::string filename2, const int startmod) {
// i --> different energies
// j --> different angles
// ---------------
// unsigned int nth_m = 360; //270;
unsigned int N = 4;
value_type two_pi = 2.0 * M_PI;
value_type tpiN = two_pi / value_type(N);
value_type piN = M_PI / value_type(N);
// unsigned int nr_m = 180;
unsigned int cnt = 0;
bool set = false;
value_type gap = 0.05;
value_type K1 = 0.45;
std::vector<matrix_type> Mcyc(nth_m);
value_type beta_m;
std::vector<value_type> gamma(nr_m), E(nr_m), PC(nr_m), nur(nr_m), nuz(nr_m);
std::vector<value_type> R(nr_m), r(nr_m);
std::vector<value_type> Bmag(nr_m), k(nr_m);
std::vector<value_type> alpha(nr_m), beta(nr_m);
//----------------------------------------------------------------------
// read files
std::ifstream infile2(filename2);
std::ifstream infile1(filename1);
if (!infile1.is_open() || !infile2.is_open()) {
std::cerr << "Couldn't open files!" << std::endl;
std::exit(0);
}
unsigned int n = 0;
value_type rN1, cN1, sN1, aN1, phN1, rN2, cN2, sN2, aN2, phN2;
while (infile1 >> rN1 >> cN1 >> sN1 >> aN1 >> phN1 &&
infile2 >> rN2 >> cN2 >> sN2 >> aN2 >> phN2) {
R[n] = rN1 * 0.01;
alpha[n] = 2.0 * std::acos(aN2 / aN1) / value_type(N);
beta[n] = std::atan(sN1 / cN1) / value_type(N);
value_type f4 = aN1 * M_PI * 0.5 / std::sin(0.5 * alpha[n] * value_type(N));
value_type f8 = aN2 * M_PI / std::sin(alpha[n] * value_type(N));
Bmag[n] = 0.5 * (f4 + f8);
n++;
}
infile1.close();
infile2.close();
//----------------------------------------------------------------------
value_type s = std::sin(piN);
value_type c = std::cos(piN);
std::vector<value_type> dalp(nr_m), dbet(nr_m), len2(nr_m), len(nr_m), t1eff(nr_m), t2eff(nr_m);
std::vector<value_type> t1(nr_m), t2(nr_m), t3(nr_m);
dalp[0] = (alpha[1] - alpha[0]) / (R[1] - R[0]);
dalp[n-1] = (alpha[n-1] - alpha[n-2]) / (R[n-1] - R[n-2]);
dbet[0] = (beta[1] - alpha[0]) / (R[1] - R[0]);
dbet[n-1] = (beta[n-1] - beta[n-2]) / (R[n-1] - R[n-2]);
for (unsigned int i = 1; i < n - 1; ++i) {
dalp[i] = R[i] * (alpha[i+1] - alpha[i-1]) / (R[i+1] - R[i-1]);
dbet[i] = R[i] * (beta[i+1] - beta[i-1]) / (R[i+1] - R[i-1]);
}
for (unsigned int i = 0; i < n; ++i) {
value_type cot = 1.0 / std::tan(0.5 * alpha[i]);
value_type eps1 = std::atan(dbet[i] - 0.5 * dalp[i]);
value_type eps2 = std::atan(dbet[i] + 0.5 * dalp[i]);
value_type phi = piN - 0.5 * alpha[i];
// bending radius
r[i] = R[i] * sin(0.5 * alpha[i]) / s;
len2[i] = R[i] * std::sin(phi);
len[i] = two_pi * r[i] + 2.0 * len2[i] * value_type(N);
value_type g1 = phi + eps1;
value_type g2 = - phi + eps2;
t1[i] = std::tan(g1);
t2[i] = std::tan(g2);
t3[i] = std::tan(phi);
value_type fac = gap / r[i] * K1;
value_type psi = fac * (1.0 + std::sin(g1) * std::sin(g1)) / std::cos(g1);
t1eff[i] = std::tan(g1 - psi);
psi = fac * (1.0 + std::sin(g2) *