From 8651496f089afa3ce52db5a60fcfd124e703f5f8 Mon Sep 17 00:00:00 2001 From: 2xB <31772910+2xB@users.noreply.github.com> Date: Wed, 5 Jun 2024 15:44:21 +0200 Subject: [PATCH 01/15] Docker: Fedora 38 => 40 --- Dockerfile | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/Dockerfile b/Dockerfile index 932c5545e..283dbbcc9 100644 --- a/Dockerfile +++ b/Dockerfile @@ -10,8 +10,7 @@ ARG KASSIOPEIA_GIT_COMMIT="" ARG KASSIOPEIA_CPUS="" # --- runtime-base --- -# NOTE: For Fedora 39, remove marked lines below -FROM fedora:38 as runtime-base +FROM fedora:40 as runtime-base ARG KASSIOPEIA_UID ARG KASSIOPEIA_USER ARG KASSIOPEIA_GID @@ -19,13 +18,6 @@ ARG KASSIOPEIA_GROUP LABEL description="Runtime base container" -# TODO REMOVE FOR FEDORA 39 -RUN dnf update -y \ - && dnf install -y --setopt=install_weak_deps=False dnf-plugins-core \ - && dnf clean all -RUN dnf copr enable thofmann/log4xx-1.x -y -# END TODO - COPY Docker/packages.runtime packages RUN dnf update -y \ && dnf install -y --setopt=install_weak_deps=False $(cat packages) \ From 85fbde81a906437781ad84b184be60e6cd59f007 Mon Sep 17 00:00:00 2001 From: 2xB <2xB@users.noreply.github.com> Date: Wed, 8 May 2024 11:21:13 +0200 Subject: [PATCH 02/15] Fix warnings about adding controls/interactions multiple times Before this commit, some controls/interactions showed warnings about being added twice even if just the same parent type was added twice. This is an alternative approach to https://github.com/KATRIN-Experiment/Kassiopeia/pull/76 , just showing the messages if relevant instead of changing their log level and content. --- Kassiopeia/Utility/Include/KSList.h | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/Kassiopeia/Utility/Include/KSList.h b/Kassiopeia/Utility/Include/KSList.h index 3ee739900..894b3e55d 100644 --- a/Kassiopeia/Utility/Include/KSList.h +++ b/Kassiopeia/Utility/Include/KSList.h @@ -1,6 +1,8 @@ #ifndef KSLIST_H_ #define KSLIST_H_ +#include + namespace Kassiopeia { @@ -148,16 +150,16 @@ template int KSList::FindElement(XType* anElement) } return tIndex; } -template int KSList::FindElementByType(XType*) +template int KSList::FindElementByType(XType* anElement) { - int tIndex = -1; - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - if (dynamic_cast(fElements[fCurrentElement]) != nullptr) { - tIndex = fCurrentElement; - break; + const std::type_info& elementType = typeid(*anElement); + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + const std::type_info& compareType = typeid(*fElements[currentElement]); + if (elementType == compareType) { + return currentElement; } } - return tIndex; + return -1; } template int KSList::RemoveElement(XType* anElement) { From 9fca81101653b265852f7332ec6af906afcedff2 Mon Sep 17 00:00:00 2001 From: 2xB <2xB@users.noreply.github.com> Date: Wed, 8 May 2024 11:21:53 +0200 Subject: [PATCH 03/15] KSList: Avoid unneeded mutable variable --- Kassiopeia/Utility/Include/KSList.h | 142 ++++++++++++++-------------- 1 file changed, 69 insertions(+), 73 deletions(-) diff --git a/Kassiopeia/Utility/Include/KSList.h b/Kassiopeia/Utility/Include/KSList.h index 894b3e55d..1bedce959 100644 --- a/Kassiopeia/Utility/Include/KSList.h +++ b/Kassiopeia/Utility/Include/KSList.h @@ -72,25 +72,24 @@ template class KSList private: XType** fElements; - mutable int fCurrentElement; int fEndElement; const int fLastElement; }; template -KSList::KSList(const int aMaxSize) : fCurrentElement(0), fEndElement(0), fLastElement(aMaxSize) +KSList::KSList(const int aMaxSize) : fEndElement(0), fLastElement(aMaxSize) { fElements = new XType*[fLastElement]; - for (fCurrentElement = 0; fCurrentElement < fLastElement; fCurrentElement++) { - fElements[fCurrentElement] = nullptr; + for (int currentElement = 0; currentElement < fLastElement; currentElement++) { + fElements[currentElement] = nullptr; } } template -KSList::KSList(const KSList& aCopy) : fCurrentElement(0), fEndElement(0), fLastElement(aCopy.fLastElement) +KSList::KSList(const KSList& aCopy) : fEndElement(0), fLastElement(aCopy.fLastElement) { fElements = new XType*[fLastElement]; - for (fCurrentElement = 0; fCurrentElement < fLastElement; fCurrentElement++) { - fElements[fCurrentElement] = aCopy.fElements[fCurrentElement]; + for (int currentElement = 0; currentElement < fLastElement; currentElement++) { + fElements[currentElement] = aCopy.fElements[currentElement]; } } template KSList::~KSList() @@ -109,9 +108,9 @@ template bool KSList::Full() const template void KSList::ClearContent() { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - if (fElements[fCurrentElement] != nullptr) { - fElements[fCurrentElement] = nullptr; + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + if (fElements[currentElement] != nullptr) { + fElements[currentElement] = nullptr; } } fEndElement = 0; @@ -119,10 +118,10 @@ template void KSList::ClearContent() } template void KSList::DeleteContent() { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - if (fElements[fCurrentElement] != nullptr) { - delete fElements[fCurrentElement]; - fElements[fCurrentElement] = nullptr; + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + if (fElements[currentElement] != nullptr) { + delete fElements[currentElement]; + fElements[currentElement] = nullptr; } } fEndElement = 0; @@ -134,21 +133,19 @@ template int KSList::AddElement(XType* anElement) if (fEndElement == fLastElement) { return -1; } - fCurrentElement = fEndElement; - fElements[fCurrentElement] = anElement; + const int currentElement = fEndElement; + fElements[currentElement] = anElement; fEndElement++; - return fCurrentElement; + return currentElement; } template int KSList::FindElement(XType* anElement) { - int tIndex = -1; - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - if (fElements[fCurrentElement] == anElement) { - tIndex = fCurrentElement; - break; + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + if (fElements[currentElement] == anElement) { + return currentElement; } } - return tIndex; + return -1; } template int KSList::FindElementByType(XType* anElement) { @@ -164,15 +161,15 @@ template int KSList::FindElementByType(XType* anElement) template int KSList::RemoveElement(XType* anElement) { int tIndex = -1; - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { if (tIndex == -1) { - if (fElements[fCurrentElement] == anElement) { - fElements[fCurrentElement] = nullptr; - tIndex = fCurrentElement; + if (fElements[currentElement] == anElement) { + fElements[currentElement] = nullptr; + tIndex = currentElement; } continue; } - fElements[fCurrentElement - 1] = fElements[fCurrentElement]; + fElements[currentElement - 1] = fElements[currentElement]; } if (tIndex != -1) { fEndElement--; @@ -182,19 +179,18 @@ template int KSList::RemoveElement(XType* anElement) } template XType* KSList::ElementAt(int& anIndex) const { - fCurrentElement = anIndex; - if (fCurrentElement >= fEndElement || fCurrentElement < 0) { + if (anIndex >= fEndElement || anIndex < 0) { return nullptr; } - return fElements[fCurrentElement]; + return fElements[anIndex]; } template template void KSList::ForEach(XReturnType (XClassType::*aMember)()) { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - (fElements[fCurrentElement]->*aMember)(); + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + (fElements[currentElement]->*aMember)(); } return; } @@ -202,8 +198,8 @@ template template void KSList::ForEach(XReturnType (XClassType::*aMember)() const) { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - (fElements[fCurrentElement]->*aMember)(); + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + (fElements[currentElement]->*aMember)(); } return; } @@ -211,8 +207,8 @@ template template void KSList::ForEach(XReturnType (XClassType::*aMember)(XArgumentType), XArgumentType anArgument) { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - (fElements[fCurrentElement]->*aMember)(anArgument); + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + (fElements[currentElement]->*aMember)(anArgument); } return; } @@ -220,8 +216,8 @@ template template void KSList::ForEach(XReturnType (XClassType::*aMember)(XArgumentType) const, XArgumentType anArgument) { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - (fElements[fCurrentElement]->*aMember)(anArgument); + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + (fElements[currentElement]->*aMember)(anArgument); } return; } @@ -231,8 +227,8 @@ template::ForEach(XReturnType (XClassType::*aMember)(XArgumentType1, XArgumentType2), XArgumentType1 anArgument1, XArgumentType2 anArgument2) { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - (fElements[fCurrentElement]->*aMember)(anArgument1, anArgument2); + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + (fElements[currentElement]->*aMember)(anArgument1, anArgument2); } return; } @@ -241,16 +237,16 @@ template::ForEach(XReturnType (XClassType::*aMember)(XArgumentType1, XArgumentType2) const, XArgumentType1 anArgument1, XArgumentType2 anArgument2) { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - (fElements[fCurrentElement]->*aMember)(anArgument1, anArgument2); + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + (fElements[currentElement]->*aMember)(anArgument1, anArgument2); } return; } template template bool KSList::ForEachUntilTrue(bool (XClassType::*aMember)()) { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - if ((fElements[fCurrentElement]->*aMember)() == true) { + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + if ((fElements[currentElement]->*aMember)() == true) { return true; } } @@ -260,8 +256,8 @@ template template bool KSList::ForEachUntilTrue(bool (XClassType::*aMember)() const) { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - if ((fElements[fCurrentElement]->*aMember)() == true) { + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + if ((fElements[currentElement]->*aMember)() == true) { return true; } } @@ -271,8 +267,8 @@ template template bool KSList::ForEachUntilTrue(bool (XClassType::*aMember)(XArgumentType), XArgumentType anArgument) { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - if ((fElements[fCurrentElement]->*aMember)(anArgument) == true) { + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + if ((fElements[currentElement]->*aMember)(anArgument) == true) { return true; } } @@ -282,8 +278,8 @@ template template bool KSList::ForEachUntilTrue(bool (XClassType::*aMember)(XArgumentType) const, XArgumentType anArgument) { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - if ((fElements[fCurrentElement]->*aMember)(anArgument) == true) { + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + if ((fElements[currentElement]->*aMember)(anArgument) == true) { return true; } } @@ -294,8 +290,8 @@ template bool KSList::ForEachUntilTrue(bool (XClassType::*aMember)(XArgumentType1, XArgumentType2), XArgumentType1 anArgument1, XArgumentType2 anArgument2) { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - if ((fElements[fCurrentElement]->*aMember)(anArgument1, anArgument2) == true) { + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + if ((fElements[currentElement]->*aMember)(anArgument1, anArgument2) == true) { return true; } } @@ -306,8 +302,8 @@ template bool KSList::ForEachUntilTrue(bool (XClassType::*aMember)(XArgumentType1, XArgumentType2) const, XArgumentType1 anArgument1, XArgumentType2 anArgument2) { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - if ((fElements[fCurrentElement]->*aMember)(anArgument1, anArgument2) == true) { + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + if ((fElements[currentElement]->*aMember)(anArgument1, anArgument2) == true) { return true; } } @@ -316,8 +312,8 @@ bool KSList::ForEachUntilTrue(bool (XClassType::*aMember)(XArgumentType1, template template bool KSList::ForEachUntilFalse(bool (XClassType::*aMember)()) { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - if ((fElements[fCurrentElement]->*aMember)() == false) { + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + if ((fElements[currentElement]->*aMember)() == false) { return false; } } @@ -327,8 +323,8 @@ template template bool KSList::ForEachUntilFalse(bool (XClassType::*aMember)() const) { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - if ((fElements[fCurrentElement]->*aMember)() == false) { + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + if ((fElements[currentElement]->*aMember)() == false) { return false; } } @@ -338,8 +334,8 @@ template template bool KSList::ForEachUntilFalse(bool (XClassType::*aMember)(XArgumentType), XArgumentType anArgument) { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - if ((fElements[fCurrentElement]->*aMember)(anArgument) == false) { + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + if ((fElements[currentElement]->*aMember)(anArgument) == false) { return false; } } @@ -349,8 +345,8 @@ template template bool KSList::ForEachUntilFalse(bool (XClassType::*aMember)(XArgumentType) const, XArgumentType anArgument) { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - if ((fElements[fCurrentElement]->*aMember)(anArgument) == false) { + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + if ((fElements[currentElement]->*aMember)(anArgument) == false) { return false; } } @@ -361,8 +357,8 @@ template bool KSList::ForEachUntilFalse(bool (XClassType::*aMember)(XArgumentType1, XArgumentType2), XArgumentType1 anArgument1, XArgumentType2 anArgument2) { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - if ((fElements[fCurrentElement]->*aMember)(anArgument1, anArgument2) == false) { + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + if ((fElements[currentElement]->*aMember)(anArgument1, anArgument2) == false) { return false; } } @@ -373,8 +369,8 @@ template bool KSList::ForEachUntilFalse(bool (XClassType::*aMember)(XArgumentType1, XArgumentType2) const, XArgumentType1 anArgument1, XArgumentType2 anArgument2) { - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - if ((fElements[fCurrentElement]->*aMember)(anArgument1, anArgument2) == false) { + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + if ((fElements[currentElement]->*aMember)(anArgument1, anArgument2) == false) { return false; } } @@ -387,8 +383,8 @@ XReturnType KSList::LargestOfAll(XReturnType (XClassType::*aMember)()) { XReturnType Largest = -1.e300; XReturnType CurrentValue; - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - CurrentValue = (fElements[fCurrentElement]->*aMember)(); + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + CurrentValue = (fElements[currentElement]->*aMember)(); if (CurrentValue > Largest) { Largest = CurrentValue; } @@ -402,8 +398,8 @@ XReturnType KSList::SmallestOfAll(XReturnType (XClassType::*aMember)()) { XReturnType Smallest = 1.e300; XReturnType CurrentValue; - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - CurrentValue = (fElements[fCurrentElement]->*aMember)(); + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + CurrentValue = (fElements[currentElement]->*aMember)(); if (CurrentValue < Smallest) { Smallest = CurrentValue; } @@ -416,8 +412,8 @@ template XReturnType KSList::SumOfAll(XReturnType (XClassType::*aMember)()) { XReturnType Sum = 0.; - for (fCurrentElement = 0; fCurrentElement < fEndElement; fCurrentElement++) { - Sum += (fElements[fCurrentElement]->*aMember)(); + for (int currentElement = 0; currentElement < fEndElement; currentElement++) { + Sum += (fElements[currentElement]->*aMember)(); } return Sum; } From ce14ed9ff41b63b7e185890cbbc0c29468637ff8 Mon Sep 17 00:00:00 2001 From: Richard Salomon Date: Tue, 21 Nov 2023 10:31:18 +0100 Subject: [PATCH 04/15] Fix compiler warnings in MagfieldCoils.cxx MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit When building the class, several compiler warnings were triggered: [build] [566/1263 43% :: 34.774] Building CXX object KEMField/Source/ExternalFields/MagfieldCoils/CMakeFiles/KEMMagfieldCoils.dir/src/MagfieldCoils.cxx.o [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx: In member function ‘void MagfieldCoils::CoilRead()’: [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx:170:58: warning: unused variable ‘v’ [-Wunused-variable] [build] 170 | double cu, Cx, Cy, Cz, alpha, beta, tu, L, Rmin,Rmax, v[3]; [build] | ^ [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx: In member function ‘void MagfieldCoils::Magfield2EllipticCoil(int, double, double, double&, double&)’: [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx:661:16: warning: unused variable ‘st’ [-Wunused-variable] [build] 661 | double sign, st, delRr, Rlow[2], Rhigh[2]; [build] | ^~ [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx: In member function ‘void MagfieldCoils::Magsource2RemoteCoil(int, double, double)’: [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx:873:44: warning: unused variable ‘st’ [-Wunused-variable] [build] 873 | double L, sigma, Zmin, Zmax, Rmin, Rmax, st; [build] | ^~ [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx: In member function ‘void MagfieldCoils::RemoteSourcepointGroup(int)’: [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx:1046:2: warning: this ‘if’ clause does not guard... [-Wmisleading-indentation] [build] 1046 | if(zAzmax) zmax=zB; [build] | ^~ [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx:1046:25: note: ...this statement, but the latter is misleadingly indented as if it were guarded by the ‘if’ [build] 1046 | if(zAzmax) zmax=zB; [build] | ^~ [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx: In member function ‘void MagfieldCoils::MagsourceMagchargeCoils()’: [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx:1079:22: warning: variable ‘L’ set but not used [-Wunused-but-set-variable] [build] 1079 | double Rmin, Rmax, L, rorem, sigma; [build] | ^ [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx:1079:25: warning: unused variable ‘rorem’ [-Wunused-variable] [build] 1079 | double Rmin, Rmax, L, rorem, sigma; [build] | ^~~~~ [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx: In member function ‘void MagfieldCoils::Magsource2CentralCoil(int, double, double)’: [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx:1190:44: warning: unused variable ‘st’ [-Wunused-variable] [build] 1190 | double L, sigma, Zmin, Zmax, Rmin, Rmax, st; [build] | ^~ [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx: In member function ‘void MagfieldCoils::MagsourceCentralCoils()’: [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx:1262:10: warning: unused variable ‘z0’ [-Wunused-variable] [build] 1262 | double z0, Rmin, Rmax, L, rorem, cu, tu; [build] | ^~ [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx:1262:14: warning: unused variable ‘Rmin’ [-Wunused-variable] [build] 1262 | double z0, Rmin, Rmax, L, rorem, cu, tu; [build] | ^~~~ [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx:1262:20: warning: unused variable ‘Rmax’ [-Wunused-variable] [build] 1262 | double z0, Rmin, Rmax, L, rorem, cu, tu; [build] | ^~~~ [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx:1262:26: warning: unused variable ‘L’ [-Wunused-variable] [build] 1262 | double z0, Rmin, Rmax, L, rorem, cu, tu; [build] | ^ [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx:1262:29: warning: unused variable ‘rorem’ [-Wunused-variable] [build] 1262 | double z0, Rmin, Rmax, L, rorem, cu, tu; [build] | ^~~~~ [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx:1262:36: warning: unused variable ‘cu’ [-Wunused-variable] [build] 1262 | double z0, Rmin, Rmax, L, rorem, cu, tu; [build] | ^~ [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx:1262:40: warning: unused variable ‘tu’ [-Wunused-variable] [build] 1262 | double z0, Rmin, Rmax, L, rorem, cu, tu; [build] | ^~ [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx: In member function ‘void MagfieldCoils::CentralSourcepointsGroup(int)’: [build] ../KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx:1431:16: warning: unused variable ‘rocen’ [-Wunused-variable] [build] 1431 | double z0, rocen; [build] | ^~~~~ --- .../MagfieldCoils/src/MagfieldCoils.cxx | 30 ++++++++++--------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx b/KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx index bc4a6c63e..8e3f62f39 100755 --- a/KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx +++ b/KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx @@ -167,7 +167,7 @@ void MagfieldCoils::CoilRead() coil[i]= new double[14]; } // Reading the coil parameters: - double cu, Cx, Cy, Cz, alpha, beta, tu, L, Rmin,Rmax, v[3]; + double cu, Cx, Cy, Cz, alpha, beta, tu, L, Rmin, Rmax; for(int i=0; i> cu >> Cx >> Cy >> Cz >> alpha >> beta >> tu >> L >> Rmin >> Rmax ; @@ -658,7 +658,7 @@ void MagfieldCoils::Magfield2EllipticCoil(int i, double z, double r, double& Bz, // the first, second and third complete elliptic integrals. double L, Zmin, Zmax, Rmin, Rmax, sigma; double R, Z, delr2, sumr2, delz2, eta, d, K, EK, PIK, S; - double sign, st, delRr, Rlow[2], Rhigh[2]; + double sign, delRr, Rlow[2], Rhigh[2]; const double mu0=4.*M_PI*1.e-7; double x[2][1001], w[2][1001]; // nodes and weights // Coil parameters: @@ -870,7 +870,7 @@ void MagfieldCoils::Magsource2RemoteCoil(int i, double z0, double rorem) { const double mu0=4.*M_PI*1.e-7; double x[1001], w[1001]; // nodes and weights - double L, sigma, Zmin, Zmax, Rmin, Rmax, st; + double L, sigma, Zmin, Zmax, Rmin, Rmax; // Coil parameters: L=coil[i][7]; // coil length Zmin=-L/2.; Zmax=L/2.; // coil endpoints relative to coil center @@ -1043,7 +1043,10 @@ void MagfieldCoils::RemoteSourcepointGroup(int g) { int i=Cin[g][c]; double L=coil[i][7]; double zA=Z[g][c]-L/2.; double zB=Z[g][c]+L/2.; // coil edges - if(zAzmax) zmax=zB; + if (zA < zmin) + zmin = zA; + if (zB > zmax) + zmax = zB; } z0rem=z0remG[g]=(zmin+zmax)/2.; // group center in group Z-coordinate system // Remote convergence radius of group --> roremG[g] : @@ -1076,7 +1079,7 @@ void MagfieldCoils::MagsourceMagchargeCoils() // i: coil index, n source constant index. { // Output to file dirname+objectname+magsource_remote.txt : - double Rmin, Rmax, L, rorem, sigma; + double Rmin, Rmax, sigma; string filename=dirname+objectname+"magsource_remote.txt"; ofstream output; output.precision(16); @@ -1089,7 +1092,6 @@ void MagfieldCoils::MagsourceMagchargeCoils() for(int i=0; i Z0cen; - double z0, rocen; + double z0; double roremC=roremG[g]; double z0min=z0remG[g]-2.*roremC; double z0max=z0remG[g]+2.*roremC; @@ -1524,7 +1525,7 @@ bool MagfieldCoils::Magfield2Remote(bool bcoil, int ig, double z, double r, doub } // We start here the series expansion: Bz=Br=0.; - int nlast; + int nlast = 0; double brem; if(bcoil==true || even==true) // even calculation { @@ -1699,7 +1700,7 @@ bool MagfieldCoils::Hfield(int i, double z, double r, double& Hz, double& Hr, do // We start here the series expansion. // Only the odd-n terms are needed in the series. double Hzplus,Hrplus,Hzplusold,Hrplusold,Heps,Hdelta; - int nlast; + int nlast=0; for(int n=1; n<=nmax-1; n+=2) { nlast=n; @@ -1805,10 +1806,10 @@ bool MagfieldCoils::Magfield2Central(bool bcoil, int ig, int j, double z, double Bz=BcenG[ig][j][0]+BcenG[ig][j][1]*rc*u; Br=-sr*BcenG[ig][j][1]/2.*rc; } + int nlast=0; if(rc<1.e-10) goto label; // // We start here the central series expansion: - int nlast; double bcen; for(int n=2; n<=nmax-1; n++) { @@ -2339,8 +2340,9 @@ double MagfieldCoils::RJ_Carlson(double x,double y,double z,double p) // (see: Press et al., Numerical Recipes, Sec. 6.11). const double ERRTOL=0.0015,TINY=1.e-20,BIG=1.e12,C1=3./14.,C2=1./3., C3=3./22.,C4=3./26.,C5=0.75*C3,C6=1.5*C4,C7=0.5*C2,C8=2.*C3; - double a,alamb,alpha,ans,ave,b,beta,delp,delx,dely,delz,ea,eb,ec,ed,ee, - fac,pt,rcx,rho,sum,sqrtx,sqrty,sqrtz,tau,xt,yt,zt; + double alamb,alpha,ans,ave,beta,delp,delx,dely,delz,ea,eb,ec,ed,ee, + fac,pt,rho,sum,sqrtx,sqrty,sqrtz,tau,xt,yt,zt; + double rcx = 0.; double b = 0.; double a = 0.; if(FMIN3(x,y,z)<0. || FMIN(FMIN(x+y,x+z),FMIN(y+z,fabs(p)))BIG) { From 918c5bad804442e51092575a104d39192581780b Mon Sep 17 00:00:00 2001 From: Richard Salomon Date: Tue, 21 Nov 2023 10:40:27 +0100 Subject: [PATCH 05/15] clang format MagfieldCoils.cxx to make it more readable --- .../MagfieldCoils/include/MagfieldCoils.h | 229 +- .../MagfieldCoils/src/MagfieldCoils.cxx | 4331 +++++++++-------- 2 files changed, 2331 insertions(+), 2229 deletions(-) diff --git a/KEMField/Source/ExternalFields/MagfieldCoils/include/MagfieldCoils.h b/KEMField/Source/ExternalFields/MagfieldCoils/include/MagfieldCoils.h index 956063319..70b95eac6 100755 --- a/KEMField/Source/ExternalFields/MagfieldCoils/include/MagfieldCoils.h +++ b/KEMField/Source/ExternalFields/MagfieldCoils/include/MagfieldCoils.h @@ -3,14 +3,14 @@ #ifndef MagfieldCoils_h #define MagfieldCoils_h -#include -#include +#include #include -#include #include #include +#include +#include +#include #include -#include using namespace std; @@ -20,108 +20,113 @@ using namespace std; class MagfieldCoils { public: - MagfieldCoils(string inputdirname, string inputobjectname, string inputcoilfilename, - int inputNelliptic, int inputnmax, double inputepstol); // constructor - MagfieldCoils(string inputdirname, string inputobjectname); - ~MagfieldCoils(); // deconstructor - bool Magfield(const double *P, double *B); - void MagfieldElliptic(const double *P, double *B); - void SetTrackingStart(); + MagfieldCoils(string inputdirname, string inputobjectname, string inputcoilfilename, int inputNelliptic, + int inputnmax, double inputepstol); // constructor + MagfieldCoils(string inputdirname, string inputobjectname); + ~MagfieldCoils(); // deconstructor + bool Magfield(const double* P, double* B); + void MagfieldElliptic(const double* P, double* B); + void SetTrackingStart(); + private: -// MEMBER FUNCTIONS: -// General: - void CoilRead(); - void CoilGroupWrite(); - void CoilGroupRead(); - void MagsourceRead(); - void DynamicMemoryAllocations(); -// Elliptic: - void Magfield2EllipticCoil(int i, double z, double r, double& Bz,double& Br); - void MagfieldEllipticCoil(int i, const double *P, double *B); -// Source Remote: - double Funrorem(int i, double z0); - void Magsource2RemoteCoil(int i, double z0, double rorem); - void MagsourceRemoteCoils(); - void MagsourceRemoteGroups(); - void RemoteSourcepointGroup(int g); -// Source Magcharge: + // MEMBER FUNCTIONS: + // General: + void CoilRead(); + void CoilGroupWrite(); + void CoilGroupRead(); + void MagsourceRead(); + void DynamicMemoryAllocations(); + // Elliptic: + void Magfield2EllipticCoil(int i, double z, double r, double& Bz, double& Br); + void MagfieldEllipticCoil(int i, const double* P, double* B); + // Source Remote: + double Funrorem(int i, double z0); + void Magsource2RemoteCoil(int i, double z0, double rorem); + void MagsourceRemoteCoils(); + void MagsourceRemoteGroups(); + void RemoteSourcepointGroup(int g); + // Source Magcharge: void MagsourceMagchargeCoils(); -// Source Central: - double Funrocen(int i, double z0); - double FunrocenG(int g, double z0); - void Magsource2CentralCoil(int i, double z0, double rocen); - void MagsourceCentralCoils(); - void CentralSourcepointsCoil(int i); - void MagsourceCentralGroups(); - void CentralSourcepointsGroup(int g); -// Magfield Remote: - bool Magfield2Remote(bool coil, int ig, double z, double r, double& Bz,double& Br,double& rc); -// Magfield Magcharge: + // Source Central: + double Funrocen(int i, double z0); + double FunrocenG(int g, double z0); + void Magsource2CentralCoil(int i, double z0, double rocen); + void MagsourceCentralCoils(); + void CentralSourcepointsCoil(int i); + void MagsourceCentralGroups(); + void CentralSourcepointsGroup(int g); + // Magfield Remote: + bool Magfield2Remote(bool coil, int ig, double z, double r, double& Bz, double& Br, double& rc); + // Magfield Magcharge: bool Magfield2Magcharge(int i, double z, double r, double& Bz, double& Br, double& rc); bool Hfield(int i, double z, double r, double& Hz, double& Hr, double& rc); -// Magfield Central: - bool Magfield2Central(bool coil, int ig, int j, double z, double r, double& Bz,double& Br,double& rc); -// Magfield Coil, Group: - bool MagfieldCoil(int i, const double *P, double *B); - bool MagfieldGroup(int g, const double *P, double *B); + // Magfield Central: + bool Magfield2Central(bool coil, int ig, int j, double z, double r, double& Bz, double& Br, double& rc); + // Magfield Coil, Group: + bool MagfieldCoil(int i, const double* P, double* B); + bool MagfieldGroup(int g, const double* P, double* B); bool SourcePointSearching(bool coil, int ig, double z, double r, int type, int& jbest, double& rcbest); -// Integration: - void GaussLegendreIntegration(int& N, double a, double b, double* x, double* w); -// Carlson elliptic: - double RF_Carlson(double x,double y,double z); - double RD_Carlson(double x,double y,double z); - double RJ_Carlson(double x,double y,double z,double p); - double RC_Carlson(double x,double y); -// Simple: - double pow2(double x); - double FMIN(double a, double b); - double FMAX(double a, double b); - double FMIN3(double a, double b, double c); - double FMAX3(double a, double b, double c); - double Mag(double* vec); - double Mag(double x, double y, double z); -// VARIABLES: - int Ncoil; // number of coils; coil indexing: i=0, ..., Ncoil-1 - double** coil; // coil parameters: coil[i][j], i=0, ..., Ncoil-1; j=0,...,13. - string dirname; // name of directory where the coil and source coefficient data files are stored - string coilfilename; // name of coil input file - string objectname; // this string identifies the Magfield object (used as starting part of the source files) - int Nelliptic; // radial numerical integration parameter for magnetic field calc. with elliptic integrals - int nmax; // number of source coefficients for fixed z: nmax+1 - double epstol; // small tolerance parameter needed for the symmetry group definitions (distance: in meter) - int Ng; // number of coil symmetry groups; group indexing: g=0, ..., Ng-1 - int* G; // coil with global index i is in symmetry group g=G[i], i=0,...,Ncoil-1 - int* Nc; // number of coils in symmetry group g: Nc[g]; local coil indexing in group g: c=0,..., Nc[g]-1 - int** Cin; // global coil indices i in symmetry group g: i=Cin[g][c],...,c=0 to Nc[g]-1 (g=0,...,Ng-1) - double** Line; // symmetry group axis line of symmetry group g (g=0,...,Ng-1): -// line point: Line[g][0], Line[g][1], Line[g][2]; line direction unit vector: Line[g][3], Line[g][4], Line[g][5]; - double** Z; // local z coordinate of coil center in symmetry group g: Z[g][c], c=0,..., Nc[g]-1; Z[g][0]=0 - // (g: group index, c: local coil index) - double *C0, *C1, *c1, *c2, *c3, *c4, *c5, *c6, *c7, *c8, *c9, *c10, *c11, *c12; - - double* Pp, *P; // Legendre polynomial P and first derivative Pp - double* Bzplus, *Brplus; // used in Bz and Br mag. field calc. -// Remote source constant variables: - double* Brem1; // 1-dim. remote source constants Brem1[n], n=2,...,nmax. - double* rorem; // remote convergence radius for coil i: rorem[i] - double** Brem; // remote source constant for coil i: Brem[i][n] (coil index i, source constant index n) - double* z0remG; // remote source point for group G: z0remG[g] - double* roremG; // remote convergence radius for group G: roremG[g] - double** BremG; // remote source constant for symmetry group g: BremG[g][n] (group index g, source constant index n) -// Magnetic charge source constant variables: - double** Vrem; // magnetic charge remote source constant for coil i: Vrem[i][n] (coil index i, source const. index n) -// Central source constant variables: - double* Bcen1; // 1-dim. central source constants Bcen1[n], n=0,...,nmax. - int* Nsp; // number of central source points for coil i: Nsp[i] (i=0,...,Ncoil-1) - double** z0cen; // central source point local z value for coil i: z0cen[i][j] (coil index i, source point index j) - double** rocen; // central convergence radius for coil i: rocen[i][j] (coil index i, source point index j) - double*** Bcen; // central source constant for coil i: Bcen[i][j][n] (coil index i, source point index j, source constant index n) - int* NspG; // number of central source points for group g: Nsp[g] (g=0,...,Ng-1) - double** z0cenG; // central source point local z value: z0cenG[g][j] (group index g, source point index j) - double** rocenG; // central convergence radius: rocenG[g][j] (group index g, source point index j) - double*** BcenG; // central source constant: BcenG[g][j][n] (group index g, source point index j, source const. index n) - int *jlast, *jlastG; // last central source point index for coil and group calculation -double rclimit; + // Integration: + void GaussLegendreIntegration(int& N, double a, double b, double* x, double* w); + // Carlson elliptic: + double RF_Carlson(double x, double y, double z); + double RD_Carlson(double x, double y, double z); + double RJ_Carlson(double x, double y, double z, double p); + double RC_Carlson(double x, double y); + // Simple: + double pow2(double x); + double FMIN(double a, double b); + double FMAX(double a, double b); + double FMIN3(double a, double b, double c); + double FMAX3(double a, double b, double c); + double Mag(double* vec); + double Mag(double x, double y, double z); + // VARIABLES: + int Ncoil; // number of coils; coil indexing: i=0, ..., Ncoil-1 + double** coil; // coil parameters: coil[i][j], i=0, ..., Ncoil-1; j=0,...,13. + string dirname; // name of directory where the coil and source coefficient data files are stored + string coilfilename; // name of coil input file + string objectname; // this string identifies the Magfield object (used as starting part of the source files) + int Nelliptic; // radial numerical integration parameter for magnetic field calc. with elliptic integrals + int nmax; // number of source coefficients for fixed z: nmax+1 + double epstol; // small tolerance parameter needed for the symmetry group definitions (distance: in meter) + int Ng; // number of coil symmetry groups; group indexing: g=0, ..., Ng-1 + int* G; // coil with global index i is in symmetry group g=G[i], i=0,...,Ncoil-1 + int* Nc; // number of coils in symmetry group g: Nc[g]; local coil indexing in group g: c=0,..., Nc[g]-1 + int** Cin; // global coil indices i in symmetry group g: i=Cin[g][c],...,c=0 to Nc[g]-1 (g=0,...,Ng-1) + double** Line; // symmetry group axis line of symmetry group g (g=0,...,Ng-1): + // line point: Line[g][0], Line[g][1], Line[g][2]; line direction unit vector: Line[g][3], Line[g][4], Line[g][5]; + double** Z; // local z coordinate of coil center in symmetry group g: Z[g][c], c=0,..., Nc[g]-1; Z[g][0]=0 + // (g: group index, c: local coil index) + double *C0, *C1, *c1, *c2, *c3, *c4, *c5, *c6, *c7, *c8, *c9, *c10, *c11, *c12; + + double *Pp, *P; // Legendre polynomial P and first derivative Pp + double *Bzplus, *Brplus; // used in Bz and Br mag. field calc. + // Remote source constant variables: + double* Brem1; // 1-dim. remote source constants Brem1[n], n=2,...,nmax. + double* rorem; // remote convergence radius for coil i: rorem[i] + double** Brem; // remote source constant for coil i: Brem[i][n] (coil index i, source constant index n) + double* z0remG; // remote source point for group G: z0remG[g] + double* roremG; // remote convergence radius for group G: roremG[g] + double** + BremG; // remote source constant for symmetry group g: BremG[g][n] (group index g, source constant index n) + // Magnetic charge source constant variables: + double** + Vrem; // magnetic charge remote source constant for coil i: Vrem[i][n] (coil index i, source const. index n) + // Central source constant variables: + double* Bcen1; // 1-dim. central source constants Bcen1[n], n=0,...,nmax. + int* Nsp; // number of central source points for coil i: Nsp[i] (i=0,...,Ncoil-1) + double** z0cen; // central source point local z value for coil i: z0cen[i][j] (coil index i, source point index j) + double** rocen; // central convergence radius for coil i: rocen[i][j] (coil index i, source point index j) + double*** + Bcen; // central source constant for coil i: Bcen[i][j][n] (coil index i, source point index j, source constant index n) + int* NspG; // number of central source points for group g: Nsp[g] (g=0,...,Ng-1) + double** z0cenG; // central source point local z value: z0cenG[g][j] (group index g, source point index j) + double** rocenG; // central convergence radius: rocenG[g][j] (group index g, source point index j) + double*** + BcenG; // central source constant: BcenG[g][j][n] (group index g, source point index j, source const. index n) + int *jlast, *jlastG; // last central source point index for coil and group calculation + double rclimit; }; @@ -129,59 +134,59 @@ double rclimit; inline double MagfieldCoils::pow2(double x) { - return x*x; + return x * x; } //////////////////////////////////////////////////////// inline double MagfieldCoils::FMIN(double a, double b) { - return ((a)<=(b)?(a):(b)); + return ((a) <= (b) ? (a) : (b)); } //////////////////////////////////////////////////////// inline double MagfieldCoils::FMAX(double a, double b) { - return ((a)>(b)?(a):(b)); + return ((a) > (b) ? (a) : (b)); } //////////////////////////////////////////////////////// inline double MagfieldCoils::FMIN3(double a, double b, double c) { - return (FMIN(a,b)<=(c)?(FMIN(a,b)):(c)); + return (FMIN(a, b) <= (c) ? (FMIN(a, b)) : (c)); } //////////////////////////////////////////////////////// inline double MagfieldCoils::FMAX3(double a, double b, double c) { - return (FMAX(a,b)>(c)?(FMAX(a,b)):(c)); + return (FMAX(a, b) > (c) ? (FMAX(a, b)) : (c)); } //////////////////////////////////////////////////////// inline double MagfieldCoils::Mag(double* vec) { - return sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]); + return sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]); } //////////////////////////////////////////////////////// inline double MagfieldCoils::Mag(double x, double y, double z) { - return sqrt(x*x+y*y+z*z); + return sqrt(x * x + y * y + z * z); } //////////////////////////////////////////////////////// inline void MagfieldCoils::SetTrackingStart() { - for(int i=0; i> Ncoil; - if (Ncoil<1) - { - puts("Message from function MagfieldCoils::CoilRead:"); - puts("Ncoil is not positive !"); - puts("Program running is stopped !!! "); - exit(1); - } -// Dynamic memory allocation for the coil parameters: - coil=new double*[Ncoil]; - for(int i=0; i> cu >> Cx >> Cy >> Cz >> alpha >> beta >> tu >> L >> Rmin >> Rmax ; -// Coil parameter determination: - coil[i][0]=cu; coil[i][1]=Cx; coil[i][2]=Cy; coil[i][3]=Cz; - coil[i][4]=alpha; coil[i][5]=beta; coil[i][6]=tu; - coil[i][7]=L; coil[i][8]=Rmin; coil[i][9]=Rmax; - coil[i][10]=cu*tu/(L*(Rmax-Rmin)); // current density - coil[i][11]=sin(beta/180.*M_PI)*sin(alpha/180.*M_PI); // coil direction unit vector comp. ux - coil[i][12]=-sin(beta/180.*M_PI)*cos(alpha/180.*M_PI); // coil direction unit vector comp. uy - coil[i][13]=cos(beta/180.*M_PI); // coil direction unit vector comp. uz - } - input.close(); -// -// Test of coil parameters: - for(int i=0; i=coil[i][9]) - { - printf("Message from function MagfieldCoils::CoilRead1:\ + input >> Ncoil; + if (Ncoil < 1) { + puts("Message from function MagfieldCoils::CoilRead:"); + puts("Ncoil is not positive !"); + puts("Program running is stopped !!! "); + exit(1); + } + // Dynamic memory allocation for the coil parameters: + coil = new double*[Ncoil]; + for (int i = 0; i < Ncoil; i++) { + coil[i] = new double[14]; + } + // Reading the coil parameters: + double cu, Cx, Cy, Cz, alpha, beta, tu, L, Rmin, Rmax; + for (int i = 0; i < Ncoil; i++) { + input >> cu >> Cx >> Cy >> Cz >> alpha >> beta >> tu >> L >> Rmin >> Rmax; + // Coil parameter determination: + coil[i][0] = cu; + coil[i][1] = Cx; + coil[i][2] = Cy; + coil[i][3] = Cz; + coil[i][4] = alpha; + coil[i][5] = beta; + coil[i][6] = tu; + coil[i][7] = L; + coil[i][8] = Rmin; + coil[i][9] = Rmax; + coil[i][10] = cu * tu / (L * (Rmax - Rmin)); // current density + coil[i][11] = sin(beta / 180. * M_PI) * sin(alpha / 180. * M_PI); // coil direction unit vector comp. ux + coil[i][12] = -sin(beta / 180. * M_PI) * cos(alpha / 180. * M_PI); // coil direction unit vector comp. uy + coil[i][13] = cos(beta / 180. * M_PI); // coil direction unit vector comp. uz + } + input.close(); + // + // Test of coil parameters: + for (int i = 0; i < Ncoil; i++) { + if (coil[i][6] <= 1.e-9 || coil[i][7] <= 1.e-9 || coil[i][8] <= 1.e-9 || coil[i][9] < 1.e-9) { + printf("Message from function MagfieldCoils::CoilRead1:\ + non-positive coil parameters: tu, Rmin, Rmax or L !!!\ + Program running is stopped !!! i= %9i \n", + i); + exit(1); + } + if (coil[i][8] >= coil[i][9]) { + printf("Message from function MagfieldCoils::CoilRead1:\ Rmin>=Rmax!!!\ - Program running is stopped !!! i= %9i \n",i); - exit(1); - } - } - -//------------------------------------------------------ -// -// Coil symmetry group calculations: -// - G=new int[Ncoil]; // coil with index i is in symmetry group g=G[i] - for(int i=0; i alpha, betac-> beta; -// If uc is opposite to u: change of current sign - for(int k=0; k<3; k++) - vec[k]=u[k]-uc[k]; - if(Mag(vec)>0.5) - coil[ic][0]*=-1.; - coil[ic][4]=alpha; coil[ic][5]=beta; - coil[ic][11]=sin(beta/180.*M_PI)*sin(alpha/180.*M_PI); // coil direction unit vector comp. ux - coil[ic][12]=-sin(beta/180.*M_PI)*cos(alpha/180.*M_PI); // coil direction unit vector comp. uy - coil[ic][13]=cos(beta/180.*M_PI); // coil direction unit vector comp. uz -// After that, all coils within a symmetry group have the same coil directions (alpha, beta). - } // end of epstol-if - } // end of symmetry group loop ic - } // end of symmetry group if - } // end of coil index loop i -// Nc, Cin, Line and Z array calculation: - Nc=new int[Ng]; - for(int g=0; g alpha, betac-> beta; + // If uc is opposite to u: change of current sign + for (int k = 0; k < 3; k++) + vec[k] = u[k] - uc[k]; + if (Mag(vec) > 0.5) + coil[ic][0] *= -1.; + coil[ic][4] = alpha; + coil[ic][5] = beta; + coil[ic][11] = + sin(beta / 180. * M_PI) * sin(alpha / 180. * M_PI); // coil direction unit vector comp. ux + coil[ic][12] = + -sin(beta / 180. * M_PI) * cos(alpha / 180. * M_PI); // coil direction unit vector comp. uy + coil[ic][13] = cos(beta / 180. * M_PI); // coil direction unit vector comp. uz + // After that, all coils within a symmetry group have the same coil directions (alpha, beta). + } // end of epstol-if + } // end of symmetry group loop ic + } // end of symmetry group if + } // end of coil index loop i + // Nc, Cin, Line and Z array calculation: + Nc = new int[Ng]; + for (int g = 0; g < Ng; g++) + Nc[g] = Nctemp[g]; // number of coils in symmetry group g + Cin = new int*[Ng]; + for (int g = 0; g < Ng; g++) { + Cin[g] = new int[Nc[g]]; + for (int c = 0; c < Nc[g]; c++) + Cin[g][c] = Cintemp[g][c]; // coil indices i in symmetry group g + } + Line = new double*[Ng]; + for (int g = 0; g < Ng; g++) { + Line[g] = new double[6]; + for (int k = 0; k < 6; k++) + Line[g][k] = Linetemp[g][k]; // line point and line direction unit vector + } + Z = new double*[Ng]; + for (int g = 0; g < Ng; g++) { + Z[g] = new double[Nc[g]]; + for (int c = 0; c < Nc[g]; c++) + Z[g][c] = Ztemp[g][c]; // local z coordinate of coil center in symmetry group g; Z[g][0]=0. + } + // Delete Nctemp, Cintemp, Linetemp, Ztemp: + delete[] Nctemp; + for (int i = 0; i < Ncoil; i++) delete[] Cintemp[i]; - delete[] Cintemp; - for(int i=0; i>Ncoil >> Ng; - if(Ncoil<1 || Ng<1) - { - printf("Message from function MagfieldCoils::CoilRead2:\ + // This function reads the coil parameters from the data file dirname+objectname+coilgroup.txt. + string filename = dirname + objectname + "coilgroup.txt"; + ifstream input; + input.open(filename.c_str()); + if (!input.is_open()) { + puts("Message from function CoilRead2 of class Magfield:"); + puts("Cannot open the coil.txt file!"); + puts("Program running is stopped !!! "); + exit(1); + } + // We read the coil parameters: + input >> Ncoil >> Ng; + if (Ncoil < 1 || Ng < 1) { + printf("Message from function MagfieldCoils::CoilRead2:\ Ncoil<1 or Ng<1 !!!\ - Computation is stopped !!! Ncoil, Ng= %9i %9i \n",Ncoil,Ng); - exit(1); - } - input >> Nelliptic >> nmax; -// Dynamic memory allocation for the coil and symmetry group parameters: - coil=new double*[Ncoil]; - for(int i=0; i> ix ; - input >> coil[i][0] >> coil[i][1] >> coil[i][2] >> coil[i][3] ; - input >> coil[i][4] >> coil[i][5] >> coil[i][6] ; - input >> coil[i][7] >> coil[i][8] >> coil[i][9] >> coil[i][10] ; - input >> coil[i][11] >> coil[i][12] >> coil[i][13] ; - } -// Dynamic memory allocation for the symmetry group arrays: - Nc=new int[Ng]; - Line=new double*[Ng]; - Cin=new int*[Ng]; - Z=new double*[Ng]; -// We read the symmetry group parameters: - int gx; - for(int g=0; g> gx >> Nc[g] ; -// Cin[g][c]: // global coil index corresponding to local coil index c in group g -// Z[g][c]: local z coordinate of coil center in symmetry group g; c=0,..., Nc[g]-1; Z[g][0]=0 -// Symmetry axis line point: Line[g][0], Line[g][1], Line[g][2]; line direction unit vector: Line[g][3], Line[g][4], Line[g][5]; - Line[g]= new double[6]; - input >> Line[g][0] >> Line[g][1] >> Line[g][2] ; - input >> Line[g][3] >> Line[g][4] >> Line[g][5] ; - Cin[g]= new int[Nc[g]]; - Z[g]= new double[Nc[g]]; - for(int c=0; c> Cin[g][c] >> Z[g][c] ; - } - - - input.close(); - - G=new int[Ncoil]; - Brem1=new double[nmax+1]; - Bcen1=new double[nmax+1]; - - return; + Computation is stopped !!! Ncoil, Ng= %9i %9i \n", + Ncoil, + Ng); + exit(1); + } + input >> Nelliptic >> nmax; + // Dynamic memory allocation for the coil and symmetry group parameters: + coil = new double*[Ncoil]; + for (int i = 0; i < Ncoil; i++) { + coil[i] = new double[14]; + } + // We read the coil parameters: + int ix; + for (int i = 0; i < Ncoil; i++) { + input >> ix; + input >> coil[i][0] >> coil[i][1] >> coil[i][2] >> coil[i][3]; + input >> coil[i][4] >> coil[i][5] >> coil[i][6]; + input >> coil[i][7] >> coil[i][8] >> coil[i][9] >> coil[i][10]; + input >> coil[i][11] >> coil[i][12] >> coil[i][13]; + } + // Dynamic memory allocation for the symmetry group arrays: + Nc = new int[Ng]; + Line = new double*[Ng]; + Cin = new int*[Ng]; + Z = new double*[Ng]; + // We read the symmetry group parameters: + int gx; + for (int g = 0; g < Ng; g++) { + input >> gx >> Nc[g]; + // Cin[g][c]: // global coil index corresponding to local coil index c in group g + // Z[g][c]: local z coordinate of coil center in symmetry group g; c=0,..., Nc[g]-1; Z[g][0]=0 + // Symmetry axis line point: Line[g][0], Line[g][1], Line[g][2]; line direction unit vector: Line[g][3], Line[g][4], Line[g][5]; + Line[g] = new double[6]; + input >> Line[g][0] >> Line[g][1] >> Line[g][2]; + input >> Line[g][3] >> Line[g][4] >> Line[g][5]; + Cin[g] = new int[Nc[g]]; + Z[g] = new double[Nc[g]]; + for (int c = 0; c < Nc[g]; c++) + input >> Cin[g][c] >> Z[g][c]; + } + + + input.close(); + + G = new int[Ncoil]; + Brem1 = new double[nmax + 1]; + Bcen1 = new double[nmax + 1]; + + return; } @@ -456,126 +484,117 @@ void MagfieldCoils::CoilGroupRead() void MagfieldCoils::MagsourceRead() { -// This function reads the source constants from the data files -// dirname+objectname+magsource_central.txt and dirname+objectname+magsource_remote.txt. -// Part 1: reading the remote source constants: - string filename=dirname+objectname+"magsource_remote.txt"; - ifstream input; - input.open(filename.c_str()); - if (!input.is_open()) - { - puts("Message from function MagfieldCoils::MagsourceRead:"); - puts("Cannot open the input source file magsource_remote.txt!"); - puts("Program running is stopped !!! "); - exit(1); - } -// Dynamic memory allocations (rorem, Brem): - rorem=new double[Ncoil]; - Brem=new double*[Ncoil]; - for(int i=0; i> ix >> rorem[i] ; - for(int n=2; n<=nmax; n++) - input >> nx >>Brem[i][n] ; - } -// Dynamic memory allocations (Vrem): - Vrem=new double*[Ncoil]; - for(int i=0; i> ix ; - for(int n=0; n<=nmax; n++) - input >> nx >>Vrem[i][n] ; - } -// Dynamic memory allocations (z0remG, roremG, BremG): - z0remG=new double[Ng]; - roremG=new double[Ng]; - BremG=new double*[Ng]; - for(int g=0; g> gx >> Ncg >> z0remG[g] >> roremG[g] ; - for(int n=2; n<=nmax; n++) - input >> nx >>BremG[g][n] ; - } - input.close(); -// -// Part 2: reading the central source coefficients: - filename=dirname+objectname+"magsource_central.txt"; - input.open(filename.c_str()); - if (!input.is_open()) - { - puts("Message from function MagfieldCoils::MagsourceRead:"); - puts("Cannot open the input source file magsource_central.txt!"); - puts("Program running is stopped !!! "); - exit(1); - } - int jx; -// Dynamic memory allocations (Nsp, z0cen, rocen, Bcen): - Nsp=new int[Ncoil]; - z0cen=new double*[Ncoil]; - rocen=new double*[Ncoil]; - Bcen=new double**[Ncoil]; -// Coil index loop: - string text; -// input >> text; - for(int i=0; i> ix >> Nsp[i]; -// Source point loop: - z0cen[i]= new double[Nsp[i]]; - rocen[i]= new double[Nsp[i]]; - Bcen[i]=new double*[Nsp[i]]; - for(int j=0; j> jx >> z0cen[i][j] >> rocen[i][j] ; - Bcen[i][j]=new double[nmax+1]; - for(int n=0; n<=nmax; n++) - input >> nx >>Bcen[i][j][n] ; - } - } -// Dynamic memory allocations (NspG, z0cenG, rocenG, BcenG): - NspG=new int[Ng]; - z0cenG=new double*[Ng]; - rocenG=new double*[Ng]; - BcenG=new double**[Ng]; -// Group index loop: -// input >> text; - for(int g=0; g> gx >> NspG[g]; -// Source point loop: - z0cenG[g]= new double[NspG[g]]; - rocenG[g]= new double[NspG[g]]; - BcenG[g]=new double*[NspG[g]]; - for(int j=0; j> jx >> z0cenG[g][j] >> rocenG[g][j] ; - BcenG[g][j]=new double[nmax+1]; - for(int n=0; n<=nmax; n++) - input >> nx >>BcenG[g][j][n] ; - } - } - input.close(); + // This function reads the source constants from the data files + // dirname+objectname+magsource_central.txt and dirname+objectname+magsource_remote.txt. + // Part 1: reading the remote source constants: + string filename = dirname + objectname + "magsource_remote.txt"; + ifstream input; + input.open(filename.c_str()); + if (!input.is_open()) { + puts("Message from function MagfieldCoils::MagsourceRead:"); + puts("Cannot open the input source file magsource_remote.txt!"); + puts("Program running is stopped !!! "); + exit(1); + } + // Dynamic memory allocations (rorem, Brem): + rorem = new double[Ncoil]; + Brem = new double*[Ncoil]; + for (int i = 0; i < Ncoil; i++) + Brem[i] = new double[nmax + 1]; + // Coil index loop: + int ix, nx; + for (int i = 0; i < Ncoil; i++) { + input >> ix >> rorem[i]; + for (int n = 2; n <= nmax; n++) + input >> nx >> Brem[i][n]; + } + // Dynamic memory allocations (Vrem): + Vrem = new double*[Ncoil]; + for (int i = 0; i < Ncoil; i++) + Vrem[i] = new double[nmax + 1]; + for (int i = 0; i < Ncoil; i++) { + input >> ix; + for (int n = 0; n <= nmax; n++) + input >> nx >> Vrem[i][n]; + } + // Dynamic memory allocations (z0remG, roremG, BremG): + z0remG = new double[Ng]; + roremG = new double[Ng]; + BremG = new double*[Ng]; + for (int g = 0; g < Ng; g++) + BremG[g] = new double[nmax + 1]; + // Group index loop: + int gx, Ncg; + for (int g = 0; g < Ng; g++) { + if (Nc[g] == 1) + continue; + input >> gx >> Ncg >> z0remG[g] >> roremG[g]; + for (int n = 2; n <= nmax; n++) + input >> nx >> BremG[g][n]; + } + input.close(); + // + // Part 2: reading the central source coefficients: + filename = dirname + objectname + "magsource_central.txt"; + input.open(filename.c_str()); + if (!input.is_open()) { + puts("Message from function MagfieldCoils::MagsourceRead:"); + puts("Cannot open the input source file magsource_central.txt!"); + puts("Program running is stopped !!! "); + exit(1); + } + int jx; + // Dynamic memory allocations (Nsp, z0cen, rocen, Bcen): + Nsp = new int[Ncoil]; + z0cen = new double*[Ncoil]; + rocen = new double*[Ncoil]; + Bcen = new double**[Ncoil]; + // Coil index loop: + string text; + // input >> text; + for (int i = 0; i < Ncoil; i++) { + input >> ix >> Nsp[i]; + // Source point loop: + z0cen[i] = new double[Nsp[i]]; + rocen[i] = new double[Nsp[i]]; + Bcen[i] = new double*[Nsp[i]]; + for (int j = 0; j < Nsp[i]; j++) { + input >> jx >> z0cen[i][j] >> rocen[i][j]; + Bcen[i][j] = new double[nmax + 1]; + for (int n = 0; n <= nmax; n++) + input >> nx >> Bcen[i][j][n]; + } + } + // Dynamic memory allocations (NspG, z0cenG, rocenG, BcenG): + NspG = new int[Ng]; + z0cenG = new double*[Ng]; + rocenG = new double*[Ng]; + BcenG = new double**[Ng]; + // Group index loop: + // input >> text; + for (int g = 0; g < Ng; g++) { + if (Nc[g] == 1) // no group calculation if the group has only 1 coil + { + z0cenG[g] = new double[1]; + rocenG[g] = new double[1]; + NspG[g] = 1; + BcenG[g] = new double*[1]; + BcenG[g][0] = new double[nmax + 1]; + continue; + } + input >> gx >> NspG[g]; + // Source point loop: + z0cenG[g] = new double[NspG[g]]; + rocenG[g] = new double[NspG[g]]; + BcenG[g] = new double*[NspG[g]]; + for (int j = 0; j < NspG[g]; j++) { + input >> jx >> z0cenG[g][j] >> rocenG[g][j]; + BcenG[g][j] = new double[nmax + 1]; + for (int n = 0; n <= nmax; n++) + input >> nx >> BcenG[g][j][n]; + } + } + input.close(); } @@ -583,64 +602,69 @@ void MagfieldCoils::MagsourceRead() void MagfieldCoils::DynamicMemoryAllocations() { -// This function makes dynamic memory allocations for the variables -// C0, C1, c0 to c12, P, Pp, Bzplus, Brplus. -// It also computes the values C0[n], C1[n], c0[n] -- c12[n] (n=2,...,nmax). - C0=new double[nmax+1]; C1=new double[nmax+1]; - c1=new double[nmax+1]; c2=new double[nmax+1]; c3=new double[nmax+1]; - c4=new double[nmax+1]; c5=new double[nmax+1]; c6=new double[nmax+1]; - c7=new double[nmax+1]; c8=new double[nmax+1]; c9=new double[nmax+1]; - c10=new double[nmax+1]; c11=new double[nmax+1]; c12=new double[nmax+1]; - - C0[1]=1.; - for(int n=2; n<=nmax; n++) - { - C0[n]=1./(1.*n); - C1[n]=1./(1.*(n-1.)); - c1[n]=(2.*n-1.)/(1.*n); - c2[n]=(n-1.)/(1.*n); - c3[n]=(2.*n-1.)/(1.*(n-1.)); - c4[n]=(1.*n)/(1.*(n-1.)); - c5[n]=(1.)/(n*1.); - c6[n]=(1.)/(n+1.); - if(n>=4) - { - int m=n-2; - double M=(m+1.)*(m+2.)*(2.*m-1.); - double Mp=(m)*(m+1.)*(2.*m-1.); - double A=(2.*m-1.)*(2.*m+1.)*(2.*m+3.); - double Ap=A; - double B=2.*m*m*(2.*m+3.)-1.; - double Bp=2.*m*(m+2.)*(2.*m-1.)-3.; - double C=m*(m-1.)*(2.*m+3.); - double Cp=m*(m+1.)*(2.*m+3.); - c7[n]=A/M; - c8[n]=B/M; - c9[n]=C/M; - c10[n]=Ap/Mp; - c11[n]=Bp/Mp; - c12[n]=Cp/Mp; - } - } -// Dynamic memory for P, Pp, Bzplus, Brplus, jlast, jlastG: - P=new double[nmax+2]; - Pp=new double[nmax+2]; - Bzplus=new double[nmax+1]; - Brplus=new double[nmax+1]; -// - jlast=new int[Ncoil]; - for(int i=0; i= 4) { + int m = n - 2; + double M = (m + 1.) * (m + 2.) * (2. * m - 1.); + double Mp = (m) * (m + 1.) * (2. * m - 1.); + double A = (2. * m - 1.) * (2. * m + 1.) * (2. * m + 3.); + double Ap = A; + double B = 2. * m * m * (2. * m + 3.) - 1.; + double Bp = 2. * m * (m + 2.) * (2. * m - 1.) - 3.; + double C = m * (m - 1.) * (2. * m + 3.); + double Cp = m * (m + 1.) * (2. * m + 3.); + c7[n] = A / M; + c8[n] = B / M; + c9[n] = C / M; + c10[n] = Ap / Mp; + c11[n] = Bp / Mp; + c12[n] = Cp / Mp; + } + } + // Dynamic memory for P, Pp, Bzplus, Brplus, jlast, jlastG: + P = new double[nmax + 2]; + Pp = new double[nmax + 2]; + Bzplus = new double[nmax + 1]; + Brplus = new double[nmax + 1]; + // + jlast = new int[Ncoil]; + for (int i = 0; i < Ncoil; i++) + jlast[i] = -2; + jlastG = new int[Ng]; + for (int g = 0; g < Ng; g++) + jlastG[g] = -2; } ///////////////////////////////////////////////////// - - // #include "MagfieldElliptic.cc" //////////////////////////////////////////////////////////////////////////////////////////// @@ -652,187 +676,186 @@ void MagfieldCoils::DynamicMemoryAllocations() void MagfieldCoils::Magfield2EllipticCoil(int i, double z, double r, double& Bz, double& Br) { -// This function computes the magnetic field components Bz and Br -// of an axially symmetric coil index i, with z axis as symmetry axis, -// at a fieldpoint with (z,r) cylindrical coordinates relative to coil center and direction, using -// the first, second and third complete elliptic integrals. - double L, Zmin, Zmax, Rmin, Rmax, sigma; - double R, Z, delr2, sumr2, delz2, eta, d, K, EK, PIK, S; - double sign, delRr, Rlow[2], Rhigh[2]; - const double mu0=4.*M_PI*1.e-7; - double x[2][1001], w[2][1001]; // nodes and weights -// Coil parameters: - L=coil[i][7]; // coil length - Zmin=-L/2.; Zmax=L/2.; // coil endpoints relative to coil center - Rmin=coil[i][8]; Rmax=coil[i][9]; // coil inner and outer radii - sigma=coil[i][10]; // current density -// Field point test: - if(fabs(z-Zmin)<1.e-8 && r>=Rmin-(Rmax-Rmin)*1.e-8 && - r<=Rmax+(Rmax-Rmin)*1.e-8) - z=Zmin-1.e-8; - if(fabs(z-Zmax)<1.e-8 && r>=Rmin-(Rmax-Rmin)*1.e-8 && - r<=Rmax+(Rmax-Rmin)*1.e-8) - z=Zmax+1.e-8; -// Radial numerical integration node definition: - int N[2]; - N[0]=N[1]=Nelliptic; -// -// R-integration limits: - int M; - if(z>Zmin && zRmin && r= Rmin - (Rmax - Rmin) * 1.e-8 && r <= Rmax + (Rmax - Rmin) * 1.e-8) + z = Zmin - 1.e-8; + if (fabs(z - Zmax) < 1.e-8 && r >= Rmin - (Rmax - Rmin) * 1.e-8 && r <= Rmax + (Rmax - Rmin) * 1.e-8) + z = Zmax + 1.e-8; + // Radial numerical integration node definition: + int N[2]; + N[0] = N[1] = Nelliptic; + // + // R-integration limits: + int M; + if (z > Zmin && z < Zmax && r > Rmin && r < Rmax) // field point inside the coil winding: integral in 2 steps + { + M = 2; + Rlow[0] = Rmin; + Rhigh[0] = r - (r - Rmin) * 1.e-12; + Rlow[1] = r + (Rmax - r) * 1.e-12; + Rhigh[1] = Rmax; + GaussLegendreIntegration(N[0], Rlow[0], Rhigh[0], x[0], w[0]); + GaussLegendreIntegration(N[1], Rlow[1], Rhigh[1], x[1], w[1]); + } + else { + M = 1; + Rlow[0] = Rmin; + Rhigh[0] = Rmax; + GaussLegendreIntegration(N[0], Rlow[0], Rhigh[0], x[0], w[0]); + } + // Integration: + Bz = Br = 0.; + for (int m = 0; m < M; m++) // integral step m-loop + { + for (int iR = 0; iR < N[m]; iR++) // radius R-loop { - if(iZ==0) - { Z=Zmax; sign=1.; } - else - { Z=Zmin; sign=-1.; } - delr2=(r-R)*(r-R); - delz2=(z-Z)*(z-Z); - sumr2=(r+R)*(r+R); - d=delr2/sumr2; - eta=(delr2+delz2)/(sumr2+delz2); - S=sqrt(sumr2+delz2); - K=RF_Carlson(0., eta, 1.); - EK=-1./3.*RD_Carlson(0., eta, 1.); - delRr=R-r; - if(d<1.e-18) - { - d=1.e-18; - if(R>r) - delRr=(r+R)*1.e-9; - else if(R r) + delRr = (r + R) * 1.e-9; + else if (R < r) + delRr = -(r + R) * 1.e-9; + } + PIK = 1. / 3. * RJ_Carlson(0., eta, 1., d); + Bzhat += sign * R * (z - Z) / (S * (R + r)) * (K + delRr / (2. * R) * PIK * (1. - d)); + Brhat += sign * R / S * (2. * EK + K); + } // end of Z-loop + Bz += w[m][iR] * Bzhat; + Br += w[m][iR] * Brhat; + // cout.precision(4); + // cout << scientific << m << " " << iR << " " << x[m][iR] << " " << w[m][iR] << endl; + } // end of radius R-loop + } // end of integral step m-loop + double C = -mu0 / M_PI * sigma; + Bz *= C; + Br *= C; + return; } /////////////////////////////////////////////////// -void MagfieldCoils::MagfieldEllipticCoil(int i, const double *P, double *B) +void MagfieldCoils::MagfieldEllipticCoil(int i, const double* P, double* B) { -// This function computes the magnetic field components B[0],B[1],B[2] -// at a field point P[0],P[1],P[2], due to the 3-dimensional coil -// with index i, using the first, second and third complete elliptic -// integrals. SI units are used (P[k] in m, B[k] in T, k=0,1,2)! -// The coil is defined by the coil parameters coil[i][j], j=0,...,10. - double z, r, Bz, Br, u[3], Ploc[3], Pr[3], w[3], C[3]; -// We define a local coordinate system: -// origo at the coil center, local z axis parallel to coil axis, in u vector direction. -// Coil direction u and center C: - u[0]=coil[i][11]; - u[1]=coil[i][12]; - u[2]=coil[i][13]; - for(int k=0;k<=2;k++) - C[k]=coil[i][1+k]; // coil center -// Local z and r coordinates of the field point P: - for(int k=0;k<=2;k++) - Ploc[k]=P[k]-C[k]; - z=Ploc[1]*u[1]+Ploc[2]*u[2]+Ploc[0]*u[0]; - for(int k=0;k<=2;k++) - Pr[k]=Ploc[k]-z*u[k]; - r=sqrt(Pr[1]*Pr[1]+Pr[2]*Pr[2]+Pr[0]*Pr[0]); -// Bz and Br calculation: - Magfield2EllipticCoil(i, z, r, Bz, Br); -// B[k] calculation from Bz, Br: - if(r<1.e-12) - for(int k=0; k<=2; k++) - B[k]=Bz*u[k]; - else - { - for(int k=0;k<=2;k++) - w[k]=Pr[k]/r; - for(int k=0;k<=2;k++) - B[k]=Bz*u[k]+Br*w[k]; - } - return; + // This function computes the magnetic field components B[0],B[1],B[2] + // at a field point P[0],P[1],P[2], due to the 3-dimensional coil + // with index i, using the first, second and third complete elliptic + // integrals. SI units are used (P[k] in m, B[k] in T, k=0,1,2)! + // The coil is defined by the coil parameters coil[i][j], j=0,...,10. + double z, r, Bz, Br, u[3], Ploc[3], Pr[3], w[3], C[3]; + // We define a local coordinate system: + // origo at the coil center, local z axis parallel to coil axis, in u vector direction. + // Coil direction u and center C: + u[0] = coil[i][11]; + u[1] = coil[i][12]; + u[2] = coil[i][13]; + for (int k = 0; k <= 2; k++) + C[k] = coil[i][1 + k]; // coil center + // Local z and r coordinates of the field point P: + for (int k = 0; k <= 2; k++) + Ploc[k] = P[k] - C[k]; + z = Ploc[1] * u[1] + Ploc[2] * u[2] + Ploc[0] * u[0]; + for (int k = 0; k <= 2; k++) + Pr[k] = Ploc[k] - z * u[k]; + r = sqrt(Pr[1] * Pr[1] + Pr[2] * Pr[2] + Pr[0] * Pr[0]); + // Bz and Br calculation: + Magfield2EllipticCoil(i, z, r, Bz, Br); + // B[k] calculation from Bz, Br: + if (r < 1.e-12) + for (int k = 0; k <= 2; k++) + B[k] = Bz * u[k]; + else { + for (int k = 0; k <= 2; k++) + w[k] = Pr[k] / r; + for (int k = 0; k <= 2; k++) + B[k] = Bz * u[k] + Br * w[k]; + } + return; } /////////////////////////////////////////////////////////////////////////////// -void MagfieldCoils::MagfieldElliptic(const double *P, double *B) +void MagfieldCoils::MagfieldElliptic(const double* P, double* B) { -// This function computes the magnetic field components B[0],B[1],B[2] -// at a field point P[0],P[1],P[2], due to all coils, using elliptic integrals. -// SI units are used (P[k] in m, B[k] in T, k=0,1,2 --> x, y, z components)! - double Bi[3]; - for(int k=0; k<=2; k++) - B[k]=0.; - for(int i=0; i x, y, z components)! + double Bi[3]; + for (int k = 0; k <= 2; k++) + B[k] = 0.; + for (int i = 0; i < Ncoil; i++) { + MagfieldEllipticCoil(i, P, Bi); + for (int k = 0; k <= 2; k++) + B[k] += Bi[k]; + } + return; +} // #include "Magfield.cc" /////////////////////////////////////////////////////////////////////////////////////////// -bool MagfieldCoils::Magfield(const double *P, double *B) +bool MagfieldCoils::Magfield(const double* P, double* B) { -// This function computes the magnetic field components B[0],B[1],B[2] -// at a field point P[0],P[1],P[2], due to all symmetry groups (coils), -// using central or remote zonal harmonic (Legendre polynomial) expansion, or elliptic integrals. -// SI units are used (P[k] in m, B[k] in T, k=0,1,2 --> x,y,z components)! -// The return value is true, if only zonal harmonic expansions are used for the magnetic field calc.; -// if elliptic integral calc. is used (even for 1 coil): the return value is false. - double Bgroup[3]; - bool magfield=true; - for(int k=0; k<3; k++) - B[k]=0.; - for(int g=0; g x,y,z components)! + // The return value is true, if only zonal harmonic expansions are used for the magnetic field calc.; + // if elliptic integral calc. is used (even for 1 coil): the return value is false. + double Bgroup[3]; + bool magfield = true; + for (int k = 0; k < 3; k++) + B[k] = 0.; + for (int g = 0; g < Ng; g++) { + bool magfieldgroup = MagfieldGroup(g, P, Bgroup); + if (magfieldgroup == false) + magfield = false; + for (int k = 0; k < 3; k++) + B[k] += Bgroup[k]; + } + return magfield; } /////////////////////////////////////////////////////////////////////////////////////// - - // #include "MagsourceRemote.cc" @@ -840,18 +863,20 @@ bool MagfieldCoils::Magfield(const double *P, double *B) // Remote convergence radius calculation: double MagfieldCoils::Funrorem(int i, double z0) -// This function computes the remote convergence radius -// Funrorem=rorem for coil i, at the coil axis source point z0 +// This function computes the remote convergence radius +// Funrorem=rorem for coil i, at the coil axis source point z0 // (defined relative to the coil center; positive z0 in coil direction u, which is determined by alpha and beta). // rorem = maximum distance of the axis point z0 from the coil winding (outer corner points). { - double L, Rmax, Zmin, Zmax; - L=coil[i][7]; Rmax=coil[i][9]; - Zmin=-L/2.; Zmax=L/2.; // coil endpoint Z values defined relative to the coil center - if(z0>0.) - return sqrt((z0-Zmin)*(z0-Zmin)+Rmax*Rmax); - else - return sqrt((z0-Zmax)*(z0-Zmax)+Rmax*Rmax); + double L, Rmax, Zmin, Zmax; + L = coil[i][7]; + Rmax = coil[i][9]; + Zmin = -L / 2.; + Zmax = L / 2.; // coil endpoint Z values defined relative to the coil center + if (z0 > 0.) + return sqrt((z0 - Zmin) * (z0 - Zmin) + Rmax * Rmax); + else + return sqrt((z0 - Zmax) * (z0 - Zmax) + Rmax * Rmax); } ///////////////////////////////////////////////////// @@ -868,56 +893,57 @@ void MagfieldCoils::Magsource2RemoteCoil(int i, double z0, double rorem) // N is small if the(Rmax-Rmin)/Rmin ratio is small, and it is large // if this ratio is large. { - const double mu0=4.*M_PI*1.e-7; - double x[1001], w[1001]; // nodes and weights - double L, sigma, Zmin, Zmax, Rmin, Rmax; -// Coil parameters: - L=coil[i][7]; // coil length - Zmin=-L/2.; Zmax=L/2.; // coil endpoints relative to coil center - Rmin=coil[i][8]; Rmax=coil[i][9]; // coil inner and outer radii - sigma=coil[i][10]; // current density -// Radial integration number N: - int N; - double ratio=(Rmax-Rmin)/Rmin; - if(ratio<0.2) - N=32; - else - N=32*ratio/0.2; - if(N>1000) - N=1000; -// Initialization of Brem1[n] and Pp[n]: - for(int n=2; n<=nmax; n++) - Brem1[n]=0.; - Pp[0]=0.; Pp[1]=1.; -// C: - double C=mu0*sigma/(2.*rorem*rorem); -// Radial integration nodes and weights: - GaussLegendreIntegration(N, Rmin, Rmax, x, w); -// Zminmax, sign: - double Zminmax[2]={Zmin, Zmax}; - int sign[2]={-1, 1}; -// Gauss-Legendre integration loop: - for(int integ=0; integ 1000) + N = 1000; + // Initialization of Brem1[n] and Pp[n]: + for (int n = 2; n <= nmax; n++) + Brem1[n] = 0.; + Pp[0] = 0.; + Pp[1] = 1.; + // C: + double C = mu0 * sigma / (2. * rorem * rorem); + // Radial integration nodes and weights: + GaussLegendreIntegration(N, Rmin, Rmax, x, w); + // Zminmax, sign: + double Zminmax[2] = {Zmin, Zmax}; + int sign[2] = {-1, 1}; + // Gauss-Legendre integration loop: + for (int integ = 0; integ < N; integ++) { + double R = x[integ]; + double R2 = R * R; + for (int iz = 0; iz <= 1; iz++) { + double Z = Zminmax[iz]; + double ros = sqrt(R2 + (Z - z0) * (Z - z0)); + double us = (Z - z0) / ros; + double a = ros / rorem; + double b = a; + double d = C * R2 * sign[iz] * w[integ]; + for (int n = 2; n <= nmax; n++) { + double usPp = us * Pp[n - 1]; + Pp[n] = 2. * usPp - Pp[n - 2] + (usPp - Pp[n - 2]) * C1[n]; + Brem1[n] += d * c6[n] * b * Pp[n]; + b *= a; + } + } + } + return; } ////////////////////////////////////////////////////////////////////////////////////// @@ -926,7 +952,7 @@ void MagfieldCoils::Magsource2RemoteCoil(int i, double z0, double rorem) // (3-dim. case) void MagfieldCoils::MagsourceRemoteCoils() -// This function computes the Brem[i][n] (i=0,...,Ncoil-1; n=2,...,nmax) +// This function computes the Brem[i][n] (i=0,...,Ncoil-1; n=2,...,nmax) // remote source constants for all coils. // Number of remote source points for a fixed coil is 1 (coil center). // The results are written into the data file dirname+objectname+magsource_remote.txt: @@ -935,46 +961,44 @@ void MagfieldCoils::MagsourceRemoteCoils() // i: coil index, // rorem[i]: remote convergence radius for coil i. { -// Output to file dirname+objectname+magsource_remote.txt : - string filename=dirname+objectname+"magsource_remote.txt"; - ofstream output; - output.precision(16); - output.open(filename.c_str()); -// Dynamic memory allocations (Brem1, rorem, Brem): - Brem1=new double[nmax+1]; - rorem=new double[Ncoil]; - Brem=new double*[Ncoil]; - for(int i=0; i z0remG[g] (in group Z-coordinate system): - double z0rem, zmin=1.e9, zmax=-1.e9; - for(int c=0; c z0remG[g] (in group Z-coordinate system): + double z0rem, zmin = 1.e9, zmax = -1.e9; + for (int c = 0; c < Nc[g]; c++) // c: local coil index in group + { + int i = Cin[g][c]; + double L = coil[i][7]; + double zA = Z[g][c] - L / 2.; + double zB = Z[g][c] + L / 2.; // coil edges if (zA < zmin) zmin = zA; if (zB > zmax) zmax = zB; - } - z0rem=z0remG[g]=(zmin+zmax)/2.; // group center in group Z-coordinate system -// Remote convergence radius of group --> roremG[g] : - double rorem=0.; - for(int c=0; crorem) rorem=roremc; - } - roremG[g]=rorem; + } + z0rem = z0remG[g] = (zmin + zmax) / 2.; // group center in group Z-coordinate system + // Remote convergence radius of group --> roremG[g] : + double rorem = 0.; + for (int c = 0; c < Nc[g]; c++) // c: local coil index in group + { + int i = Cin[g][c]; + double roremc = Funrorem(i, z0rem - Z[g][c]); // remote source point relative to coil center + if (roremc > rorem) + rorem = roremc; + } + roremG[g] = rorem; } - - - // #include "MagsourceMagcharge.cc" ///////////////////////////////////////////////////// @@ -1071,104 +1093,106 @@ void MagfieldCoils::RemoteSourcepointGroup(int g) // Magnetic charge remote source constant calculation for all coils void MagfieldCoils::MagsourceMagchargeCoils() -// This function computes the Vrem[i][n] (i=0,...,Ncoil-1; n=0,...,nmax) +// This function computes the Vrem[i][n] (i=0,...,Ncoil-1; n=0,...,nmax) // magnetic charge remote source constants for all coils. // The results are written into the data file dirname+objectname+magsource_remote.txt: // Ncoil: number of coils, // nmax: number of source coefficients for a fixed coil is nmax+1 (n=0,...,nmax), // i: coil index, n source constant index. { -// Output to file dirname+objectname+magsource_remote.txt : - double Rmin, Rmax, sigma; - string filename=dirname+objectname+"magsource_remote.txt"; - ofstream output; - output.precision(16); - output.open(filename.c_str(), ios::app); -// Dynamic memory allocations (Vrem): - Vrem=new double*[Ncoil]; - for(int i=0; i0 && n%2==0) // even positive n - P[n]=-(n-1.)/double(n)*P[n-2]; - if(n%2==1) - Vrem[i][n]=0.; - else - Vrem[i][n]=sigma*P[n]*Rmax*Rmax/(2.*(n+2.)*(n+3.))*(1.-ration3); - if(n%2==1) - output << scientific << setprecision(1) << setw(15) << n << setw(28) < 0 && n % 2 == 0) // even positive n + P[n] = -(n - 1.) / double(n) * P[n - 2]; + if (n % 2 == 1) + Vrem[i][n] = 0.; + else + Vrem[i][n] = sigma * P[n] * Rmax * Rmax / (2. * (n + 2.) * (n + 3.)) * (1. - ration3); + if (n % 2 == 1) + output << scientific << setprecision(1) << setw(15) << n << setw(28) << Vrem[i][n] << endl; + else + output << scientific << setprecision(16) << setw(15) << n << setw(28) << Vrem[i][n] << endl; + if (ration3 < 1.e-20) + ration3 = 0.; + else + ration3 *= ratio; + } + } + output.close(); + return; } - - // #include "MagsourceCentral.cc" // Central convergence radius calculation: double MagfieldCoils::Funrocen(int i, double z0) -// This function computes the effective central convergence radius +// This function computes the effective central convergence radius // Funrocen=rocen for coil i, at the axis source point z0 // (defined relative to the coil center; positive z0 in coil direction u). // rocen = minimum distance of the axis point z0 from the coil inner corner points // (Zmin,Rmin) and (Zmax,Rmin). { - double L, Zmin, Zmax,Rmin; - L=coil[i][7]; Rmin=coil[i][8]; - Zmin=-L/2.; Zmax=L/2.; // coil endpoint Z values defined relative to the coil center - if(z0<0.) - return sqrt((z0-Zmin)*(z0-Zmin)+Rmin*Rmin); - else - return sqrt((z0-Zmax)*(z0-Zmax)+Rmin*Rmin); + double L, Zmin, Zmax, Rmin; + L = coil[i][7]; + Rmin = coil[i][8]; + Zmin = -L / 2.; + Zmax = L / 2.; // coil endpoint Z values defined relative to the coil center + if (z0 < 0.) + return sqrt((z0 - Zmin) * (z0 - Zmin) + Rmin * Rmin); + else + return sqrt((z0 - Zmax) * (z0 - Zmax) + Rmin * Rmin); } ///////////////////////////////////////////////////// double MagfieldCoils::FunrocenG(int g, double z0) -// This function computes the central convergence radius +// This function computes the central convergence radius // FunrocenG=rocenmin for the symmetry group g, at the axis source point z0 // (defined in the group Z coordinate system). // rocenmin = minimum distance of the axis point z0 from the windings of the group coils. { - double rocenmin=1.e9; - for(int c=0; c=Zmax) - rocen=sqrt((Z0-Zmax)*(Z0-Zmax)+Rmin*Rmin); + double rocenmin = 1.e9; + for (int c = 0; c < Nc[g]; c++) // c: local coil index in group + { + int i = Cin[g][c]; // i: global coil index + double L, Zmin, Zmax, Rmin, Z0; + L = coil[i][7]; + Rmin = coil[i][8]; + Zmin = -L / 2.; + Zmax = L / 2.; // coil endpoint Z values defined relative to the coil center + Z0 = z0 - Z[g][c]; // source point relative to coil center + double rocen; + if (Z0 <= Zmin) + rocen = sqrt((Z0 - Zmin) * (Z0 - Zmin) + Rmin * Rmin); + else if (Z0 >= Zmax) + rocen = sqrt((Z0 - Zmax) * (Z0 - Zmax) + Rmin * Rmin); else - rocen=Rmin; - if(rocen1000) - N=1000; -// Initialization of Bcen1[n] and Pp[n]: - for(int n=0; n<=nmax; n++) - Bcen1[n]=0.; - Pp[0]=0.; Pp[1]=1.; -// C: - double C=-0.5*mu0*sigma*Rocen; -// Radial integration nodes and weights: - GaussLegendreIntegration(N, Rmin, Rmax, x, w); -// Zminmax, sign: - double Zminmax[2]={Zmin, Zmax}; - const int sign[2]={-1, 1}; -// Gauss-Legendre integration loop: - for(int integ=0; integ 1000) + N = 1000; + // Initialization of Bcen1[n] and Pp[n]: + for (int n = 0; n <= nmax; n++) + Bcen1[n] = 0.; + Pp[0] = 0.; + Pp[1] = 1.; + // C: + double C = -0.5 * mu0 * sigma * Rocen; + // Radial integration nodes and weights: + GaussLegendreIntegration(N, Rmin, Rmax, x, w); + // Zminmax, sign: + double Zminmax[2] = {Zmin, Zmax}; + const int sign[2] = {-1, 1}; + // Gauss-Legendre integration loop: + for (int integ = 0; integ < N; integ++) { + double R = x[integ]; + double R2 = R * R; + for (int iz = 0; iz <= 1; iz++) { + double Z = Zminmax[iz]; + double ros = sqrt(R2 + (Z - z0) * (Z - z0)); + double us = (Z - z0) / ros; + double a = Rocen / ros; + double b = a; + double d = 1. / (ros * ros * ros); + double D = C * R2 * sign[iz] * d * w[integ]; + Bcen1[0] += 0.5 * mu0 * sigma * sign[iz] * us * w[integ]; + Bcen1[1] += D; + for (int n = 2; n <= nmax; n++) { + double usPp = us * Pp[n - 1]; + Pp[n] = 2. * usPp - Pp[n - 2] + (usPp - Pp[n - 2]) * C1[n]; + Bcen1[n] += D * C0[n] * b * Pp[n]; + b *= a; + } + } + } + return; } ////////////////////////////////////////////////////////////////////////////////////// @@ -1261,46 +1286,43 @@ void MagfieldCoils::MagsourceCentralCoils() // (defined relative to the coil center; positive z0 in coil direction u), // rocen: central convergence radius for a fixed source point and coil. { -// Output to file dirname+objectname+magsource_central.txt : - string filename=dirname+objectname+"magsource_central.txt"; - ofstream output; - output.precision(16); - output.open(filename.c_str()); -// Dynamic memory allocations (Bcen1, Nsp, z0cen, rocen, Bcen): - Bcen1=new double[nmax+1]; - Nsp=new int[Ncoil]; - z0cen=new double*[Ncoil]; - rocen=new double*[Ncoil]; - Bcen=new double**[Ncoil]; -// Coil index loop: -// output << "coils" << endl; - for(int i=0; i Z0cen; Z0cen.push_back(0.); double z0; - double roremC=rorem[i]; - double z0max=2.5*roremC; - int nsp=0; // number of central source points on positive z side - int k=0; - do - { - k+=1; - double rocen1=Funrocen(i, Z0cen[k-1]); - double del1=rocen1/4.; - z0=Z0cen[k-1]+del1; - rocen1=Funrocen(i, z0); - double del2=rocen1/4.; - double del=FMIN(del1, del2); - z0=Z0cen[k-1]+del; - Z0cen.push_back(z0); - if(z0 Z0cen; double z0; - double roremC=roremG[g]; - double z0min=z0remG[g]-2.*roremC; - double z0max=z0remG[g]+2.*roremC; - Z0cen.push_back(z0min); // j=0 - int nsp=1; // number of central source points in group g - int j=0; // --> z0min - do - { - j+=1; - double rocen1=FunrocenG(g, Z0cen[j-1]); - double del1=rocen1/4.; - z0=Z0cen[j-1]+del1; - rocen1=FunrocenG(g, z0); - double del2=rocen1/4.; - double del=FMIN(del1, del2); - z0=Z0cen[j-1]+del; - Z0cen.push_back(z0); - if(z0 z0min + do { + j += 1; + double rocen1 = FunrocenG(g, Z0cen[j - 1]); + double del1 = rocen1 / 4.; + z0 = Z0cen[j - 1] + del1; + rocen1 = FunrocenG(g, z0); + double del2 = rocen1 / 4.; + double del = FMIN(del1, del2); + z0 = Z0cen[j - 1] + del; + Z0cen.push_back(z0); + if (z0 < z0max) + nsp += 1; + } while (z0 < z0max); + // Number of central source points for group g: + NspG[g] = nsp; + // Central source points and conv. radii for group g (source point loop): + z0cenG[g] = new double[nsp]; + rocenG[g] = new double[nsp]; + // + for (int j = 0; j < nsp; j++) { + z0 = z0cenG[g][j] = Z0cen[j]; + rocenG[g][j] = FunrocenG(g, z0); } } - // #include "MagfieldRemote.cc" @@ -1473,10 +1485,10 @@ void MagfieldCoils::CentralSourcepointsGroup(int g) // Magnetic field calculation with remote zonal harmonic expansion // (local axisymmetric case) -bool MagfieldCoils::Magfield2Remote(bool bcoil, int ig, double z, double r, double& Bz,double& Br,double& rc) +bool MagfieldCoils::Magfield2Remote(bool bcoil, int ig, double z, double r, double& Bz, double& Br, double& rc) // This function computes the magnetic field components Bz and Br at fieldpoint z and r by remote expansion, // either for a coil or for a symmetry group. -// It computes also the remote convergence ratio rc=Rorem/ro. +// It computes also the remote convergence ratio rc=Rorem/ro. // bcoil=true: coil calculation; bcoil=false: symmetry group calculation. // ig: coil index for bcoil=true, and group index for bcoil=false. // bcoil=true: Brem[ig][n], n=2,...,nmax remote source constants for coils are used. @@ -1491,118 +1503,120 @@ bool MagfieldCoils::Magfield2Remote(bool bcoil, int ig, double z, double r, doub // // z0, Rorem, ro, u, sr, rc, rcn: { - double z0, Rorem; // source point and conv. radius - if(bcoil==true) - { z0=0.; Rorem=rorem[ig]; } - else - { z0=z0remG[ig]; Rorem=roremG[ig]; } - double delz=z-z0; - double ro=sqrt(r*r+delz*delz); - if(ro<1.e-9) // field point is very close to the source point --> rc>1 ! - ro=1.e-9; - double u=delz/ro; - double sr=r/ro; - rc=Rorem/ro; // convergence ratio - double rcn=rc*rc; -// If rc>rclimit: the Legendre polynomial series is too slow, or not convergent (for rc>1) !!! - if(rc>rclimit) - { - return false; - } -// First 2 terms of Legendre polynomial P and its derivative Pp (P-primed) - P[0]=1.; P[1]=u; - Pp[0]=0.; Pp[1]=1.; -// - bool even=false; // true: only the even n terms are used - if(bcoil==false) // group calculation - { - double sum=0., sumodd=0.; - for(int n=2; n<=12; n++) - sum+=fabs(BremG[ig][n]); - for(int n=3; n<=12; n+=2) - sumodd+=fabs(BremG[ig][n]); - if(sumodd5) - { - double Beps=1.e-15*(fabs(Bz)+fabs(Br)); - double Bdelta=fabs(Bzplus[n])+fabs(Brplus[n])+fabs(Bzplus[n-1])+fabs(Brplus[n-1])+ - fabs(Bzplus[n-2])+fabs(Brplus[n-2])+fabs(Bzplus[n-3])+fabs(Brplus[n-3]); - if(Bdelta rc>1 ! + ro = 1.e-9; + double u = delz / ro; + double sr = r / ro; + rc = Rorem / ro; // convergence ratio + double rcn = rc * rc; + // If rc>rclimit: the Legendre polynomial series is too slow, or not convergent (for rc>1) !!! + if (rc > rclimit) { + return false; + } + // First 2 terms of Legendre polynomial P and its derivative Pp (P-primed) + P[0] = 1.; + P[1] = u; + Pp[0] = 0.; + Pp[1] = 1.; + // + bool even = false; // true: only the even n terms are used + if (bcoil == false) // group calculation + { + double sum = 0., sumodd = 0.; + for (int n = 2; n <= 12; n++) + sum += fabs(BremG[ig][n]); + for (int n = 3; n <= 12; n += 2) + sumodd += fabs(BremG[ig][n]); + if (sumodd < sum * 1.e-8) + even = true; + } + // We start here the series expansion: + Bz = Br = 0.; + int nlast = 0; + double brem; + if (bcoil == true || even == true) // even calculation + { + // Only the even terms are needed in the series + // (because for the middle source point in the coil centre + // Brem[ig][n]=0 for odd n) + for (int n = 2; n <= 4; n++) { + if (bcoil == true) + brem = Brem[ig][n]; + else + brem = BremG[ig][n]; + // P[n]=c1[n]*u*P[n-1]-c2[n]*P[n-2]; + // Pp[n]=c3[n]*u*Pp[n-1]-c4[n]*Pp[n-2]; + double uPp = u * Pp[n - 1]; + double uP = u * P[n - 1]; + P[n] = 2. * uP - P[n - 2] - (uP - P[n - 2]) * C0[n]; + Pp[n] = 2. * uPp - Pp[n - 2] + (uPp - Pp[n - 2]) * C1[n]; + rcn = rcn * rc; + Bzplus[n] = brem * rcn * P[n]; + Brplus[n] = sr * brem * c5[n] * rcn * Pp[n]; + Bz += Bzplus[n]; + Br += Brplus[n]; + } + double u2 = u * u; + double rc2 = rc * rc; + for (int n = 6; n <= nmax - 1; n += 2) { + if (bcoil == true) + brem = Brem[ig][n]; + else + brem = BremG[ig][n]; + nlast = n; + P[n] = (c7[n] * u2 - c8[n]) * P[n - 2] - c9[n] * P[n - 4]; + Pp[n] = (c10[n] * u2 - c11[n]) * Pp[n - 2] - c12[n] * Pp[n - 4]; + rcn = rcn * rc2; + Bzplus[n] = brem * rcn * P[n]; + Brplus[n] = sr * brem * c5[n] * rcn * Pp[n]; + Bz += Bzplus[n]; + Br += Brplus[n]; + double Beps = 1.e-15 * (fabs(Bz) + fabs(Br)); + double Bdelta = fabs(Bzplus[n]) + fabs(Brplus[n]) + fabs(Bzplus[n - 2]) + fabs(Brplus[n - 2]); + if (Bdelta < Beps || Bdelta < 1.e-20) + break; + } + } + else // even + odd group calculation: both even and odd terms are needed in the series + { + for (int n = 2; n <= nmax - 1; n++) { + brem = BremG[ig][n]; + nlast = n; + P[n] = c1[n] * u * P[n - 1] - c2[n] * P[n - 2]; + Pp[n] = c3[n] * u * Pp[n - 1] - c4[n] * Pp[n - 2]; + rcn = rcn * rc; + Bzplus[n] = brem * rcn * P[n]; + Brplus[n] = sr * brem * c5[n] * rcn * Pp[n]; + Bz += Bzplus[n]; + Br += Brplus[n]; + if (n > 5) { + double Beps = 1.e-15 * (fabs(Bz) + fabs(Br)); + double Bdelta = fabs(Bzplus[n]) + fabs(Brplus[n]) + fabs(Bzplus[n - 1]) + fabs(Brplus[n - 1]) + + fabs(Bzplus[n - 2]) + fabs(Brplus[n - 2]) + fabs(Bzplus[n - 3]) + fabs(Brplus[n - 3]); + if (Bdelta < Beps || Bdelta < 1.e-20) + break; + } } - } - } - if(nlast>=nmax-2) - return false; - else - return true; + } + if (nlast >= nmax - 2) + return false; + else + return true; } - - ///////////////////////////////////////////////// - - - - // #include "MagfieldMagcharge.cc" @@ -1614,130 +1628,125 @@ bool MagfieldCoils::Magfield2Remote(bool bcoil, int ig, double z, double r, doub bool MagfieldCoils::Magfield2Magcharge(int i, double z, double r, double& Bz, double& Br, double& rc) { -// This function computes the magnetic field components Bz and Br at fieldpoint z and r by remote -// magnetic charge expansion, for coil i. -// It computes also the remote convergence ratio rc. -// The Vrem[i][n], n=0,...,nmax magnetic charge remote source constants for coil i are used. -// z: it is relative to the coil center (z=0). -// nmax: maximal index of the source constants (maximum of n). -// SI units are used ! -// If the convergence of the Legendre-series expansion is too slow: -// rc>rclimit, and the return value is false -// (in this case one should use another computation method !!!) -// - const double mu0=4.*M_PI*1.e-7; - double Hzmax, Hrmax, Hzmin, Hrmin, Hz, Hr, Mz, rcmin, rcmax; - double L, Zmin, Zmax, Rmin, Rmax, sigma; - L=coil[i][7]; Rmin=coil[i][8]; Rmax=coil[i][9]; // coil length, inner and outer radius - sigma=coil[i][10]; // current density - Zmin=-L/2.; Zmax=L/2.; // coil endpoints relative to coil center - bool hfieldmax=Hfield(i, z-Zmax, r, Hzmax, Hrmax, rcmax); - bool hfieldmin=Hfield(i, z-Zmin, r, Hzmin, Hrmin, rcmin); - rc=FMAX(rcmin, rcmax); -// If rc>rclimit: the Legendre polynomial series is too slow, or not convergent (for rc>1) !!! - if(hfieldmax==false || hfieldmin==false) - { - return false; - } - Hz=Hzmax-Hzmin; - Hr=Hrmax-Hrmin; - if(z>=Zmin && z<=Zmax && r<=Rmin) - Mz=sigma*(Rmax-Rmin); - else if(z>=Zmin && z<=Zmax && r>Rmin && r<=Rmax) - Mz=sigma*(Rmax-r); - else - Mz=0.; - Bz=mu0*(Mz+Hz); - Br=mu0*Hr; -// cout.precision(15); -// cout << scientific <=rclimit: the Legendre polynomial series is too slow, or not convergent (for rc>1) !!! - if(rc>=rclimit) - { - return false; - } -// First 2 terms of Legendre polynomial P and its derivative Pp (P-primed) - P[0]=1.; P[1]=u; - Pp[0]=0.; Pp[1]=1.; -// First 2 terms of the series: - if(bcoil==true) - { - Bz=Bcen[ig][j][0]+Bcen[ig][j][1]*rc*u; - Br=-sr*Bcen[ig][j][1]/2.*rc; - } - else - { - Bz=BcenG[ig][j][0]+BcenG[ig][j][1]*rc*u; - Br=-sr*BcenG[ig][j][1]/2.*rc; - } - int nlast=0; - if(rc<1.e-10) goto label; -// -// We start here the central series expansion: - double bcen; - for(int n=2; n<=nmax-1; n++) - { - if(bcoil==true) - bcen=Bcen[ig][j][n]; - else - bcen=BcenG[ig][j][n]; - nlast=n; - rcn*=rc; - double uPp=u*Pp[n-1]; double uP=u*P[n-1]; - P[n]=2.*uP-P[n-2] - (uP - P[n-2])*C0[n]; - Pp[n]=2.*uPp-Pp[n-2] +(uPp - Pp[n-2])*C1[n]; -// P[n]=c1[n]*u*P[n-1]-c2[n]*P[n-2]; -// Pp[n]=c3[n]*u*Pp[n-1]-c4[n]*Pp[n-2]; - Bzplus[n]=bcen*rcn*P[n]; - Brplus[n]=-sr*bcen*c6[n]*rcn*Pp[n]; - Bz+=Bzplus[n]; Br+=Brplus[n]; - if(n>5) - { - double Beps=1.e-15*(fabs(Bz)+fabs(Br)); - double Bdelta=fabs(Bzplus[n])+fabs(Brplus[n])+fabs(Bzplus[n-1])+fabs(Brplus[n-1])+ - fabs(Bzplus[n-2])+fabs(Brplus[n-2])+fabs(Bzplus[n-3])+fabs(Brplus[n-3]); - if(Bdelta=nmax-1) - { - rc=1.; - return false; - } -label: ; -// Effective central correction for Bz: - if(bcoil==true) - { - double Lhalf=coil[ig][7]*0.5; - if(z>-Lhalf && zRmax) - Bz+=-mu0*sigma*(Rmax-Rmin); - else if(r>Rmin && r=rclimit: the Legendre polynomial series is too slow, or not convergent (for rc>1) !!! + if (rc >= rclimit) { + return false; + } + // First 2 terms of Legendre polynomial P and its derivative Pp (P-primed) + P[0] = 1.; + P[1] = u; + Pp[0] = 0.; + Pp[1] = 1.; + // First 2 terms of the series: + if (bcoil == true) { + Bz = Bcen[ig][j][0] + Bcen[ig][j][1] * rc * u; + Br = -sr * Bcen[ig][j][1] / 2. * rc; + } + else { + Bz = BcenG[ig][j][0] + BcenG[ig][j][1] * rc * u; + Br = -sr * BcenG[ig][j][1] / 2. * rc; + } + int nlast = 0; + if (rc < 1.e-10) + goto label; + // + // We start here the central series expansion: + double bcen; + for (int n = 2; n <= nmax - 1; n++) { + if (bcoil == true) + bcen = Bcen[ig][j][n]; + else + bcen = BcenG[ig][j][n]; + nlast = n; + rcn *= rc; + double uPp = u * Pp[n - 1]; + double uP = u * P[n - 1]; + P[n] = 2. * uP - P[n - 2] - (uP - P[n - 2]) * C0[n]; + Pp[n] = 2. * uPp - Pp[n - 2] + (uPp - Pp[n - 2]) * C1[n]; + // P[n]=c1[n]*u*P[n-1]-c2[n]*P[n-2]; + // Pp[n]=c3[n]*u*Pp[n-1]-c4[n]*Pp[n-2]; + Bzplus[n] = bcen * rcn * P[n]; + Brplus[n] = -sr * bcen * c6[n] * rcn * Pp[n]; + Bz += Bzplus[n]; + Br += Brplus[n]; + if (n > 5) { + double Beps = 1.e-15 * (fabs(Bz) + fabs(Br)); + double Bdelta = fabs(Bzplus[n]) + fabs(Brplus[n]) + fabs(Bzplus[n - 1]) + fabs(Brplus[n - 1]) + + fabs(Bzplus[n - 2]) + fabs(Brplus[n - 2]) + fabs(Bzplus[n - 3]) + fabs(Brplus[n - 3]); + if (Bdelta < Beps || Bdelta < 1.e-20) + break; + } + } + if (nlast >= nmax - 1) { + rc = 1.; + return false; + } +label:; + // Effective central correction for Bz: + if (bcoil == true) { + double Lhalf = coil[ig][7] * 0.5; + if (z > -Lhalf && z < Lhalf) { + double Rmin = coil[ig][8], Rmax = coil[ig][9]; + double sigma = coil[ig][10]; // current density + if (r > Rmax) + Bz += -mu0 * sigma * (Rmax - Rmin); + else if (r > Rmin && r < Rmax) + Bz += -mu0 * sigma * (r - Rmin); + } + } + return true; } - - ///////////////////////////////////////////////// - - - // #include "MagfieldCoil.cc" -bool MagfieldCoils::MagfieldCoil(int i, const double *P, double *B) +bool MagfieldCoils::MagfieldCoil(int i, const double* P, double* B) { -// This function computes the magnetic field components B[0],B[1],B[2] -// at a field point P[0],P[1],P[2], due to the 3-dimensional coil -// with index i, using central, remote or magnetic charge zonal harmonic expansions, -// or elliptic integrals. -// SI units are used (P[k] in m, B[k] in T, k=0,1,2 --> components x, y, z)! -// The coil is defined by the coil parameters coil[i][j], j=0,...,13. -// The return value is true, if only zonal harmonic expansions are used for the magnetic field calc.; -// if elliptic integral calc. is used: the return value is false. + // This function computes the magnetic field components B[0],B[1],B[2] + // at a field point P[0],P[1],P[2], due to the 3-dimensional coil + // with index i, using central, remote or magnetic charge zonal harmonic expansions, + // or elliptic integrals. + // SI units are used (P[k] in m, B[k] in T, k=0,1,2 --> components x, y, z)! + // The coil is defined by the coil parameters coil[i][j], j=0,...,13. + // The return value is true, if only zonal harmonic expansions are used for the magnetic field calc.; + // if elliptic integral calc. is used: the return value is false. + // + // We start now the magnetic field calculation. + // First we define the local coordinate system of the coil, and the + // z and r local coordinates of the field point P. + // Local coordinate system: + // origo at coil center, + // local z axis parallel to coil axis, in u vector direction. + bool magfieldcoil = true; + double C[3], u[3]; + for (int k = 0; k < 3; k++) + C[k] = coil[i][1 + k]; // coil center + u[0] = coil[i][11]; // coil direction unit vector component ux + u[1] = coil[i][12]; // coil direction unit vector component uy + u[2] = coil[i][13]; // coil direction unit vector component uz + // Local cylindrical z and r coordinates of the field point P (relative to coil center and direction): + double Ploc[3], Pr[3], w[3]; + for (int k = 0; k < 3; k++) + Ploc[k] = P[k] - C[k]; + double z = Ploc[1] * u[1] + Ploc[2] * u[2] + Ploc[0] * u[0]; + for (int k = 0; k < 3; k++) + Pr[k] = Ploc[k] - z * u[k]; + double r = sqrt(Pr[1] * Pr[1] + Pr[2] * Pr[2] + Pr[0] * Pr[0]); + double r2 = r * r; + // + // Starting now Bz and Br calculation. + bool sps; + int jbest; + double rcbest; + // --------------------------- + // Step 1. Probably most of the coils are far from the field point, therefore + // we try first the remote series expansion + // (with remote source point at the coil center) + double ro2 = r2 + z * z; + double rorem2 = rorem[i] * rorem[i]; + double Bz, Br, rc; + if (ro2 > rorem2 * 2.) { + bool magfield = Magfield2Remote(true, i, z, r, Bz, Br, rc); + if (magfield == true) + goto labelend; + } + // Step 2. We try the central zonal method in the beginning of tracking, + // by searching all central source points. + if (jlast[i] == -2) { + jlast[i] = -1; + sps = SourcePointSearching(true, i, z, r, 0, jbest, rcbest); + if (sps == true) { + bool magfield = Magfield2Central(true, i, jbest, z, r, Bz, Br, rc); + if (magfield == true) { + magfieldcoil = magfield; + goto labelend; + } + } + } + // --------------------------- + // Step 3. The field point is close to the coil. We try next the central + // Legendre polynomial expansion. If some source point for this coil + // has been already used (jlast[i]>-1), we search the central source + // point j close to the old source point jlast[i]. + sps = false; + if (jlast[i] > -1) + sps = SourcePointSearching(true, i, z, r, 1, jbest, rcbest); + if (jlast[i] > -1 && sps == false) + sps = SourcePointSearching(true, i, z, r, 2, jbest, rcbest); + if (sps == true) { + bool magfield = Magfield2Central(true, i, jbest, z, r, Bz, Br, rc); + if (magfield == true) + goto labelend; + } + // --------------------------- + // Step 4. We try again the remote series expansion + // (with remote source point at the coil center) + if (ro2 > rorem2 * 1.0204) { + bool magfield = Magfield2Remote(true, i, z, r, Bz, Br, rc); + if (magfield == true) + goto labelend; + } + // --------------------------- + // --------------------------- + // Step 5. We try now the magnetic charge method + { + double L = coil[i][7]; + double Rmax = coil[i][9]; // coil length and outer radius + double Zmin = -L / 2.; + double Zmax = L / 2.; // coil endpoints relative to coil center + double romax2 = r2 + (z - Zmax) * (z - Zmax); + double romin2 = r2 + (z - Zmin) * (z - Zmin); + rorem2 = Rmax * Rmax; + if (romax2 > rorem2 * 1.0204 && romin2 > rorem2 * 1.0204) { + bool magfield = Magfield2Magcharge(i, z, r, Bz, Br, rc); + if (magfield == true) + goto labelend; + } + } + // --------------------------- + // Step 6. We try again the central zonal method, by searching now all central source points. + sps = SourcePointSearching(true, i, z, r, 0, jbest, rcbest); + if (sps == true) { + bool magfield = Magfield2Central(true, i, jbest, z, r, Bz, Br, rc); + if (magfield == true) { + magfieldcoil = magfield; + goto labelend; + } + } + // ----------------------------- + // Step 7. Unfortunately, no appropriate central, remote or magnetic charge expansion + // was found. We have to use elliptic integrals: + // cout << "step 6 " << endl; + Magfield2EllipticCoil(i, z, r, Bz, Br); + magfieldcoil = false; // -// We start now the magnetic field calculation. -// First we define the local coordinate system of the coil, and the -// z and r local coordinates of the field point P. -// Local coordinate system: -// origo at coil center, -// local z axis parallel to coil axis, in u vector direction. - bool magfieldcoil=true; - double C[3], u[3]; - for(int k=0; k<3; k++) - C[k]=coil[i][1+k]; // coil center - u[0]=coil[i][11]; // coil direction unit vector component ux - u[1]=coil[i][12]; // coil direction unit vector component uy - u[2]=coil[i][13]; // coil direction unit vector component uz -// Local cylindrical z and r coordinates of the field point P (relative to coil center and direction): - double Ploc[3], Pr[3], w[3]; - for(int k=0; k<3; k++) - Ploc[k]=P[k]-C[k]; - double z=Ploc[1]*u[1]+Ploc[2]*u[2]+Ploc[0]*u[0]; - for(int k=0; k<3; k++) - Pr[k]=Ploc[k]-z*u[k]; - double r=sqrt(Pr[1]*Pr[1]+Pr[2]*Pr[2]+Pr[0]*Pr[0]); - double r2=r*r; -// -// Starting now Bz and Br calculation. - bool sps; - int jbest; double rcbest; -// --------------------------- -// Step 1. Probably most of the coils are far from the field point, therefore -// we try first the remote series expansion -// (with remote source point at the coil center) - double ro2=r2+z*z; - double rorem2=rorem[i]*rorem[i]; - double Bz, Br, rc; - if(ro2>rorem2*2.) - { - bool magfield=Magfield2Remote(true, i, z, r, Bz, Br, rc); - if(magfield==true) - goto labelend; - } -// Step 2. We try the central zonal method in the beginning of tracking, -// by searching all central source points. - if(jlast[i]==-2) - { - jlast[i]=-1; - sps=SourcePointSearching(true, i, z, r, 0, jbest, rcbest); - if(sps==true) - { - bool magfield=Magfield2Central(true, i, jbest, z, r, Bz, Br, rc); - if(magfield==true) - { - magfieldcoil=magfield; - goto labelend; - } - } - } -// --------------------------- -// Step 3. The field point is close to the coil. We try next the central -// Legendre polynomial expansion. If some source point for this coil -// has been already used (jlast[i]>-1), we search the central source -// point j close to the old source point jlast[i]. - sps=false; - if(jlast[i]>-1) - sps=SourcePointSearching(true, i, z, r, 1, jbest, rcbest); - if(jlast[i]>-1 && sps==false) - sps=SourcePointSearching(true, i, z, r, 2, jbest, rcbest); - if(sps==true) - { - bool magfield=Magfield2Central(true, i, jbest, z, r, Bz, Br, rc); - if(magfield==true) - goto labelend; - } -// --------------------------- -// Step 4. We try again the remote series expansion -// (with remote source point at the coil center) - if(ro2>rorem2*1.0204) - { - bool magfield=Magfield2Remote(true, i, z, r, Bz, Br, rc); - if(magfield==true) - goto labelend; - } -// --------------------------- -// --------------------------- -// Step 5. We try now the magnetic charge method - { - double L=coil[i][7]; double Rmax=coil[i][9]; // coil length and outer radius - double Zmin=-L/2.; double Zmax=L/2.; // coil endpoints relative to coil center - double romax2=r2+(z-Zmax)*(z-Zmax); - double romin2=r2+(z-Zmin)*(z-Zmin); - rorem2=Rmax*Rmax; - if(romax2>rorem2*1.0204 && romin2>rorem2*1.0204) - { - bool magfield=Magfield2Magcharge(i, z, r, Bz, Br, rc); - if(magfield==true) - goto labelend; - } - } -// --------------------------- -// Step 6. We try again the central zonal method, by searching now all central source points. - sps=SourcePointSearching(true, i, z, r, 0, jbest, rcbest); - if(sps==true) - { - bool magfield=Magfield2Central(true, i, jbest, z, r, Bz, Br, rc); - if(magfield==true) - { - magfieldcoil=magfield; - goto labelend; - } - } -// ----------------------------- -// Step 7. Unfortunately, no appropriate central, remote or magnetic charge expansion -// was found. We have to use elliptic integrals: -// cout << "step 6 " << endl; - Magfield2EllipticCoil(i, z, r, Bz, Br); - magfieldcoil=false; -// // B[k],k=0,1,2 calculation from Bz, Br: -labelend: ; - if(r<1.e-15 || fabs(Br) components x, y, z)! -// The return value is true, if only zonal harmonic expansions are used for the magnetic field calc.; -// if elliptic integral calc. is also used: the return value is false. -// If the group has only 1 coil --> coil calculation: - if(Nc[g]==1) - return MagfieldCoil(Cin[g][0], P, B); + // This function computes the magnetic field components B[0],B[1],B[2] + // at a field point P[0],P[1],P[2], due to the symmetry group + // with index g, using central or remote zonal harmonic expansions, + // or single coil calculations. + // SI units are used (P[k] in m, B[k] in T, k=0,1,2 --> components x, y, z)! + // The return value is true, if only zonal harmonic expansions are used for the magnetic field calc.; + // if elliptic integral calc. is also used: the return value is false. + // If the group has only 1 coil --> coil calculation: + if (Nc[g] == 1) + return MagfieldCoil(Cin[g][0], P, B); + // + // We start now the magnetic field calculation. + // First we define the local coordinate system of the group, and the + // z and r local coordinates of the field point P. + // Local coordinate system: + // origo at Z[g][0], local z axis parallel to group axis, in u vector direction. + bool magfieldgroup = true; + double C[3], u[3], w[3]; + for (int k = 0; k < 3; k++) + C[k] = Line[g][k]; // group origo + for (int k = 0; k < 3; k++) + u[k] = Line[g][3 + k]; // group direction unit vector + // Local cylindrical z and r coordinates of the field point P (relative to group origo and direction): + double Ploc[3], Pr[3]; + for (int k = 0; k < 3; k++) + Ploc[k] = P[k] - C[k]; + double z = Ploc[1] * u[1] + Ploc[2] * u[2] + Ploc[0] * u[0]; + for (int k = 0; k < 3; k++) + Pr[k] = Ploc[k] - z * u[k]; + double r = sqrt(Pr[1] * Pr[1] + Pr[2] * Pr[2] + Pr[0] * Pr[0]); + double r2 = r * r; + // + // Starting now Bz and Br calculation. + bool sps; + int jbest; + double rcbest; + // --------------------------- + // Step 1. We try first the remote series expansion + // (with remote source point at the group center) + double ro2 = r2 + (z - z0remG[g]) * (z - z0remG[g]); + double rorem2 = roremG[g] * roremG[g]; + double Bz, Br, rc; + if (ro2 > rorem2 * 2.) { + bool magfield = Magfield2Remote(false, g, z, r, Bz, Br, rc); + if (magfield == true) + goto labelend; + } + // --------------------------- + // Step 2. We try the central zonal method in the beginning of tracking, + // by searching all central source points. + if (jlastG[g] == -2) { + jlastG[g] = -1; + sps = SourcePointSearching(false, g, z, r, 0, jbest, rcbest); + if (sps == true) { + bool magfield = Magfield2Central(false, g, jbest, z, r, Bz, Br, rc); + if (magfield == true) + goto labelend; + } + } + // --------------------------- + // Step 3. The field point is close to the group. We try next the central + // Legendre polynomial expansion. If some source point for this group + // has been already used (jlastG[g]>-1), we search the central source + // point j close to the old source point jlastG[g]. + sps = false; + if (jlastG[g] > -1) + sps = SourcePointSearching(false, g, z, r, 1, jbest, rcbest); + if (jlastG[g] > -1 && sps == false) + sps = SourcePointSearching(false, g, z, r, 2, jbest, rcbest); + if (sps == true) { + bool magfield = Magfield2Central(false, g, jbest, z, r, Bz, Br, rc); + if (magfield == true) + goto labelend; + } + // --------------------------- + // Step 4. We try again the remote series expansion + // (with remote source point at the group center) + if (ro2 > rorem2 * 1.0204) { + bool magfield = Magfield2Remote(false, g, z, r, Bz, Br, rc); + if (magfield == true) + goto labelend; + } + // --------------------------- + // Step 5. We try again the central zonal method, by searching now all central source points. + sps = SourcePointSearching(false, g, z, r, 0, jbest, rcbest); + if (sps == true) { + bool magfield = Magfield2Central(false, g, jbest, z, r, Bz, Br, rc); + if (magfield == true) + goto labelend; + } + // ----------------------------- + // Step 6. Unfortunately, no appropriate central or remote expansion + // was found. We have to use coil calculations: + { + B[0] = B[1] = B[2] = 0.; + double Bcoil[3]; + bool magfieldcoil; + for (int c = 0; c < Nc[g]; c++) { + int i = Cin[g][c]; + magfieldcoil = MagfieldCoil(i, P, Bcoil); + for (int k = 0; k < 3; k++) + B[k] += Bcoil[k]; + if (magfieldcoil == false) + magfieldgroup = false; + } + return magfieldgroup; + } // -// We start now the magnetic field calculation. -// First we define the local coordinate system of the group, and the -// z and r local coordinates of the field point P. -// Local coordinate system: -// origo at Z[g][0], local z axis parallel to group axis, in u vector direction. - bool magfieldgroup=true; - double C[3], u[3], w[3]; - for(int k=0; k<3; k++) - C[k]=Line[g][k]; // group origo - for(int k=0; k<3; k++) - u[k]=Line[g][3+k]; // group direction unit vector -// Local cylindrical z and r coordinates of the field point P (relative to group origo and direction): - double Ploc[3], Pr[3]; - for(int k=0; k<3; k++) - Ploc[k]=P[k]-C[k]; - double z=Ploc[1]*u[1]+Ploc[2]*u[2]+Ploc[0]*u[0]; - for(int k=0; k<3; k++) - Pr[k]=Ploc[k]-z*u[k]; - double r=sqrt(Pr[1]*Pr[1]+Pr[2]*Pr[2]+Pr[0]*Pr[0]); - double r2=r*r; -// -// Starting now Bz and Br calculation. - bool sps; - int jbest; double rcbest; -// --------------------------- -// Step 1. We try first the remote series expansion -// (with remote source point at the group center) - double ro2=r2+(z-z0remG[g])*(z-z0remG[g]); - double rorem2=roremG[g]*roremG[g]; - double Bz, Br, rc; - if(ro2>rorem2*2.) - { - bool magfield=Magfield2Remote(false, g, z, r, Bz, Br, rc); - if(magfield==true) - goto labelend; - } -// --------------------------- -// Step 2. We try the central zonal method in the beginning of tracking, -// by searching all central source points. - if(jlastG[g]==-2) - { - jlastG[g]=-1; - sps=SourcePointSearching(false, g, z, r, 0, jbest, rcbest); - if(sps==true) - { - bool magfield=Magfield2Central(false, g, jbest, z, r, Bz, Br, rc); - if(magfield==true) - goto labelend; - } - } -// --------------------------- -// Step 3. The field point is close to the group. We try next the central -// Legendre polynomial expansion. If some source point for this group -// has been already used (jlastG[g]>-1), we search the central source -// point j close to the old source point jlastG[g]. - sps=false; - if(jlastG[g]>-1) - sps=SourcePointSearching(false, g, z, r, 1, jbest, rcbest); - if(jlastG[g]>-1 && sps==false) - sps=SourcePointSearching(false, g, z, r, 2, jbest, rcbest); - if(sps==true) - { - bool magfield=Magfield2Central(false, g, jbest, z, r, Bz, Br, rc); - if(magfield==true) - goto labelend; - } -// --------------------------- -// Step 4. We try again the remote series expansion -// (with remote source point at the group center) - if(ro2>rorem2*1.0204) - { - bool magfield=Magfield2Remote(false, g, z, r, Bz, Br, rc); - if(magfield==true) - goto labelend; - } -// --------------------------- -// Step 5. We try again the central zonal method, by searching now all central source points. - sps=SourcePointSearching(false, g, z, r, 0, jbest, rcbest); - if(sps==true) - { - bool magfield=Magfield2Central(false, g, jbest, z, r, Bz, Br, rc); - if(magfield==true) - goto labelend; - } -// ----------------------------- -// Step 6. Unfortunately, no appropriate central or remote expansion -// was found. We have to use coil calculations: - { - B[0]=B[1]=B[2]=0.; - double Bcoil[3]; - bool magfieldcoil; - for(int c=0; c-1) - { rcbest=sqrt(rcmin2); jlast[i]=jbest; return true; } - else - return false; - } // end of coil=true, type=0 case - else // only the 3 or the 11 source points near the last source point jlast are searched - { - if(jlast[i]<0) - return false; - int delj; - if(type==1) - delj=1; - else - delj=5; - int jmin=jlast[i]-delj; if(jmin<0) jmin=0; - int jmax=jlast[i]+delj; if(jmax>Nsp[i]-1) jmax=Nsp[i]-1; - rcmin2=rclimit*rclimit; jbest=-1; - for(int j=jmin; j<=jmax; j++) - { - delz=z-z0cen[i][j]; ro2=r2+delz*delz; rocen2=rocen[i][j]*rocen[i][j]; - if(ro2-1) - { rcbest=sqrt(rcmin2); jlast[i]=jbest; return true; } - else - return false; - } // end of coil=true, type=1 or 2 case - } // end of coil case - else // group case - { - int g=ig; - if(type==0) // all central source points of the group are searched - { - rcmin2=rclimit*rclimit; jbest=-1; - for(int j=0; j-1) - { rcbest=sqrt(rcmin2); jlastG[g]=jbest; return true; } - else - return false; - } // end of coil=false, type=0 case - else // only the 3 or the 11 source points near the last source point jlast are searched - { - if(jlastG[g]<0) - return false; - int delj; - if(type==1) - delj=1; - else - delj=5; - int jmin=jlastG[g]-delj; if(jmin<0) jmin=0; - int jmax=jlastG[g]+delj; if(jmax>NspG[g]-1) jmax=NspG[g]-1; - rcmin2=rclimit*rclimit; jbest=-1; - for(int j=jmin; j<=jmax; j++) - { - delz=z-z0cenG[g][j]; ro2=r2+delz*delz; rocen2=rocenG[g][j]*rocenG[g][j]; - if(ro2-1) - { rcbest=sqrt(rcmin2); jlastG[g]=jbest; return true; } - else - return false; - } // end of coil=false, type=1 or 2 case - } // end of group case + // This function searches the optimal central source point either + // for a coil (if coil==true) or for a group (if coil==false). + // Coil or group index: ig. + // Field point is defined by z and r (relative to the coil center or in the + // group Z-system, and to the coil or group direction). + // type=0: all central source points of the coil or group are searched; + // type=1: only the 3 source points near the last source point jlast are searched; + // type=2: only the 11 source points near the last source point jlast are searched. + // If a good central source point is found: the function return value is true; + // if no good central source point is found: the function return value is false. + // The best source point index found by this function: jbest (should be used only with true return value); + // the corresponding best convergence ratio (for source point jbest): rcbest (only for true return value). + double rcmin2, delz, ro2, rocen2; + double r2 = r * r; + if (coil == true) // coil case + { + int i = ig; + if (type == 0) // all central source points of the coil are searched + { + rcmin2 = rclimit * rclimit; + jbest = -1; + for (int j = 0; j < Nsp[i]; j++) { + delz = z - z0cen[i][j]; + ro2 = r2 + delz * delz; + rocen2 = rocen[i][j] * rocen[i][j]; + if (ro2 < rocen2 * rcmin2) { + rcmin2 = ro2 / rocen2; + jbest = j; + } + } + if (jbest > -1) { + rcbest = sqrt(rcmin2); + jlast[i] = jbest; + return true; + } + else + return false; + } // end of coil=true, type=0 case + else // only the 3 or the 11 source points near the last source point jlast are searched + { + if (jlast[i] < 0) + return false; + int delj; + if (type == 1) + delj = 1; + else + delj = 5; + int jmin = jlast[i] - delj; + if (jmin < 0) + jmin = 0; + int jmax = jlast[i] + delj; + if (jmax > Nsp[i] - 1) + jmax = Nsp[i] - 1; + rcmin2 = rclimit * rclimit; + jbest = -1; + for (int j = jmin; j <= jmax; j++) { + delz = z - z0cen[i][j]; + ro2 = r2 + delz * delz; + rocen2 = rocen[i][j] * rocen[i][j]; + if (ro2 < rocen2 * rcmin2) { + rcmin2 = ro2 / rocen2; + jbest = j; + } + } + if (jbest > -1) { + rcbest = sqrt(rcmin2); + jlast[i] = jbest; + return true; + } + else + return false; + } // end of coil=true, type=1 or 2 case + } // end of coil case + else // group case + { + int g = ig; + if (type == 0) // all central source points of the group are searched + { + rcmin2 = rclimit * rclimit; + jbest = -1; + for (int j = 0; j < NspG[g]; j++) { + delz = z - z0cenG[g][j]; + ro2 = r2 + delz * delz; + rocen2 = rocenG[g][j] * rocenG[g][j]; + if (ro2 < rocen2 * rcmin2) { + rcmin2 = ro2 / rocen2; + jbest = j; + } + } + if (jbest > -1) { + rcbest = sqrt(rcmin2); + jlastG[g] = jbest; + return true; + } + else + return false; + } // end of coil=false, type=0 case + else // only the 3 or the 11 source points near the last source point jlast are searched + { + if (jlastG[g] < 0) + return false; + int delj; + if (type == 1) + delj = 1; + else + delj = 5; + int jmin = jlastG[g] - delj; + if (jmin < 0) + jmin = 0; + int jmax = jlastG[g] + delj; + if (jmax > NspG[g] - 1) + jmax = NspG[g] - 1; + rcmin2 = rclimit * rclimit; + jbest = -1; + for (int j = jmin; j <= jmax; j++) { + delz = z - z0cenG[g][j]; + ro2 = r2 + delz * delz; + rocen2 = rocenG[g][j] * rocenG[g][j]; + if (ro2 < rocen2 * rcmin2) { + rcmin2 = ro2 / rocen2; + jbest = j; + } + } + if (jbest > -1) { + rcbest = sqrt(rcmin2); + jlastG[g] = jbest; + return true; + } + else + return false; + } // end of coil=false, type=1 or 2 case + } // end of group case } - - - - // #include "CarlsonEllipticIntegrals.cc" ///////////////////////////////////////////////////// @@ -2264,167 +2287,194 @@ bool MagfieldCoils::SourcePointSearching(bool coil, int ig, double z, double r, ///////////////////////////////////////////////////// -double MagfieldCoils::RF_Carlson(double x,double y,double z) +double MagfieldCoils::RF_Carlson(double x, double y, double z) { -// This function computes Carlson's elliptic integral of the first kind: -// R_F(x,y,z). x, y, z must be nonnegative, and at most one can be zero -// (see: Press et al., Numerical Recipes, Sec. 6.11). - const double ERRTOL=0.002,TINY=1.e-38,BIG=1.e38,C1=1./24.,C2=0.1, - C3=3./44.,C4=1./14.,THIRD=1./3.; - double alamb,ave,delx,dely,delz,e2,e3,sqrtx,sqrty,sqrtz,xt,yt,zt; - if(FMIN3(x,y,z)<0. || FMIN3(x+y,x+z,y+z)BIG) - { - puts("Message from function MagfieldCoils::RF_Carlson: invalid arguments !!!"); - puts("Program running is stopped !!! "); - exit(1); - } - xt=x; yt=y; zt=z; - do - { - sqrtx=sqrt(xt); sqrty=sqrt(yt); sqrtz=sqrt(zt); - alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz; - xt=0.25*(xt+alamb); yt=0.25*(yt+alamb); zt=0.25*(zt+alamb); - ave=THIRD*(xt+yt+zt); - delx=(ave-xt)/ave; dely=(ave-yt)/ave; delz=(ave-zt)/ave; - } - while (FMAX3(fabs(delx),fabs(dely),fabs(delz))>ERRTOL); - e2=delx*dely-delz*delz; - e3=delx*dely*delz; - return (1.+(C1*e2-C2-C3*e3)*e2+C4*e3)/sqrt(ave); + // This function computes Carlson's elliptic integral of the first kind: + // R_F(x,y,z). x, y, z must be nonnegative, and at most one can be zero + // (see: Press et al., Numerical Recipes, Sec. 6.11). + const double ERRTOL = 0.002, TINY = 1.e-38, BIG = 1.e38, C1 = 1. / 24., C2 = 0.1, C3 = 3. / 44., C4 = 1. / 14., + THIRD = 1. / 3.; + double alamb, ave, delx, dely, delz, e2, e3, sqrtx, sqrty, sqrtz, xt, yt, zt; + if (FMIN3(x, y, z) < 0. || FMIN3(x + y, x + z, y + z) < TINY || FMAX3(x, y, z) > BIG) { + puts("Message from function MagfieldCoils::RF_Carlson: invalid arguments !!!"); + puts("Program running is stopped !!! "); + exit(1); + } + xt = x; + yt = y; + zt = z; + do { + sqrtx = sqrt(xt); + sqrty = sqrt(yt); + sqrtz = sqrt(zt); + alamb = sqrtx * (sqrty + sqrtz) + sqrty * sqrtz; + xt = 0.25 * (xt + alamb); + yt = 0.25 * (yt + alamb); + zt = 0.25 * (zt + alamb); + ave = THIRD * (xt + yt + zt); + delx = (ave - xt) / ave; + dely = (ave - yt) / ave; + delz = (ave - zt) / ave; + } while (FMAX3(fabs(delx), fabs(dely), fabs(delz)) > ERRTOL); + e2 = delx * dely - delz * delz; + e3 = delx * dely * delz; + return (1. + (C1 * e2 - C2 - C3 * e3) * e2 + C4 * e3) / sqrt(ave); } /////////////////////////////////////////////////////////////// -double MagfieldCoils::RD_Carlson(double x,double y,double z) +double MagfieldCoils::RD_Carlson(double x, double y, double z) { -// This function computes Carlson's elliptic integral of the second kind: -// R_D(x,y,z). x and y must be nonnegative, and at most one can be zero. -// z must be positive -// (see: Press et al., Numerical Recipes, Sec. 6.11). - const double ERRTOL=0.0015,TINY=1.e-25,BIG=1.e22,C1=3./14.,C2=1./6., - C3=9./22.,C4=3./26.,C5=0.25*C3,C6=1.5*C4; - double alamb,ave,delx,dely,delz,ea,eb,ec,ed,ee,fac,sum, - sqrtx,sqrty,sqrtz,xt,yt,zt; - if(FMIN(x,y)<0. || FMIN(x+y,z)BIG) - { - puts("Message from function MagfieldCoils::RD_Carlson: invalid arguments !!!"); - puts("Program running is stopped !!! "); - exit(1); - } - xt=x; yt=y; zt=z; - sum=0.; fac=1.; - do - { - sqrtx=sqrt(xt); sqrty=sqrt(yt); sqrtz=sqrt(zt); - alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz; - sum+=fac/(sqrtz*(zt+alamb)); - fac=0.25*fac; - xt=0.25*(xt+alamb); yt=0.25*(yt+alamb); zt=0.25*(zt+alamb); - ave=0.2*(xt+yt+3.*zt); - delx=(ave-xt)/ave; dely=(ave-yt)/ave; delz=(ave-zt)/ave; - } - while (FMAX3(fabs(delx),fabs(dely),fabs(delz))>ERRTOL); - ea=delx*dely; eb=delz*delz; - ec=ea-eb; ed=ea-6.*eb; ee=ed+ec+ec; - return 3.*sum+fac*(1.+ed*(-C1+C5*ed-C6*delz*ee)+ - delz*(C2*ee+delz*(-C3*ec+delz*C4*ea)))/(ave*sqrt(ave)); + // This function computes Carlson's elliptic integral of the second kind: + // R_D(x,y,z). x and y must be nonnegative, and at most one can be zero. + // z must be positive + // (see: Press et al., Numerical Recipes, Sec. 6.11). + const double ERRTOL = 0.0015, TINY = 1.e-25, BIG = 1.e22, C1 = 3. / 14., C2 = 1. / 6., C3 = 9. / 22., C4 = 3. / 26., + C5 = 0.25 * C3, C6 = 1.5 * C4; + double alamb, ave, delx, dely, delz, ea, eb, ec, ed, ee, fac, sum, sqrtx, sqrty, sqrtz, xt, yt, zt; + if (FMIN(x, y) < 0. || FMIN(x + y, z) < TINY || FMAX3(x, y, z) > BIG) { + puts("Message from function MagfieldCoils::RD_Carlson: invalid arguments !!!"); + puts("Program running is stopped !!! "); + exit(1); + } + xt = x; + yt = y; + zt = z; + sum = 0.; + fac = 1.; + do { + sqrtx = sqrt(xt); + sqrty = sqrt(yt); + sqrtz = sqrt(zt); + alamb = sqrtx * (sqrty + sqrtz) + sqrty * sqrtz; + sum += fac / (sqrtz * (zt + alamb)); + fac = 0.25 * fac; + xt = 0.25 * (xt + alamb); + yt = 0.25 * (yt + alamb); + zt = 0.25 * (zt + alamb); + ave = 0.2 * (xt + yt + 3. * zt); + delx = (ave - xt) / ave; + dely = (ave - yt) / ave; + delz = (ave - zt) / ave; + } while (FMAX3(fabs(delx), fabs(dely), fabs(delz)) > ERRTOL); + ea = delx * dely; + eb = delz * delz; + ec = ea - eb; + ed = ea - 6. * eb; + ee = ed + ec + ec; + return 3. * sum + + fac * (1. + ed * (-C1 + C5 * ed - C6 * delz * ee) + delz * (C2 * ee + delz * (-C3 * ec + delz * C4 * ea))) / + (ave * sqrt(ave)); } /////////////////////////////////////////////////////////////// -double MagfieldCoils::RJ_Carlson(double x,double y,double z,double p) +double MagfieldCoils::RJ_Carlson(double x, double y, double z, double p) { -// This function computes Carlson's elliptic integral of the third kind: -// R_J(x,y,z,p). x, y and z must be nonnegative, and at most one can be zero. -// p must be nonzero. If p<0, the Cauchy principal value is returned. -// (see: Press et al., Numerical Recipes, Sec. 6.11). - const double ERRTOL=0.0015,TINY=1.e-20,BIG=1.e12,C1=3./14.,C2=1./3., - C3=3./22.,C4=3./26.,C5=0.75*C3,C6=1.5*C4,C7=0.5*C2,C8=2.*C3; - double alamb,alpha,ans,ave,beta,delp,delx,dely,delz,ea,eb,ec,ed,ee, - fac,pt,rho,sum,sqrtx,sqrty,sqrtz,tau,xt,yt,zt; - double rcx = 0.; double b = 0.; double a = 0.; - if(FMIN3(x,y,z)<0. || FMIN(FMIN(x+y,x+z),FMIN(y+z,fabs(p)))BIG) - { - puts("Message from function MagfieldCoils::RJ_Carlson: invalid arguments !!!"); - puts("Program running is stopped !!! "); - exit(1); - } - sum=0.; fac=1.; - if(p>0.) - { xt=x; yt=y; zt=z; pt=p; } - else - { - xt=FMIN3(x,y,z); - zt=FMAX3(x,y,z); - yt=x+y+z-xt-zt; - a=1./(yt-p); - b=a*(zt-yt)*(yt-xt); - pt=yt+b; - rho=xt*zt/yt; - tau=p*pt/yt; - rcx=RC_Carlson(rho,tau); - } - do - { - sqrtx=sqrt(xt); sqrty=sqrt(yt); sqrtz=sqrt(zt); - alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz; - alpha=pow2(pt*(sqrtx+sqrty+sqrtz)+sqrtx*sqrty*sqrtz); - beta=pt*pow2(pt+alamb); - sum+=fac*RC_Carlson(alpha,beta);; - fac=0.25*fac; - xt=0.25*(xt+alamb); yt=0.25*(yt+alamb); zt=0.25*(zt+alamb); - pt=0.25*(pt+alamb); - ave=0.2*(xt+yt+zt+2.*pt); - delx=(ave-xt)/ave; dely=(ave-yt)/ave; delz=(ave-zt)/ave; - delp=(ave-pt)/ave; - } - while (FMAX(FMAX(fabs(delx),fabs(dely)),FMAX(fabs(delz),fabs(delp))) - >ERRTOL); - ea=delx*(dely+delz)+dely*delz; - eb=delx*dely*delz; - ec=delp*delp; - ed=ea-3.*ec; - ee=eb+2.*delp*(ea-ec); - ans=3.*sum+fac*(1.+ed*(-C1+C5*ed-C6*ee)+eb*(C7+delp*(-C8+delp*C4))+ - delp*ea*(C2-delp*C3)-C2*delp*ec)/(ave*sqrt(ave)); - if(p<0.) - ans=a*(b*ans+3.*(rcx-RF_Carlson(xt,yt,zt))); - return ans; + // This function computes Carlson's elliptic integral of the third kind: + // R_J(x,y,z,p). x, y and z must be nonnegative, and at most one can be zero. + // p must be nonzero. If p<0, the Cauchy principal value is returned. + // (see: Press et al., Numerical Recipes, Sec. 6.11). + const double ERRTOL = 0.0015, TINY = 1.e-20, BIG = 1.e12, C1 = 3. / 14., C2 = 1. / 3., C3 = 3. / 22., C4 = 3. / 26., + C5 = 0.75 * C3, C6 = 1.5 * C4, C7 = 0.5 * C2, C8 = 2. * C3; + double alamb, alpha, ans, ave, beta, delp, delx, dely, delz, ea, eb, ec, ed, ee, fac, pt, rho, sum, sqrtx, sqrty, + sqrtz, tau, xt, yt, zt; + double rcx = 0.; + double b = 0.; + double a = 0.; + if (FMIN3(x, y, z) < 0. || FMIN(FMIN(x + y, x + z), FMIN(y + z, fabs(p))) < TINY || + FMAX(FMAX(x, y), FMAX(z, fabs(p))) > BIG) { + puts("Message from function MagfieldCoils::RJ_Carlson: invalid arguments !!!"); + puts("Program running is stopped !!! "); + exit(1); + } + sum = 0.; + fac = 1.; + if (p > 0.) { + xt = x; + yt = y; + zt = z; + pt = p; + } + else { + xt = FMIN3(x, y, z); + zt = FMAX3(x, y, z); + yt = x + y + z - xt - zt; + a = 1. / (yt - p); + b = a * (zt - yt) * (yt - xt); + pt = yt + b; + rho = xt * zt / yt; + tau = p * pt / yt; + rcx = RC_Carlson(rho, tau); + } + do { + sqrtx = sqrt(xt); + sqrty = sqrt(yt); + sqrtz = sqrt(zt); + alamb = sqrtx * (sqrty + sqrtz) + sqrty * sqrtz; + alpha = pow2(pt * (sqrtx + sqrty + sqrtz) + sqrtx * sqrty * sqrtz); + beta = pt * pow2(pt + alamb); + sum += fac * RC_Carlson(alpha, beta); + ; + fac = 0.25 * fac; + xt = 0.25 * (xt + alamb); + yt = 0.25 * (yt + alamb); + zt = 0.25 * (zt + alamb); + pt = 0.25 * (pt + alamb); + ave = 0.2 * (xt + yt + zt + 2. * pt); + delx = (ave - xt) / ave; + dely = (ave - yt) / ave; + delz = (ave - zt) / ave; + delp = (ave - pt) / ave; + } while (FMAX(FMAX(fabs(delx), fabs(dely)), FMAX(fabs(delz), fabs(delp))) > ERRTOL); + ea = delx * (dely + delz) + dely * delz; + eb = delx * dely * delz; + ec = delp * delp; + ed = ea - 3. * ec; + ee = eb + 2. * delp * (ea - ec); + ans = 3. * sum + fac * + (1. + ed * (-C1 + C5 * ed - C6 * ee) + eb * (C7 + delp * (-C8 + delp * C4)) + + delp * ea * (C2 - delp * C3) - C2 * delp * ec) / + (ave * sqrt(ave)); + if (p < 0.) + ans = a * (b * ans + 3. * (rcx - RF_Carlson(xt, yt, zt))); + return ans; } /////////////////////////////////////////////////////////////// -double MagfieldCoils::RC_Carlson(double x,double y) +double MagfieldCoils::RC_Carlson(double x, double y) { -// This function computes Carlson's degenerate elliptic integral: -// R_C(x,y). x must be nonnegative, and y must be nonzero. -// If y<0, the Cauchy principal value is returned. -// (see: Press et al., Numerical Recipes, Sec. 6.11). - const double ERRTOL=0.001,TINY=1.e-38,BIG=1.e38,SQRTNY=1.e-19,TNBG=TINY*BIG, - COMP1=2.236/SQRTNY,COMP2=TNBG*TNBG/25.,THIRD=1./3., - C1=0.3,C2=1./7.,C3=0.375,C4=9./22.; - double alamb,ave,s,w,xt,yt; - if(x<0. || y==0. || (x+fabs(y))BIG || - (y<-COMP1 && x>0. && x0.) - { xt=x; yt=y; w=1.;} - else - { xt=x-y; yt=-y; w=sqrt(x)/sqrt(xt); } - do - { - alamb=2.*sqrt(xt)*sqrt(yt)+yt; - xt=0.25*(xt+alamb); yt=0.25*(yt+alamb); - ave=THIRD*(xt+yt+yt); - s=(yt-ave)/ave; - } - while (fabs(s)>ERRTOL); - return w*(1.+s*s*(C1+s*(C2+s*(C3+s*C4))))/sqrt(ave); + // This function computes Carlson's degenerate elliptic integral: + // R_C(x,y). x must be nonnegative, and y must be nonzero. + // If y<0, the Cauchy principal value is returned. + // (see: Press et al., Numerical Recipes, Sec. 6.11). + const double ERRTOL = 0.001, TINY = 1.e-38, BIG = 1.e38, SQRTNY = 1.e-19, TNBG = TINY * BIG, COMP1 = 2.236 / SQRTNY, + COMP2 = TNBG * TNBG / 25., THIRD = 1. / 3., C1 = 0.3, C2 = 1. / 7., C3 = 0.375, C4 = 9. / 22.; + double alamb, ave, s, w, xt, yt; + if (x < 0. || y == 0. || (x + fabs(y)) < TINY || x + fabs(y) > BIG || (y < -COMP1 && x > 0. && x < COMP2)) { + puts("Message from function MagfieldCoils::RC_Carlson: invalid arguments !!!"); + puts("Program running is stopped !!! "); + exit(1); + } + if (y > 0.) { + xt = x; + yt = y; + w = 1.; + } + else { + xt = x - y; + yt = -y; + w = sqrt(x) / sqrt(xt); + } + do { + alamb = 2. * sqrt(xt) * sqrt(yt) + yt; + xt = 0.25 * (xt + alamb); + yt = 0.25 * (yt + alamb); + ave = THIRD * (xt + yt + yt); + s = (yt - ave) / ave; + } while (fabs(s) > ERRTOL); + return w * (1. + s * s * (C1 + s * (C2 + s * (C3 + s * C4)))) / sqrt(ave); } // #include "GaussLegendreIntegration.cc" @@ -2439,112 +2489,162 @@ void MagfieldCoils::GaussLegendreIntegration(int& N, double a, double b, double* // This code can change the input integer N! // Dimension of arrays x and w has to be minimum N! { - const double x2[1]={0.5773502691896257}; - const double w2[1]={1.}; - const double x3[2]={0.,0.7745966692414834}; - const double w3[2]={0.8888888888888888,0.5555555555555556}; - const double x4[2]={0.3399810435848563,0.8611363115940526}; - const double w4[2]={0.6521451548625461,0.3478548451374538}; - const double x6[3]={0.2386191860831969,0.6612093864662645,0.9324695142031521}; - const double w6[3]={0.4679139345726910,0.3607615730481386,0.1713244923791704}; - const double x8[4]={0.1834346424956498,0.5255324099163290,0.7966664774136267,0.9602898564975363 }; - const double w8[4]={0.3626837833783620,0.3137066458778873,0.2223810344533745,0.1012285362903763}; - const double x16[8]={0.09501250983763744,0.28160355077925891,0.45801677765722739,0.61787624440264375, - 0.75540440835500303,0.86563120238783174,0.94457502307323258 ,0.98940093499164993}; - const double w16[8]={0.189450610455068496,0.182603415044923589,0.169156519395002532,0.149595988816576731, - 0.124628971255533872,0.095158511682492785,0.062253523938647892,0.027152459411754095}; - const double x32[16]={0.048307665687738316,0.144471961582796493,0.239287362252137075 ,0.331868602282127650, - 0.421351276130635345,0.506899908932229390,0.587715757240762329,0.663044266930215201, - 0.732182118740289680,0.794483795967942407,0.849367613732569970,0.896321155766052124, - 0.934906075937739689,0.964762255587506431,0.985611511545268335,0.997263861849481564}; - const double w32[16]={0.09654008851472780056,0.09563872007927485942,0.09384439908080456564,0.09117387869576388471, - 0.08765209300440381114,0.08331192422694675522,0.07819389578707030647,0.07234579410884850625, - 0.06582222277636184684,0.05868409347853554714,0.05099805926237617619,0.04283589802222680057, - 0.03427386291302143313,0.02539206530926205956,0.01627439473090567065,0.00701861000947009660}; - static double X[33][32]; - static double W[33][32]; - static bool start=true; - int m, n; - if(start==true) - { - m=2; - X[m][0]=-x2[0]; X[m][1]=x2[0]; - W[m][0]=w2[0]; W[m][1]=w2[0]; - m=3; - X[m][0]=-x3[1]; X[m][1]=x3[0]; X[m][2]=x3[1]; - W[m][0]=w3[1]; W[m][1]=w3[0]; W[m][2]=w3[1]; - m=4; n=m/2; - for(int i=0; i=5 && N<=6) - N=6; - else if(N>=7 && N<=11) - N=8; - else if(N>=12 && N<=20) - N=16; - else if(N>=21 && N<=45) - N=32; - else if(N>=46) - { - int M=N/32+1; N=32*M; + const double x2[1] = {0.5773502691896257}; + const double w2[1] = {1.}; + const double x3[2] = {0., 0.7745966692414834}; + const double w3[2] = {0.8888888888888888, 0.5555555555555556}; + const double x4[2] = {0.3399810435848563, 0.8611363115940526}; + const double w4[2] = {0.6521451548625461, 0.3478548451374538}; + const double x6[3] = {0.2386191860831969, 0.6612093864662645, 0.9324695142031521}; + const double w6[3] = {0.4679139345726910, 0.3607615730481386, 0.1713244923791704}; + const double x8[4] = {0.1834346424956498, 0.5255324099163290, 0.7966664774136267, 0.9602898564975363}; + const double w8[4] = {0.3626837833783620, 0.3137066458778873, 0.2223810344533745, 0.1012285362903763}; + const double x16[8] = {0.09501250983763744, + 0.28160355077925891, + 0.45801677765722739, + 0.61787624440264375, + 0.75540440835500303, + 0.86563120238783174, + 0.94457502307323258, + 0.98940093499164993}; + const double w16[8] = {0.189450610455068496, + 0.182603415044923589, + 0.169156519395002532, + 0.149595988816576731, + 0.124628971255533872, + 0.095158511682492785, + 0.062253523938647892, + 0.027152459411754095}; + const double x32[16] = {0.048307665687738316, + 0.144471961582796493, + 0.239287362252137075, + 0.331868602282127650, + 0.421351276130635345, + 0.506899908932229390, + 0.587715757240762329, + 0.663044266930215201, + 0.732182118740289680, + 0.794483795967942407, + 0.849367613732569970, + 0.896321155766052124, + 0.934906075937739689, + 0.964762255587506431, + 0.985611511545268335, + 0.997263861849481564}; + const double w32[16] = {0.09654008851472780056, + 0.09563872007927485942, + 0.09384439908080456564, + 0.09117387869576388471, + 0.08765209300440381114, + 0.08331192422694675522, + 0.07819389578707030647, + 0.07234579410884850625, + 0.06582222277636184684, + 0.05868409347853554714, + 0.05099805926237617619, + 0.04283589802222680057, + 0.03427386291302143313, + 0.02539206530926205956, + 0.01627439473090567065, + 0.00701861000947009660}; + static double X[33][32]; + static double W[33][32]; + static bool start = true; + int m, n; + if (start == true) { + m = 2; + X[m][0] = -x2[0]; + X[m][1] = x2[0]; + W[m][0] = w2[0]; + W[m][1] = w2[0]; + m = 3; + X[m][0] = -x3[1]; + X[m][1] = x3[0]; + X[m][2] = x3[1]; + W[m][0] = w3[1]; + W[m][1] = w3[0]; + W[m][2] = w3[1]; + m = 4; + n = m / 2; + for (int i = 0; i < n; i++) { + X[m][i] = -x4[n - 1 - i]; + W[m][i] = w4[n - 1 - i]; + X[m][n + i] = x4[i]; + W[m][n + i] = w4[i]; + } + m = 6; + n = m / 2; + for (int i = 0; i < n; i++) { + X[m][i] = -x6[n - 1 - i]; + W[m][i] = w6[n - 1 - i]; + X[m][n + i] = x6[i]; + W[m][n + i] = w6[i]; + } + m = 8; + n = m / 2; + for (int i = 0; i < n; i++) { + X[m][i] = -x8[n - 1 - i]; + W[m][i] = w8[n - 1 - i]; + X[m][n + i] = x8[i]; + W[m][n + i] = w8[i]; + } + m = 16; + n = m / 2; + for (int i = 0; i < n; i++) { + X[m][i] = -x16[n - 1 - i]; + W[m][i] = w16[n - 1 - i]; + X[m][n + i] = x16[i]; + W[m][n + i] = w16[i]; + } + m = 32; + n = m / 2; + for (int i = 0; i < n; i++) { + X[m][i] = -x32[n - 1 - i]; + W[m][i] = w32[n - 1 - i]; + X[m][n + i] = x32[i]; + W[m][n + i] = w32[i]; + } + start = false; } -// - double Cp=(b+a)/2.; double Cm=(b-a)/2.; - if(N<=32) - { - for(int i=0; i= 5 && N <= 6) + N = 6; + else if (N >= 7 && N <= 11) + N = 8; + else if (N >= 12 && N <= 20) + N = 16; + else if (N >= 21 && N <= 45) + N = 32; + else if (N >= 46) { + int M = N / 32 + 1; + N = 32 * M; } - else - { - int M=N/32; // number of subintervals - double d=(b-a)/M; // length of subinterval - for(int m=0; m Date: Fri, 1 Dec 2023 14:41:08 +0100 Subject: [PATCH 06/15] Remove obsolete compiler option for MagFieldCoils --- KEMField/Source/ExternalFields/MagfieldCoils/CMakeLists.txt | 2 -- 1 file changed, 2 deletions(-) diff --git a/KEMField/Source/ExternalFields/MagfieldCoils/CMakeLists.txt b/KEMField/Source/ExternalFields/MagfieldCoils/CMakeLists.txt index bfd5bdb4b..a492cd90e 100644 --- a/KEMField/Source/ExternalFields/MagfieldCoils/CMakeLists.txt +++ b/KEMField/Source/ExternalFields/MagfieldCoils/CMakeLists.txt @@ -39,8 +39,6 @@ foreach( BASENAME ${MAGFIELDCOILS_SOURCE_BASENAMES} ) list( APPEND MAGFIELDCOILS_SOURCEFILES ${MAGFIELDCOILS_SOURCE_PATH}/${BASENAME} ) endforeach( BASENAME ) -set_property(SOURCE ${MAGFIELDCOILS_SOURCEFILES} APPEND_STRING PROPERTY COMPILE_OPTIONS -Wno-error) # FIXME - add_library (KEMMagfieldCoils SHARED ${MAGFIELDCOILS_SOURCEFILES} ${MAGFIELDCOILS_HEADERFILES}) target_include_directories(KEMMagfieldCoils From 07553ad1eb216297eb0737fae8f9ca739052a2f4 Mon Sep 17 00:00:00 2001 From: Richard Salomon Date: Fri, 8 Dec 2023 16:58:04 +0100 Subject: [PATCH 07/15] Refactoring MagfieldCoils member variables to be in line with general naming scheme --- .../MagfieldCoils/include/MagfieldCoils.h | 78 +- .../MagfieldCoils/src/MagfieldCoils.cxx | 1270 ++++++++--------- 2 files changed, 674 insertions(+), 674 deletions(-) diff --git a/KEMField/Source/ExternalFields/MagfieldCoils/include/MagfieldCoils.h b/KEMField/Source/ExternalFields/MagfieldCoils/include/MagfieldCoils.h index 70b95eac6..70d75ec43 100755 --- a/KEMField/Source/ExternalFields/MagfieldCoils/include/MagfieldCoils.h +++ b/KEMField/Source/ExternalFields/MagfieldCoils/include/MagfieldCoils.h @@ -82,51 +82,51 @@ class MagfieldCoils double Mag(double* vec); double Mag(double x, double y, double z); // VARIABLES: - int Ncoil; // number of coils; coil indexing: i=0, ..., Ncoil-1 - double** coil; // coil parameters: coil[i][j], i=0, ..., Ncoil-1; j=0,...,13. - string dirname; // name of directory where the coil and source coefficient data files are stored - string coilfilename; // name of coil input file - string objectname; // this string identifies the Magfield object (used as starting part of the source files) - int Nelliptic; // radial numerical integration parameter for magnetic field calc. with elliptic integrals - int nmax; // number of source coefficients for fixed z: nmax+1 - double epstol; // small tolerance parameter needed for the symmetry group definitions (distance: in meter) - int Ng; // number of coil symmetry groups; group indexing: g=0, ..., Ng-1 - int* G; // coil with global index i is in symmetry group g=G[i], i=0,...,Ncoil-1 - int* Nc; // number of coils in symmetry group g: Nc[g]; local coil indexing in group g: c=0,..., Nc[g]-1 - int** Cin; // global coil indices i in symmetry group g: i=Cin[g][c],...,c=0 to Nc[g]-1 (g=0,...,Ng-1) - double** Line; // symmetry group axis line of symmetry group g (g=0,...,Ng-1): + int fNcoil; // number of coils; coil indexing: i=0, ..., Ncoil-1 + double** fcoil; // coil parameters: coil[i][j], i=0, ..., Ncoil-1; j=0,...,13. + string fdirname; // name of directory where the coil and source coefficient data files are stored + string fcoilfilename; // name of coil input file + string fobjectname; // this string identifies the Magfield object (used as starting part of the source files) + int fNelliptic; // radial numerical integration parameter for magnetic field calc. with elliptic integrals + int fnmax; // number of source coefficients for fixed z: nmax+1 + double fepstol; // small tolerance parameter needed for the symmetry group definitions (distance: in meter) + int fNg; // number of coil symmetry groups; group indexing: g=0, ..., Ng-1 + int* fG; // coil with global index i is in symmetry group g=G[i], i=0,...,Ncoil-1 + int* fNc; // number of coils in symmetry group g: Nc[g]; local coil indexing in group g: c=0,..., Nc[g]-1 + int** fCin; // global coil indices i in symmetry group g: i=Cin[g][c],...,c=0 to Nc[g]-1 (g=0,...,Ng-1) + double** fLine; // symmetry group axis line of symmetry group g (g=0,...,Ng-1): // line point: Line[g][0], Line[g][1], Line[g][2]; line direction unit vector: Line[g][3], Line[g][4], Line[g][5]; - double** Z; // local z coordinate of coil center in symmetry group g: Z[g][c], c=0,..., Nc[g]-1; Z[g][0]=0 + double** fZ; // local z coordinate of coil center in symmetry group g: Z[g][c], c=0,..., Nc[g]-1; Z[g][0]=0 // (g: group index, c: local coil index) - double *C0, *C1, *c1, *c2, *c3, *c4, *c5, *c6, *c7, *c8, *c9, *c10, *c11, *c12; + double *fC0, *fC1, *fc1, *fc2, *fc3, *fc4, *fc5, *fc6, *fc7, *fc8, *fc9, *fc10, *fc11, *fc12; - double *Pp, *P; // Legendre polynomial P and first derivative Pp - double *Bzplus, *Brplus; // used in Bz and Br mag. field calc. + double *fPp, *fP; // Legendre polynomial P and first derivative Pp + double *fBzplus, *fBrplus; // used in Bz and Br mag. field calc. // Remote source constant variables: - double* Brem1; // 1-dim. remote source constants Brem1[n], n=2,...,nmax. - double* rorem; // remote convergence radius for coil i: rorem[i] - double** Brem; // remote source constant for coil i: Brem[i][n] (coil index i, source constant index n) - double* z0remG; // remote source point for group G: z0remG[g] - double* roremG; // remote convergence radius for group G: roremG[g] + double* fBrem1; // 1-dim. remote source constants Brem1[n], n=2,...,nmax. + double* frorem; // remote convergence radius for coil i: rorem[i] + double** fBrem; // remote source constant for coil i: Brem[i][n] (coil index i, source constant index n) + double* fz0remG; // remote source point for group G: z0remG[g] + double* froremG; // remote convergence radius for group G: roremG[g] double** - BremG; // remote source constant for symmetry group g: BremG[g][n] (group index g, source constant index n) + fBremG; // remote source constant for symmetry group g: BremG[g][n] (group index g, source constant index n) // Magnetic charge source constant variables: double** - Vrem; // magnetic charge remote source constant for coil i: Vrem[i][n] (coil index i, source const. index n) + fVrem; // magnetic charge remote source constant for coil i: Vrem[i][n] (coil index i, source const. index n) // Central source constant variables: - double* Bcen1; // 1-dim. central source constants Bcen1[n], n=0,...,nmax. - int* Nsp; // number of central source points for coil i: Nsp[i] (i=0,...,Ncoil-1) - double** z0cen; // central source point local z value for coil i: z0cen[i][j] (coil index i, source point index j) - double** rocen; // central convergence radius for coil i: rocen[i][j] (coil index i, source point index j) + double* fBcen1; // 1-dim. central source constants Bcen1[n], n=0,...,nmax. + int* fNsp; // number of central source points for coil i: Nsp[i] (i=0,...,Ncoil-1) + double** fz0cen; // central source point local z value for coil i: z0cen[i][j] (coil index i, source point index j) + double** frocen; // central convergence radius for coil i: rocen[i][j] (coil index i, source point index j) double*** - Bcen; // central source constant for coil i: Bcen[i][j][n] (coil index i, source point index j, source constant index n) - int* NspG; // number of central source points for group g: Nsp[g] (g=0,...,Ng-1) - double** z0cenG; // central source point local z value: z0cenG[g][j] (group index g, source point index j) - double** rocenG; // central convergence radius: rocenG[g][j] (group index g, source point index j) + fBcen; // central source constant for coil i: Bcen[i][j][n] (coil index i, source point index j, source constant index n) + int* fNspG; // number of central source points for group g: Nsp[g] (g=0,...,Ng-1) + double** fz0cenG; // central source point local z value: z0cenG[g][j] (group index g, source point index j) + double** frocenG; // central convergence radius: rocenG[g][j] (group index g, source point index j) double*** - BcenG; // central source constant: BcenG[g][j][n] (group index g, source point index j, source const. index n) - int *jlast, *jlastG; // last central source point index for coil and group calculation - double rclimit; + fBcenG; // central source constant: BcenG[g][j][n] (group index g, source point index j, source const. index n) + int *fjlast, *fjlastG; // last central source point index for coil and group calculation + double frclimit; }; @@ -183,10 +183,10 @@ inline double MagfieldCoils::Mag(double x, double y, double z) inline void MagfieldCoils::SetTrackingStart() { - for (int i = 0; i < Ncoil; i++) - jlast[i] = -1; - for (int g = 0; g < Ng; g++) - jlastG[g] = -1; + for (int i = 0; i < fNcoil; i++) + fjlast[i] = -1; + for (int g = 0; g < fNg; g++) + fjlastG[g] = -1; } diff --git a/KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx b/KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx index 069d946cd..304c76aea 100755 --- a/KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx +++ b/KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx @@ -8,19 +8,19 @@ MagfieldCoils::MagfieldCoils(string inputdirname, string inputobjectname, string // a, coil parameters are read from file coilfilename; coil parameters are tested. // b, symmetry groups are computed; coil and group parameters are written into file dirname+objectname+coil.txt // c, source constants are computed, and written into data files defined by dirname+objectname - dirname = inputdirname; // name of directory where the source coefficient files are stored - objectname = + fdirname = inputdirname; // name of directory where the source coefficient files are stored + fobjectname = inputobjectname; // this string identifies the Magfield object (used as starting part of the source files) - coilfilename = inputcoilfilename; // name of coil input file - Nelliptic = + fcoilfilename = inputcoilfilename; // name of coil input file + fNelliptic = inputNelliptic; // radial numerical integration parameter for magnetic field calc. with elliptic integrals - nmax = inputnmax; // number of source coefficients for fixed source point: nmax+1 - epstol = inputepstol; // small tolerance parameter needed for the symmetry group definitions + fnmax = inputnmax; // number of source coefficients for fixed source point: nmax+1 + fepstol = inputepstol; // small tolerance parameter needed for the symmetry group definitions // Recommended values: inputNelliptic=32, inputnmax=500, inputepstol=1.e-8. // CoilRead(); DynamicMemoryAllocations(); - rclimit = 0.99; + frclimit = 0.99; CoilGroupWrite(); // // The following functions compute all the source constants that are needed @@ -45,14 +45,14 @@ MagfieldCoils::MagfieldCoils(string inputdirname, string inputobjectname) // In this constructor: // a, coil and group parameters are read from file dirname+objectname+coil.txt // b, source constants are read from data files defined by dirname+objectname - dirname = inputdirname; // name of directory where the source coefficient files are stored - objectname = + fdirname = inputdirname; // name of directory where the source coefficient files are stored + fobjectname = inputobjectname; // this string identifies the Magfield object (used as starting part of the source files) // CoilGroupRead(); MagsourceRead(); DynamicMemoryAllocations(); - rclimit = 0.99; + frclimit = 0.99; } @@ -62,79 +62,79 @@ MagfieldCoils::~MagfieldCoils() { // One-dim arrays: - delete[] C0; - delete[] C1; - delete[] c1; - delete[] c2; - delete[] c3; - delete[] c4; - delete[] c5; - delete[] c6; - delete[] c7; - delete[] c8; - delete[] c9; - delete[] c10; - delete[] c11; - delete[] c12; - delete[] P; - delete[] Pp; - delete[] Bzplus; - delete[] Brplus; - delete[] jlast; - delete[] jlastG; - delete[] G; - delete[] Brem1; - delete[] Bcen1; - delete[] rorem; - delete[] roremG; - delete[] z0remG; + delete[] fC0; + delete[] fC1; + delete[] fc1; + delete[] fc2; + delete[] fc3; + delete[] fc4; + delete[] fc5; + delete[] fc6; + delete[] fc7; + delete[] fc8; + delete[] fc9; + delete[] fc10; + delete[] fc11; + delete[] fc12; + delete[] fP; + delete[] fPp; + delete[] fBzplus; + delete[] fBrplus; + delete[] fjlast; + delete[] fjlastG; + delete[] fG; + delete[] fBrem1; + delete[] fBcen1; + delete[] frorem; + delete[] froremG; + delete[] fz0remG; // Two-dim arrays: - for (int i = 0; i < Ncoil; i++) { - delete[] coil[i]; - delete[] Brem[i]; - delete[] Vrem[i]; - delete[] z0cen[i]; - delete[] rocen[i]; - } - delete[] coil; - delete[] z0cen; - delete[] rocen; - delete[] Brem; - delete[] Vrem; - - for (int g = 0; g < Ng; g++) { - delete[] Cin[g]; - delete[] Line[g]; - delete[] Z[g]; - delete[] BremG[g]; - delete[] z0cenG[g]; - delete[] rocenG[g]; - } - delete[] Cin; - delete[] Line; - delete[] Z; - delete[] BremG; - delete[] z0cenG; - delete[] rocenG; + for (int i = 0; i < fNcoil; i++) { + delete[] fcoil[i]; + delete[] fBrem[i]; + delete[] fVrem[i]; + delete[] fz0cen[i]; + delete[] frocen[i]; + } + delete[] fcoil; + delete[] fz0cen; + delete[] frocen; + delete[] fBrem; + delete[] fVrem; + + for (int g = 0; g < fNg; g++) { + delete[] fCin[g]; + delete[] fLine[g]; + delete[] fZ[g]; + delete[] fBremG[g]; + delete[] fz0cenG[g]; + delete[] frocenG[g]; + } + delete[] fCin; + delete[] fLine; + delete[] fZ; + delete[] fBremG; + delete[] fz0cenG; + delete[] frocenG; // Three-dim arrays: - for (int i = 0; i < Ncoil; i++) { - for (int j = 0; j < Nsp[i]; j++) - delete[] Bcen[i][j]; - delete[] Bcen[i]; - } - delete[] Bcen; - for (int g = 0; g < Ng; g++) { - for (int j = 0; j < NspG[g]; j++) - delete[] BcenG[g][j]; - delete[] BcenG[g]; - } - delete[] BcenG; - delete[] Nsp; - delete[] NspG; - delete[] Nc; + for (int i = 0; i < fNcoil; i++) { + for (int j = 0; j < fNsp[i]; j++) + delete[] fBcen[i][j]; + delete[] fBcen[i]; + } + delete[] fBcen; + for (int g = 0; g < fNg; g++) { + for (int j = 0; j < fNspG[g]; j++) + delete[] fBcenG[g][j]; + delete[] fBcenG[g]; + } + delete[] fBcenG; + delete[] fNsp; + delete[] fNspG; + delete[] fNc; } @@ -148,8 +148,8 @@ void MagfieldCoils::CoilRead() // The coils have rectangular shape cross section with various // symmetry axes. // The data in the file defined by string coilfilename are: - // First line: number of coils (Ncoil). - // Then there are Ncoil number of lines; each line contains: + // First line: number of coils (fNcoil). + // Then there are fNcoil number of lines; each line contains: // cu Cx Cy Cz alpha beta tu L Rmin Rmax // cu: current of coil (A) // Cx: x component of coil center (m) @@ -165,63 +165,63 @@ void MagfieldCoils::CoilRead() // Rmax: outer radius of coil (m). // Positive current: current flow has right-handed screw with the (ux,uy,uz) coil direction vector // (e.g. beta=0, cu positive: magnetic field in +z direction). - // The coil parameters are written into the array coil[Ncoil][14]. - // i: coil index (i=0,...,Ncoil-1) - // coil[i][0]=cu, coil[i][1]=Cx, coil[i][2]=Cy, coil[i][3]=Cz, coil[i][4]=alpha, coil[i][5]=beta, - // coil[i][6]=tu, coil[i][7]=L, coil[i][8]=Rmin, coil[i][9]=Rmax. + // The fcoil parameters are written into the array fcoil[fNcoil][14]. + // i: fcoil index (i=0,...,fNcoil-1) + // fcoil[i][0]=cu, fcoil[i][1]=Cx, fcoil[i][2]=Cy, fcoil[i][3]=Cz, fcoil[i][4]=alpha, fcoil[i][5]=beta, + // fcoil[i][6]=tu, fcoil[i][7]=L, fcoil[i][8]=Rmin, fcoil[i][9]=Rmax. // SI units are used here (A for cu, and m for Cx, Cy, Cz, L, Rmin, Rmax)! ifstream input; - input.open(coilfilename.c_str()); + input.open(fcoilfilename.c_str()); if (!input.is_open()) { puts("Message from function CoilRead1 of class Magfield:"); puts("Cannot open the input source file!"); puts("Program running is stopped !!! "); exit(1); } - input >> Ncoil; - if (Ncoil < 1) { + input >> fNcoil; + if (fNcoil < 1) { puts("Message from function MagfieldCoils::CoilRead:"); - puts("Ncoil is not positive !"); + puts("fNcoil is not positive !"); puts("Program running is stopped !!! "); exit(1); } - // Dynamic memory allocation for the coil parameters: - coil = new double*[Ncoil]; - for (int i = 0; i < Ncoil; i++) { - coil[i] = new double[14]; + // Dynamic memory allocation for the fcoil parameters: + fcoil = new double*[fNcoil]; + for (int i = 0; i < fNcoil; i++) { + fcoil[i] = new double[14]; } // Reading the coil parameters: double cu, Cx, Cy, Cz, alpha, beta, tu, L, Rmin, Rmax; - for (int i = 0; i < Ncoil; i++) { + for (int i = 0; i < fNcoil; i++) { input >> cu >> Cx >> Cy >> Cz >> alpha >> beta >> tu >> L >> Rmin >> Rmax; // Coil parameter determination: - coil[i][0] = cu; - coil[i][1] = Cx; - coil[i][2] = Cy; - coil[i][3] = Cz; - coil[i][4] = alpha; - coil[i][5] = beta; - coil[i][6] = tu; - coil[i][7] = L; - coil[i][8] = Rmin; - coil[i][9] = Rmax; - coil[i][10] = cu * tu / (L * (Rmax - Rmin)); // current density - coil[i][11] = sin(beta / 180. * M_PI) * sin(alpha / 180. * M_PI); // coil direction unit vector comp. ux - coil[i][12] = -sin(beta / 180. * M_PI) * cos(alpha / 180. * M_PI); // coil direction unit vector comp. uy - coil[i][13] = cos(beta / 180. * M_PI); // coil direction unit vector comp. uz + fcoil[i][0] = cu; + fcoil[i][1] = Cx; + fcoil[i][2] = Cy; + fcoil[i][3] = Cz; + fcoil[i][4] = alpha; + fcoil[i][5] = beta; + fcoil[i][6] = tu; + fcoil[i][7] = L; + fcoil[i][8] = Rmin; + fcoil[i][9] = Rmax; + fcoil[i][10] = cu * tu / (L * (Rmax - Rmin)); // current density + fcoil[i][11] = sin(beta / 180. * M_PI) * sin(alpha / 180. * M_PI); // coil direction unit vector comp. ux + fcoil[i][12] = -sin(beta / 180. * M_PI) * cos(alpha / 180. * M_PI); // coil direction unit vector comp. uy + fcoil[i][13] = cos(beta / 180. * M_PI); // coil direction unit vector comp. uz } input.close(); // // Test of coil parameters: - for (int i = 0; i < Ncoil; i++) { - if (coil[i][6] <= 1.e-9 || coil[i][7] <= 1.e-9 || coil[i][8] <= 1.e-9 || coil[i][9] < 1.e-9) { + for (int i = 0; i < fNcoil; i++) { + if (fcoil[i][6] <= 1.e-9 || fcoil[i][7] <= 1.e-9 || fcoil[i][8] <= 1.e-9 || fcoil[i][9] < 1.e-9) { printf("Message from function MagfieldCoils::CoilRead1:\ non-positive coil parameters: tu, Rmin, Rmax or L !!!\ Program running is stopped !!! i= %9i \n", i); exit(1); } - if (coil[i][8] >= coil[i][9]) { + if (fcoil[i][8] >= fcoil[i][9]) { printf("Message from function MagfieldCoils::CoilRead1:\ Rmin>=Rmax!!!\ Program running is stopped !!! i= %9i \n", @@ -234,37 +234,37 @@ void MagfieldCoils::CoilRead() // // Coil symmetry group calculations: // - G = new int[Ncoil]; // coil with index i is in symmetry group g=G[i] - for (int i = 0; i < Ncoil; i++) - G[i] = -1; + fG = new int[fNcoil]; // coil with index i is in symmetry group g=G[i] + for (int i = 0; i < fNcoil; i++) + fG[i] = -1; int g = -1; // Temporary arrays: - int* Nctemp = new int[Ncoil]; // max. number of groups=Ncoil - int** Cintemp = new int*[Ncoil]; - for (int i = 0; i < Ncoil; i++) - Cintemp[i] = new int[Ncoil]; // max. number of coils in a group=Ncoil - double** Linetemp = new double*[Ncoil]; - for (int i = 0; i < Ncoil; i++) + int* Nctemp = new int[fNcoil]; // max. number of groups=fNcoil + int** Cintemp = new int*[fNcoil]; + for (int i = 0; i < fNcoil; i++) + Cintemp[i] = new int[fNcoil]; // max. number of coils in a group=fNcoil + double** Linetemp = new double*[fNcoil]; + for (int i = 0; i < fNcoil; i++) Linetemp[i] = new double[6]; - double** Ztemp = new double*[Ncoil]; - for (int i = 0; i < Ncoil; i++) - Ztemp[i] = new double[Ncoil]; // max. number of coils in a group=Ncoil + double** Ztemp = new double*[fNcoil]; + for (int i = 0; i < fNcoil; i++) + Ztemp[i] = new double[fNcoil]; // max. number of coils in a group=fNcoil // Starting the symmetry group calculation: - for (int i = 0; i < Ncoil; i++) { + for (int i = 0; i < fNcoil; i++) { int c; - if (G[i] == -1) // starting a new symmetry group + if (fG[i] == -1) // starting a new symmetry group { g += 1; - Ng = g + 1; + fNg = g + 1; Nctemp[g] = 1; // Coil i defines the symmetry axis of the new symmetry group. // The symmetry axis line is defined by the coil i center C[3] and the coil i direction unit vector u[3]; // a point of the line: C[k]+z*u[k], k=0,1,2 (z: line point parameter). double C[3], u[3], alpha, beta; for (int k = 0; k < 3; k++) - C[k] = coil[i][1 + k]; // coil center - alpha = coil[i][4]; - beta = coil[i][5]; + C[k] = fcoil[i][1 + k]; // coil center + alpha = fcoil[i][4]; + beta = fcoil[i][5]; u[0] = sin(beta / 180. * M_PI) * sin(alpha / 180. * M_PI); u[1] = -sin(beta / 180. * M_PI) * cos(alpha / 180. * M_PI); u[2] = cos(beta / 180. * M_PI); @@ -277,24 +277,24 @@ void MagfieldCoils::CoilRead() Cintemp[g][0] = i; // z coordinate of first coil (local coil index c=0) in symmetry group g: Z[g][0]=0.: Ztemp[g][0] = 0.; - // If Ncoil=1: end of the i-loop - if (Ncoil == 1) + // If fNcoil=1: end of the i-loop + if (fNcoil == 1) break; // Now we start to collect the coils which belong to this symmetry group: - for (int ic = i + 1; ic < Ncoil; ic++) // ic is here global coil index + for (int ic = i + 1; ic < fNcoil; ic++) // ic is here global coil index { double vecA[3], vecB[3], dA, dB, ZA, ZB, A[3], B[3], L, vec[3], uc[3], alphac, betac; // Coil ic end points: A[3], B[3] - alphac = coil[ic][4]; - betac = coil[ic][5]; + alphac = fcoil[ic][4]; + betac = fcoil[ic][5]; uc[0] = sin(betac / 180. * M_PI) * sin(alphac / 180. * M_PI); uc[1] = -sin(betac / 180. * M_PI) * cos(alphac / 180. * M_PI); uc[2] = cos(betac / 180. * M_PI); - L = coil[ic][7]; + L = fcoil[ic][7]; ZA = ZB = 0.; for (int k = 0; k < 3; k++) { - A[k] = coil[ic][1 + k] - uc[k] * L / 2.; - B[k] = coil[ic][1 + k] + uc[k] * L / 2.; // coil ic endpoints + A[k] = fcoil[ic][1 + k] - uc[k] * L / 2.; + B[k] = fcoil[ic][1 + k] + uc[k] * L / 2.; // coil ic endpoints ZA += (A[k] - C[k]) * u[k]; ZB += (B[k] - C[k]) * u[k]; } @@ -304,9 +304,9 @@ void MagfieldCoils::CoilRead() } dA = Mag(vecA); dB = Mag(vecB); - if (dA < epstol && dB < epstol) // coil ic is put to this symmetry group + if (dA < fepstol && dB < fepstol) // coil ic is put to this symmetry group { - G[ic] = g; + fG[ic] = g; Nctemp[g] += 1; c = Nctemp[g] - 1; Cintemp[g][c] = ic; @@ -316,57 +316,57 @@ void MagfieldCoils::CoilRead() A[k] = C[k] + ZA * u[k]; B[k] = C[k] + ZB * u[k]; vec[k] = B[k] - A[k]; - coil[ic][1 + k] = (A[k] + B[k]) / 2.; + fcoil[ic][1 + k] = (A[k] + B[k]) / 2.; } // alphac-> alpha, betac-> beta; // If uc is opposite to u: change of current sign for (int k = 0; k < 3; k++) vec[k] = u[k] - uc[k]; if (Mag(vec) > 0.5) - coil[ic][0] *= -1.; - coil[ic][4] = alpha; - coil[ic][5] = beta; - coil[ic][11] = + fcoil[ic][0] *= -1.; + fcoil[ic][4] = alpha; + fcoil[ic][5] = beta; + fcoil[ic][11] = sin(beta / 180. * M_PI) * sin(alpha / 180. * M_PI); // coil direction unit vector comp. ux - coil[ic][12] = + fcoil[ic][12] = -sin(beta / 180. * M_PI) * cos(alpha / 180. * M_PI); // coil direction unit vector comp. uy - coil[ic][13] = cos(beta / 180. * M_PI); // coil direction unit vector comp. uz + fcoil[ic][13] = cos(beta / 180. * M_PI); // coil direction unit vector comp. uz // After that, all coils within a symmetry group have the same coil directions (alpha, beta). } // end of epstol-if } // end of symmetry group loop ic } // end of symmetry group if } // end of coil index loop i // Nc, Cin, Line and Z array calculation: - Nc = new int[Ng]; - for (int g = 0; g < Ng; g++) - Nc[g] = Nctemp[g]; // number of coils in symmetry group g - Cin = new int*[Ng]; - for (int g = 0; g < Ng; g++) { - Cin[g] = new int[Nc[g]]; - for (int c = 0; c < Nc[g]; c++) - Cin[g][c] = Cintemp[g][c]; // coil indices i in symmetry group g - } - Line = new double*[Ng]; - for (int g = 0; g < Ng; g++) { - Line[g] = new double[6]; + fNc = new int[fNg]; + for (int g = 0; g < fNg; g++) + fNc[g] = Nctemp[g]; // number of coils in symmetry group g + fCin = new int*[fNg]; + for (int g = 0; g < fNg; g++) { + fCin[g] = new int[fNc[g]]; + for (int c = 0; c < fNc[g]; c++) + fCin[g][c] = Cintemp[g][c]; // coil indices i in symmetry group g + } + fLine = new double*[fNg]; + for (int g = 0; g < fNg; g++) { + fLine[g] = new double[6]; for (int k = 0; k < 6; k++) - Line[g][k] = Linetemp[g][k]; // line point and line direction unit vector + fLine[g][k] = Linetemp[g][k]; // line point and line direction unit vector } - Z = new double*[Ng]; - for (int g = 0; g < Ng; g++) { - Z[g] = new double[Nc[g]]; - for (int c = 0; c < Nc[g]; c++) - Z[g][c] = Ztemp[g][c]; // local z coordinate of coil center in symmetry group g; Z[g][0]=0. + fZ = new double*[fNg]; + for (int g = 0; g < fNg; g++) { + fZ[g] = new double[fNc[g]]; + for (int c = 0; c < fNc[g]; c++) + fZ[g][c] = Ztemp[g][c]; // local z coordinate of coil center in symmetry group g; Z[g][0]=0. } // Delete Nctemp, Cintemp, Linetemp, Ztemp: delete[] Nctemp; - for (int i = 0; i < Ncoil; i++) + for (int i = 0; i < fNcoil; i++) delete[] Cintemp[i]; delete[] Cintemp; - for (int i = 0; i < Ncoil; i++) + for (int i = 0; i < fNcoil; i++) delete[] Linetemp[i]; delete[] Linetemp; - for (int i = 0; i < Ncoil; i++) + for (int i = 0; i < fNcoil; i++) delete[] Ztemp[i]; delete[] Ztemp; @@ -378,30 +378,30 @@ void MagfieldCoils::CoilRead() void MagfieldCoils::CoilGroupWrite() { // This function writes the coil and symmetry group parameters into the data file dirname+objectname+coilgroup.txt. - string filename = dirname + objectname + "coilgroup.txt"; + string filename = fdirname + fobjectname + "coilgroup.txt"; ofstream output; output.precision(16); output.open(filename.c_str()); - output << setw(15) << Ncoil << setw(15) << Ng << endl; // number of coils and groups - output << setw(15) << Nelliptic << setw(15) << nmax << endl; - for (int i = 0; i < Ncoil; i++) { + output << setw(15) << fNcoil << setw(15) << fNg << endl; // number of coils and groups + output << setw(15) << fNelliptic << setw(15) << fnmax << endl; + for (int i = 0; i < fNcoil; i++) { output << setw(10) << i << endl; - output << scientific << setw(26) << coil[i][0] << setw(26) << coil[i][1] << setw(26) << coil[i][2] << setw(26) - << coil[i][3] << endl; - output << scientific << setw(26) << coil[i][4] << setw(26) << coil[i][5] << setw(26) << coil[i][6] << endl; - output << scientific << setw(26) << coil[i][7] << setw(26) << coil[i][8] << setw(26) << coil[i][9] << setw(26) - << coil[i][10] << endl; - output << scientific << setw(26) << coil[i][11] << setw(26) << coil[i][12] << setw(26) << coil[i][13] << endl; - } - for (int g = 0; g < Ng; g++) { - output << setw(12) << g << setw(12) << Nc[g] << endl; // number of coils in symmetry group g + output << scientific << setw(26) << fcoil[i][0] << setw(26) << fcoil[i][1] << setw(26) << fcoil[i][2] << setw(26) + << fcoil[i][3] << endl; + output << scientific << setw(26) << fcoil[i][4] << setw(26) << fcoil[i][5] << setw(26) << fcoil[i][6] << endl; + output << scientific << setw(26) << fcoil[i][7] << setw(26) << fcoil[i][8] << setw(26) << fcoil[i][9] << setw(26) + << fcoil[i][10] << endl; + output << scientific << setw(26) << fcoil[i][11] << setw(26) << fcoil[i][12] << setw(26) << fcoil[i][13] << endl; + } + for (int g = 0; g < fNg; g++) { + output << setw(12) << g << setw(12) << fNc[g] << endl; // number of coils in symmetry group g // Cin[g][c]: // global coil index corresponding to local coil index c in group g // Z[g][c]: local z coordinate of coil center in symmetry group g; c=0,..., Nc[g]-1; Z[g][0]=0 // Symmetry axis line point: Line[g][0], Line[g][1], Line[g][2]; line direction unit vector: Line[g][3], Line[g][4], Line[g][5]; - output << scientific << setw(26) << Line[g][0] << setw(26) << Line[g][1] << setw(26) << Line[g][2] << endl; - output << scientific << setw(26) << Line[g][3] << setw(26) << Line[g][4] << setw(26) << Line[g][5] << endl; - for (int c = 0; c < Nc[g]; c++) - output << setw(12) << Cin[g][c] << scientific << setw(26) << Z[g][c] << endl; + output << scientific << setw(26) << fLine[g][0] << setw(26) << fLine[g][1] << setw(26) << fLine[g][2] << endl; + output << scientific << setw(26) << fLine[g][3] << setw(26) << fLine[g][4] << setw(26) << fLine[g][5] << endl; + for (int c = 0; c < fNc[g]; c++) + output << setw(12) << fCin[g][c] << scientific << setw(26) << fZ[g][c] << endl; } output.close(); return; @@ -414,7 +414,7 @@ void MagfieldCoils::CoilGroupWrite() void MagfieldCoils::CoilGroupRead() { // This function reads the coil parameters from the data file dirname+objectname+coilgroup.txt. - string filename = dirname + objectname + "coilgroup.txt"; + string filename = fdirname + fobjectname + "coilgroup.txt"; ifstream input; input.open(filename.c_str()); if (!input.is_open()) { @@ -424,57 +424,57 @@ void MagfieldCoils::CoilGroupRead() exit(1); } // We read the coil parameters: - input >> Ncoil >> Ng; - if (Ncoil < 1 || Ng < 1) { + input >> fNcoil >> fNg; + if (fNcoil < 1 || fNg < 1) { printf("Message from function MagfieldCoils::CoilRead2:\ - Ncoil<1 or Ng<1 !!!\ - Computation is stopped !!! Ncoil, Ng= %9i %9i \n", - Ncoil, - Ng); + fNcoil<1 or Ng<1 !!!\ + Computation is stopped !!! fNcoil, Ng= %9i %9i \n", + fNcoil, + fNg); exit(1); } - input >> Nelliptic >> nmax; + input >> fNelliptic >> fnmax; // Dynamic memory allocation for the coil and symmetry group parameters: - coil = new double*[Ncoil]; - for (int i = 0; i < Ncoil; i++) { - coil[i] = new double[14]; + fcoil = new double*[fNcoil]; + for (int i = 0; i < fNcoil; i++) { + fcoil[i] = new double[14]; } // We read the coil parameters: int ix; - for (int i = 0; i < Ncoil; i++) { + for (int i = 0; i < fNcoil; i++) { input >> ix; - input >> coil[i][0] >> coil[i][1] >> coil[i][2] >> coil[i][3]; - input >> coil[i][4] >> coil[i][5] >> coil[i][6]; - input >> coil[i][7] >> coil[i][8] >> coil[i][9] >> coil[i][10]; - input >> coil[i][11] >> coil[i][12] >> coil[i][13]; + input >> fcoil[i][0] >> fcoil[i][1] >> fcoil[i][2] >> fcoil[i][3]; + input >> fcoil[i][4] >> fcoil[i][5] >> fcoil[i][6]; + input >> fcoil[i][7] >> fcoil[i][8] >> fcoil[i][9] >> fcoil[i][10]; + input >> fcoil[i][11] >> fcoil[i][12] >> fcoil[i][13]; } // Dynamic memory allocation for the symmetry group arrays: - Nc = new int[Ng]; - Line = new double*[Ng]; - Cin = new int*[Ng]; - Z = new double*[Ng]; + fNc = new int[fNg]; + fLine = new double*[fNg]; + fCin = new int*[fNg]; + fZ = new double*[fNg]; // We read the symmetry group parameters: int gx; - for (int g = 0; g < Ng; g++) { - input >> gx >> Nc[g]; + for (int g = 0; g < fNg; g++) { + input >> gx >> fNc[g]; // Cin[g][c]: // global coil index corresponding to local coil index c in group g // Z[g][c]: local z coordinate of coil center in symmetry group g; c=0,..., Nc[g]-1; Z[g][0]=0 // Symmetry axis line point: Line[g][0], Line[g][1], Line[g][2]; line direction unit vector: Line[g][3], Line[g][4], Line[g][5]; - Line[g] = new double[6]; - input >> Line[g][0] >> Line[g][1] >> Line[g][2]; - input >> Line[g][3] >> Line[g][4] >> Line[g][5]; - Cin[g] = new int[Nc[g]]; - Z[g] = new double[Nc[g]]; - for (int c = 0; c < Nc[g]; c++) - input >> Cin[g][c] >> Z[g][c]; + fLine[g] = new double[6]; + input >> fLine[g][0] >> fLine[g][1] >> fLine[g][2]; + input >> fLine[g][3] >> fLine[g][4] >> fLine[g][5]; + fCin[g] = new int[fNc[g]]; + fZ[g] = new double[fNc[g]]; + for (int c = 0; c < fNc[g]; c++) + input >> fCin[g][c] >> fZ[g][c]; } input.close(); - G = new int[Ncoil]; - Brem1 = new double[nmax + 1]; - Bcen1 = new double[nmax + 1]; + fG = new int[fNcoil]; + fBrem1 = new double[fnmax + 1]; + fBcen1 = new double[fnmax + 1]; return; } @@ -487,7 +487,7 @@ void MagfieldCoils::MagsourceRead() // This function reads the source constants from the data files // dirname+objectname+magsource_central.txt and dirname+objectname+magsource_remote.txt. // Part 1: reading the remote source constants: - string filename = dirname + objectname + "magsource_remote.txt"; + string filename = fdirname + fobjectname + "magsource_remote.txt"; ifstream input; input.open(filename.c_str()); if (!input.is_open()) { @@ -497,45 +497,45 @@ void MagfieldCoils::MagsourceRead() exit(1); } // Dynamic memory allocations (rorem, Brem): - rorem = new double[Ncoil]; - Brem = new double*[Ncoil]; - for (int i = 0; i < Ncoil; i++) - Brem[i] = new double[nmax + 1]; + frorem = new double[fNcoil]; + fBrem = new double*[fNcoil]; + for (int i = 0; i < fNcoil; i++) + fBrem[i] = new double[fnmax + 1]; // Coil index loop: int ix, nx; - for (int i = 0; i < Ncoil; i++) { - input >> ix >> rorem[i]; - for (int n = 2; n <= nmax; n++) - input >> nx >> Brem[i][n]; + for (int i = 0; i < fNcoil; i++) { + input >> ix >> frorem[i]; + for (int n = 2; n <= fnmax; n++) + input >> nx >> fBrem[i][n]; } // Dynamic memory allocations (Vrem): - Vrem = new double*[Ncoil]; - for (int i = 0; i < Ncoil; i++) - Vrem[i] = new double[nmax + 1]; - for (int i = 0; i < Ncoil; i++) { + fVrem = new double*[fNcoil]; + for (int i = 0; i < fNcoil; i++) + fVrem[i] = new double[fnmax + 1]; + for (int i = 0; i < fNcoil; i++) { input >> ix; - for (int n = 0; n <= nmax; n++) - input >> nx >> Vrem[i][n]; + for (int n = 0; n <= fnmax; n++) + input >> nx >> fVrem[i][n]; } // Dynamic memory allocations (z0remG, roremG, BremG): - z0remG = new double[Ng]; - roremG = new double[Ng]; - BremG = new double*[Ng]; - for (int g = 0; g < Ng; g++) - BremG[g] = new double[nmax + 1]; + fz0remG = new double[fNg]; + froremG = new double[fNg]; + fBremG = new double*[fNg]; + for (int g = 0; g < fNg; g++) + fBremG[g] = new double[fnmax + 1]; // Group index loop: int gx, Ncg; - for (int g = 0; g < Ng; g++) { - if (Nc[g] == 1) + for (int g = 0; g < fNg; g++) { + if (fNc[g] == 1) continue; - input >> gx >> Ncg >> z0remG[g] >> roremG[g]; - for (int n = 2; n <= nmax; n++) - input >> nx >> BremG[g][n]; + input >> gx >> Ncg >> fz0remG[g] >> froremG[g]; + for (int n = 2; n <= fnmax; n++) + input >> nx >> fBremG[g][n]; } input.close(); // // Part 2: reading the central source coefficients: - filename = dirname + objectname + "magsource_central.txt"; + filename = fdirname + fobjectname + "magsource_central.txt"; input.open(filename.c_str()); if (!input.is_open()) { puts("Message from function MagfieldCoils::MagsourceRead:"); @@ -545,53 +545,53 @@ void MagfieldCoils::MagsourceRead() } int jx; // Dynamic memory allocations (Nsp, z0cen, rocen, Bcen): - Nsp = new int[Ncoil]; - z0cen = new double*[Ncoil]; - rocen = new double*[Ncoil]; - Bcen = new double**[Ncoil]; + fNsp = new int[fNcoil]; + fz0cen = new double*[fNcoil]; + frocen = new double*[fNcoil]; + fBcen = new double**[fNcoil]; // Coil index loop: string text; // input >> text; - for (int i = 0; i < Ncoil; i++) { - input >> ix >> Nsp[i]; + for (int i = 0; i < fNcoil; i++) { + input >> ix >> fNsp[i]; // Source point loop: - z0cen[i] = new double[Nsp[i]]; - rocen[i] = new double[Nsp[i]]; - Bcen[i] = new double*[Nsp[i]]; - for (int j = 0; j < Nsp[i]; j++) { - input >> jx >> z0cen[i][j] >> rocen[i][j]; - Bcen[i][j] = new double[nmax + 1]; - for (int n = 0; n <= nmax; n++) - input >> nx >> Bcen[i][j][n]; + fz0cen[i] = new double[fNsp[i]]; + frocen[i] = new double[fNsp[i]]; + fBcen[i] = new double*[fNsp[i]]; + for (int j = 0; j < fNsp[i]; j++) { + input >> jx >> fz0cen[i][j] >> frocen[i][j]; + fBcen[i][j] = new double[fnmax + 1]; + for (int n = 0; n <= fnmax; n++) + input >> nx >> fBcen[i][j][n]; } } // Dynamic memory allocations (NspG, z0cenG, rocenG, BcenG): - NspG = new int[Ng]; - z0cenG = new double*[Ng]; - rocenG = new double*[Ng]; - BcenG = new double**[Ng]; + fNspG = new int[fNg]; + fz0cenG = new double*[fNg]; + frocenG = new double*[fNg]; + fBcenG = new double**[fNg]; // Group index loop: // input >> text; - for (int g = 0; g < Ng; g++) { - if (Nc[g] == 1) // no group calculation if the group has only 1 coil + for (int g = 0; g < fNg; g++) { + if (fNc[g] == 1) // no group calculation if the group has only 1 coil { - z0cenG[g] = new double[1]; - rocenG[g] = new double[1]; - NspG[g] = 1; - BcenG[g] = new double*[1]; - BcenG[g][0] = new double[nmax + 1]; + fz0cenG[g] = new double[1]; + frocenG[g] = new double[1]; + fNspG[g] = 1; + fBcenG[g] = new double*[1]; + fBcenG[g][0] = new double[fnmax + 1]; continue; } - input >> gx >> NspG[g]; + input >> gx >> fNspG[g]; // Source point loop: - z0cenG[g] = new double[NspG[g]]; - rocenG[g] = new double[NspG[g]]; - BcenG[g] = new double*[NspG[g]]; - for (int j = 0; j < NspG[g]; j++) { - input >> jx >> z0cenG[g][j] >> rocenG[g][j]; - BcenG[g][j] = new double[nmax + 1]; - for (int n = 0; n <= nmax; n++) - input >> nx >> BcenG[g][j][n]; + fz0cenG[g] = new double[fNspG[g]]; + frocenG[g] = new double[fNspG[g]]; + fBcenG[g] = new double*[fNspG[g]]; + for (int j = 0; j < fNspG[g]; j++) { + input >> jx >> fz0cenG[g][j] >> frocenG[g][j]; + fBcenG[g][j] = new double[fnmax + 1]; + for (int n = 0; n <= fnmax; n++) + input >> nx >> fBcenG[g][j][n]; } } input.close(); @@ -605,31 +605,31 @@ void MagfieldCoils::DynamicMemoryAllocations() // This function makes dynamic memory allocations for the variables // C0, C1, c0 to c12, P, Pp, Bzplus, Brplus. // It also computes the values C0[n], C1[n], c0[n] -- c12[n] (n=2,...,nmax). - C0 = new double[nmax + 1]; - C1 = new double[nmax + 1]; - c1 = new double[nmax + 1]; - c2 = new double[nmax + 1]; - c3 = new double[nmax + 1]; - c4 = new double[nmax + 1]; - c5 = new double[nmax + 1]; - c6 = new double[nmax + 1]; - c7 = new double[nmax + 1]; - c8 = new double[nmax + 1]; - c9 = new double[nmax + 1]; - c10 = new double[nmax + 1]; - c11 = new double[nmax + 1]; - c12 = new double[nmax + 1]; - - C0[1] = 1.; - for (int n = 2; n <= nmax; n++) { - C0[n] = 1. / (1. * n); - C1[n] = 1. / (1. * (n - 1.)); - c1[n] = (2. * n - 1.) / (1. * n); - c2[n] = (n - 1.) / (1. * n); - c3[n] = (2. * n - 1.) / (1. * (n - 1.)); - c4[n] = (1. * n) / (1. * (n - 1.)); - c5[n] = (1.) / (n * 1.); - c6[n] = (1.) / (n + 1.); + fC0 = new double[fnmax + 1]; + fC1 = new double[fnmax + 1]; + fc1 = new double[fnmax + 1]; + fc2 = new double[fnmax + 1]; + fc3 = new double[fnmax + 1]; + fc4 = new double[fnmax + 1]; + fc5 = new double[fnmax + 1]; + fc6 = new double[fnmax + 1]; + fc7 = new double[fnmax + 1]; + fc8 = new double[fnmax + 1]; + fc9 = new double[fnmax + 1]; + fc10 = new double[fnmax + 1]; + fc11 = new double[fnmax + 1]; + fc12 = new double[fnmax + 1]; + + fC0[1] = 1.; + for (int n = 2; n <= fnmax; n++) { + fC0[n] = 1. / (1. * n); + fC1[n] = 1. / (1. * (n - 1.)); + fc1[n] = (2. * n - 1.) / (1. * n); + fc2[n] = (n - 1.) / (1. * n); + fc3[n] = (2. * n - 1.) / (1. * (n - 1.)); + fc4[n] = (1. * n) / (1. * (n - 1.)); + fc5[n] = (1.) / (n * 1.); + fc6[n] = (1.) / (n + 1.); if (n >= 4) { int m = n - 2; double M = (m + 1.) * (m + 2.) * (2. * m - 1.); @@ -640,26 +640,26 @@ void MagfieldCoils::DynamicMemoryAllocations() double Bp = 2. * m * (m + 2.) * (2. * m - 1.) - 3.; double C = m * (m - 1.) * (2. * m + 3.); double Cp = m * (m + 1.) * (2. * m + 3.); - c7[n] = A / M; - c8[n] = B / M; - c9[n] = C / M; - c10[n] = Ap / Mp; - c11[n] = Bp / Mp; - c12[n] = Cp / Mp; + fc7[n] = A / M; + fc8[n] = B / M; + fc9[n] = C / M; + fc10[n] = Ap / Mp; + fc11[n] = Bp / Mp; + fc12[n] = Cp / Mp; } } // Dynamic memory for P, Pp, Bzplus, Brplus, jlast, jlastG: - P = new double[nmax + 2]; - Pp = new double[nmax + 2]; - Bzplus = new double[nmax + 1]; - Brplus = new double[nmax + 1]; + fP = new double[fnmax + 2]; + fPp = new double[fnmax + 2]; + fBzplus = new double[fnmax + 1]; + fBrplus = new double[fnmax + 1]; // - jlast = new int[Ncoil]; - for (int i = 0; i < Ncoil; i++) - jlast[i] = -2; - jlastG = new int[Ng]; - for (int g = 0; g < Ng; g++) - jlastG[g] = -2; + fjlast = new int[fNcoil]; + for (int i = 0; i < fNcoil; i++) + fjlast[i] = -2; + fjlastG = new int[fNg]; + for (int g = 0; g < fNg; g++) + fjlastG[g] = -2; } ///////////////////////////////////////////////////// @@ -686,12 +686,12 @@ void MagfieldCoils::Magfield2EllipticCoil(int i, double z, double r, double& Bz, const double mu0 = 4. * M_PI * 1.e-7; double x[2][1001], w[2][1001]; // nodes and weights // Coil parameters: - L = coil[i][7]; // coil length + L = fcoil[i][7]; // coil length Zmin = -L / 2.; Zmax = L / 2.; // coil endpoints relative to coil center - Rmin = coil[i][8]; - Rmax = coil[i][9]; // coil inner and outer radii - sigma = coil[i][10]; // current density + Rmin = fcoil[i][8]; + Rmax = fcoil[i][9]; // coil inner and outer radii + sigma = fcoil[i][10]; // current density // Field point test: if (fabs(z - Zmin) < 1.e-8 && r >= Rmin - (Rmax - Rmin) * 1.e-8 && r <= Rmax + (Rmax - Rmin) * 1.e-8) z = Zmin - 1.e-8; @@ -699,7 +699,7 @@ void MagfieldCoils::Magfield2EllipticCoil(int i, double z, double r, double& Bz, z = Zmax + 1.e-8; // Radial numerical integration node definition: int N[2]; - N[0] = N[1] = Nelliptic; + N[0] = N[1] = fNelliptic; // // R-integration limits: int M; @@ -783,11 +783,11 @@ void MagfieldCoils::MagfieldEllipticCoil(int i, const double* P, double* B) // We define a local coordinate system: // origo at the coil center, local z axis parallel to coil axis, in u vector direction. // Coil direction u and center C: - u[0] = coil[i][11]; - u[1] = coil[i][12]; - u[2] = coil[i][13]; + u[0] = fcoil[i][11]; + u[1] = fcoil[i][12]; + u[2] = fcoil[i][13]; for (int k = 0; k <= 2; k++) - C[k] = coil[i][1 + k]; // coil center + C[k] = fcoil[i][1 + k]; // coil center // Local z and r coordinates of the field point P: for (int k = 0; k <= 2; k++) Ploc[k] = P[k] - C[k]; @@ -820,7 +820,7 @@ void MagfieldCoils::MagfieldElliptic(const double* P, double* B) double Bi[3]; for (int k = 0; k <= 2; k++) B[k] = 0.; - for (int i = 0; i < Ncoil; i++) { + for (int i = 0; i < fNcoil; i++) { MagfieldEllipticCoil(i, P, Bi); for (int k = 0; k <= 2; k++) B[k] += Bi[k]; @@ -845,7 +845,7 @@ bool MagfieldCoils::Magfield(const double* P, double* B) bool magfield = true; for (int k = 0; k < 3; k++) B[k] = 0.; - for (int g = 0; g < Ng; g++) { + for (int g = 0; g < fNg; g++) { bool magfieldgroup = MagfieldGroup(g, P, Bgroup); if (magfieldgroup == false) magfield = false; @@ -869,8 +869,8 @@ double MagfieldCoils::Funrorem(int i, double z0) // rorem = maximum distance of the axis point z0 from the coil winding (outer corner points). { double L, Rmax, Zmin, Zmax; - L = coil[i][7]; - Rmax = coil[i][9]; + L = fcoil[i][7]; + Rmax = fcoil[i][9]; Zmin = -L / 2.; Zmax = L / 2.; // coil endpoint Z values defined relative to the coil center if (z0 > 0.) @@ -897,12 +897,12 @@ void MagfieldCoils::Magsource2RemoteCoil(int i, double z0, double rorem) double x[1001], w[1001]; // nodes and weights double L, sigma, Zmin, Zmax, Rmin, Rmax; // Coil parameters: - L = coil[i][7]; // coil length + L = fcoil[i][7]; // coil length Zmin = -L / 2.; Zmax = L / 2.; // coil endpoints relative to coil center - Rmin = coil[i][8]; - Rmax = coil[i][9]; // coil inner and outer radii - sigma = coil[i][10]; // current density + Rmin = fcoil[i][8]; + Rmax = fcoil[i][9]; // coil inner and outer radii + sigma = fcoil[i][10]; // current density // Radial integration number N: int N; double ratio = (Rmax - Rmin) / Rmin; @@ -913,10 +913,10 @@ void MagfieldCoils::Magsource2RemoteCoil(int i, double z0, double rorem) if (N > 1000) N = 1000; // Initialization of Brem1[n] and Pp[n]: - for (int n = 2; n <= nmax; n++) - Brem1[n] = 0.; - Pp[0] = 0.; - Pp[1] = 1.; + for (int n = 2; n <= fnmax; n++) + fBrem1[n] = 0.; + fPp[0] = 0.; + fPp[1] = 1.; // C: double C = mu0 * sigma / (2. * rorem * rorem); // Radial integration nodes and weights: @@ -935,10 +935,10 @@ void MagfieldCoils::Magsource2RemoteCoil(int i, double z0, double rorem) double a = ros / rorem; double b = a; double d = C * R2 * sign[iz] * w[integ]; - for (int n = 2; n <= nmax; n++) { - double usPp = us * Pp[n - 1]; - Pp[n] = 2. * usPp - Pp[n - 2] + (usPp - Pp[n - 2]) * C1[n]; - Brem1[n] += d * c6[n] * b * Pp[n]; + for (int n = 2; n <= fnmax; n++) { + double usPp = us * fPp[n - 1]; + fPp[n] = 2. * usPp - fPp[n - 2] + (usPp - fPp[n - 2]) * fC1[n]; + fBrem1[n] += d * fc6[n] * b * fPp[n]; b *= a; } } @@ -952,41 +952,41 @@ void MagfieldCoils::Magsource2RemoteCoil(int i, double z0, double rorem) // (3-dim. case) void MagfieldCoils::MagsourceRemoteCoils() -// This function computes the Brem[i][n] (i=0,...,Ncoil-1; n=2,...,nmax) +// This function computes the Brem[i][n] (i=0,...,fNcoil-1; n=2,...,nmax) // remote source constants for all coils. // Number of remote source points for a fixed coil is 1 (coil center). // The results are written into the data file dirname+objectname+magsource_remote.txt: -// Ncoil: number of coils, +// fNcoil: number of coils, // nmax: number of source coefficients for a fixed coil is nmax+1 (n=0,...,nmax), // i: coil index, // rorem[i]: remote convergence radius for coil i. { // Output to file dirname+objectname+magsource_remote.txt : - string filename = dirname + objectname + "magsource_remote.txt"; + string filename = fdirname + fobjectname + "magsource_remote.txt"; ofstream output; output.precision(16); output.open(filename.c_str()); // Dynamic memory allocations (Brem1, rorem, Brem): - Brem1 = new double[nmax + 1]; - rorem = new double[Ncoil]; - Brem = new double*[Ncoil]; - for (int i = 0; i < Ncoil; i++) - Brem[i] = new double[nmax + 1]; + fBrem1 = new double[fnmax + 1]; + frorem = new double[fNcoil]; + fBrem = new double*[fNcoil]; + for (int i = 0; i < fNcoil; i++) + fBrem[i] = new double[fnmax + 1]; // Coil index loop: - for (int i = 0; i < Ncoil; i++) { + for (int i = 0; i < fNcoil; i++) { // Brem1[n] calc.: - rorem[i] = Funrorem(i, 0.); - Magsource2RemoteCoil(i, 0., rorem[i]); + frorem[i] = Funrorem(i, 0.); + Magsource2RemoteCoil(i, 0., frorem[i]); // Output to file dirname+objectname+magsource_remote.txt (rorem, Brem1[n]) - output << scientific << setw(7) << i << setw(28) << rorem[i] << endl; - for (int n = 2; n <= nmax; n++) { - if (fabs(Brem1[n]) < 1.e-30) - Brem1[n] = 0.; + output << scientific << setw(7) << i << setw(28) << frorem[i] << endl; + for (int n = 2; n <= fnmax; n++) { + if (fabs(fBrem1[n]) < 1.e-30) + fBrem1[n] = 0.; if (n % 2 == 1) - output << scientific << setprecision(1) << setw(15) << n << setw(28) << Brem1[n] << endl; + output << scientific << setprecision(1) << setw(15) << n << setw(28) << fBrem1[n] << endl; else - output << scientific << setprecision(16) << setw(15) << n << setw(28) << Brem1[n] << endl; - Brem[i][n] = Brem1[n]; + output << scientific << setprecision(16) << setw(15) << n << setw(28) << fBrem1[n] << endl; + fBrem[i][n] = fBrem1[n]; } } output.close(); @@ -1009,40 +1009,40 @@ void MagfieldCoils::MagsourceRemoteGroups() // roremG[g]: remote convergence radius for group g at the remote source point z0remG[g]. { // Output to file dirname+objectname+magsource_remote.txt : - string filename = dirname + objectname + "magsource_remote.txt"; + string filename = fdirname + fobjectname + "magsource_remote.txt"; ofstream output; output.precision(16); output.open(filename.c_str(), ios::app); // Dynamic memory allocations (z0remG, roremG, BremG): - z0remG = new double[Ng]; - roremG = new double[Ng]; - BremG = new double*[Ng]; - for (int g = 0; g < Ng; g++) - BremG[g] = new double[nmax + 1]; + fz0remG = new double[fNg]; + froremG = new double[fNg]; + fBremG = new double*[fNg]; + for (int g = 0; g < fNg; g++) + fBremG[g] = new double[fnmax + 1]; // Group index loop: - for (int g = 0; g < Ng; g++) { - if (Nc[g] == 1) // no group calculation if the group has only 1 coil + for (int g = 0; g < fNg; g++) { + if (fNc[g] == 1) // no group calculation if the group has only 1 coil continue; // z0remG[g], roremG[g] calculation: RemoteSourcepointGroup(g); // Output to file dirname+objectname+magsource_remote.txt (g, Nc[g], ) - output << scientific << setw(10) << g << setw(10) << Nc[g] << setw(26) << z0remG[g] << setw(26) << roremG[g] + output << scientific << setw(10) << g << setw(10) << fNc[g] << setw(26) << fz0remG[g] << setw(26) << froremG[g] << endl; // BremG[g][n] calculation by local coil index loop: - for (int n = 2; n <= nmax; n++) - BremG[g][n] = 0.; + for (int n = 2; n <= fnmax; n++) + fBremG[g][n] = 0.; // Local coil index loop: - for (int c = 0; c < Nc[g]; c++) // c: local coil index in group + for (int c = 0; c < fNc[g]; c++) // c: local coil index in group { - int i = Cin[g][c]; - Magsource2RemoteCoil(i, z0remG[g] - Z[g][c], roremG[g]); - for (int n = 2; n <= nmax; n++) - BremG[g][n] += Brem1[n]; + int i = fCin[g][c]; + Magsource2RemoteCoil(i, fz0remG[g] - fZ[g][c], froremG[g]); + for (int n = 2; n <= fnmax; n++) + fBremG[g][n] += fBrem1[n]; } - for (int n = 2; n <= nmax; n++) { - if (fabs(BremG[g][n]) < 1.e-30) - BremG[g][n] = 0.; - output << scientific << setw(9) << n << setw(26) << BremG[g][n] << endl; + for (int n = 2; n <= fnmax; n++) { + if (fabs(fBremG[g][n]) < 1.e-30) + fBremG[g][n] = 0.; + output << scientific << setw(9) << n << setw(26) << fBremG[g][n] << endl; } } output.close(); @@ -1061,28 +1061,28 @@ void MagfieldCoils::RemoteSourcepointGroup(int g) // origo at Z[g][0], local z axis parallel to group symmetry axis, in u vector direction. // Remote source point of group --> z0remG[g] (in group Z-coordinate system): double z0rem, zmin = 1.e9, zmax = -1.e9; - for (int c = 0; c < Nc[g]; c++) // c: local coil index in group + for (int c = 0; c < fNc[g]; c++) // c: local coil index in group { - int i = Cin[g][c]; - double L = coil[i][7]; - double zA = Z[g][c] - L / 2.; - double zB = Z[g][c] + L / 2.; // coil edges + int i = fCin[g][c]; + double L = fcoil[i][7]; + double zA = fZ[g][c] - L / 2.; + double zB = fZ[g][c] + L / 2.; // coil edges if (zA < zmin) zmin = zA; if (zB > zmax) zmax = zB; } - z0rem = z0remG[g] = (zmin + zmax) / 2.; // group center in group Z-coordinate system + z0rem = fz0remG[g] = (zmin + zmax) / 2.; // group center in group Z-coordinate system // Remote convergence radius of group --> roremG[g] : double rorem = 0.; - for (int c = 0; c < Nc[g]; c++) // c: local coil index in group + for (int c = 0; c < fNc[g]; c++) // c: local coil index in group { - int i = Cin[g][c]; - double roremc = Funrorem(i, z0rem - Z[g][c]); // remote source point relative to coil center + int i = fCin[g][c]; + double roremc = Funrorem(i, z0rem - fZ[g][c]); // remote source point relative to coil center if (roremc > rorem) rorem = roremc; } - roremG[g] = rorem; + froremG[g] = rorem; } @@ -1093,44 +1093,44 @@ void MagfieldCoils::RemoteSourcepointGroup(int g) // Magnetic charge remote source constant calculation for all coils void MagfieldCoils::MagsourceMagchargeCoils() -// This function computes the Vrem[i][n] (i=0,...,Ncoil-1; n=0,...,nmax) +// This function computes the Vrem[i][n] (i=0,...,fNcoil-1; n=0,...,nmax) // magnetic charge remote source constants for all coils. // The results are written into the data file dirname+objectname+magsource_remote.txt: -// Ncoil: number of coils, +// fNcoil: number of coils, // nmax: number of source coefficients for a fixed coil is nmax+1 (n=0,...,nmax), // i: coil index, n source constant index. { // Output to file dirname+objectname+magsource_remote.txt : double Rmin, Rmax, sigma; - string filename = dirname + objectname + "magsource_remote.txt"; + string filename = fdirname + fobjectname + "magsource_remote.txt"; ofstream output; output.precision(16); output.open(filename.c_str(), ios::app); // Dynamic memory allocations (Vrem): - Vrem = new double*[Ncoil]; - for (int i = 0; i < Ncoil; i++) - Vrem[i] = new double[nmax + 1]; + fVrem = new double*[fNcoil]; + for (int i = 0; i < fNcoil; i++) + fVrem[i] = new double[fnmax + 1]; // Coil index loop: - for (int i = 0; i < Ncoil; i++) { + for (int i = 0; i < fNcoil; i++) { // Coil i parameters: - Rmin = coil[i][8]; - Rmax = coil[i][9]; // inner and outer radius - sigma = coil[i][10]; // current density + Rmin = fcoil[i][8]; + Rmax = fcoil[i][9]; // inner and outer radius + sigma = fcoil[i][10]; // current density output << scientific << setw(7) << i << endl; // coil index - P[0] = 1.; + fP[0] = 1.; double ratio = Rmin / Rmax; double ration3 = ratio * ratio * ratio; - for (int n = 0; n <= nmax; n++) { + for (int n = 0; n <= fnmax; n++) { if (n > 0 && n % 2 == 0) // even positive n - P[n] = -(n - 1.) / double(n) * P[n - 2]; + fP[n] = -(n - 1.) / double(n) * fP[n - 2]; if (n % 2 == 1) - Vrem[i][n] = 0.; + fVrem[i][n] = 0.; else - Vrem[i][n] = sigma * P[n] * Rmax * Rmax / (2. * (n + 2.) * (n + 3.)) * (1. - ration3); + fVrem[i][n] = sigma * fP[n] * Rmax * Rmax / (2. * (n + 2.) * (n + 3.)) * (1. - ration3); if (n % 2 == 1) - output << scientific << setprecision(1) << setw(15) << n << setw(28) << Vrem[i][n] << endl; + output << scientific << setprecision(1) << setw(15) << n << setw(28) << fVrem[i][n] << endl; else - output << scientific << setprecision(16) << setw(15) << n << setw(28) << Vrem[i][n] << endl; + output << scientific << setprecision(16) << setw(15) << n << setw(28) << fVrem[i][n] << endl; if (ration3 < 1.e-20) ration3 = 0.; else @@ -1154,8 +1154,8 @@ double MagfieldCoils::Funrocen(int i, double z0) // (Zmin,Rmin) and (Zmax,Rmin). { double L, Zmin, Zmax, Rmin; - L = coil[i][7]; - Rmin = coil[i][8]; + L = fcoil[i][7]; + Rmin = fcoil[i][8]; Zmin = -L / 2.; Zmax = L / 2.; // coil endpoint Z values defined relative to the coil center if (z0 < 0.) @@ -1173,15 +1173,15 @@ double MagfieldCoils::FunrocenG(int g, double z0) // rocenmin = minimum distance of the axis point z0 from the windings of the group coils. { double rocenmin = 1.e9; - for (int c = 0; c < Nc[g]; c++) // c: local coil index in group + for (int c = 0; c < fNc[g]; c++) // c: local coil index in group { - int i = Cin[g][c]; // i: global coil index + int i = fCin[g][c]; // i: global coil index double L, Zmin, Zmax, Rmin, Z0; - L = coil[i][7]; - Rmin = coil[i][8]; + L = fcoil[i][7]; + Rmin = fcoil[i][8]; Zmin = -L / 2.; Zmax = L / 2.; // coil endpoint Z values defined relative to the coil center - Z0 = z0 - Z[g][c]; // source point relative to coil center + Z0 = z0 - fZ[g][c]; // source point relative to coil center double rocen; if (Z0 <= Zmin) rocen = sqrt((Z0 - Zmin) * (Z0 - Zmin) + Rmin * Rmin); @@ -1215,12 +1215,12 @@ void MagfieldCoils::Magsource2CentralCoil(int i, double z0, double Rocen) double x[1001], w[1001]; // nodes and weights double L, sigma, Zmin, Zmax, Rmin, Rmax; // Coil parameters: - L = coil[i][7]; + L = fcoil[i][7]; Zmin = -L / 2.; Zmax = L / 2.; // coil endpoints relative to coil center - Rmin = coil[i][8]; - Rmax = coil[i][9]; - sigma = coil[i][10]; // current density + Rmin = fcoil[i][8]; + Rmax = fcoil[i][9]; + sigma = fcoil[i][10]; // current density // Radial integration number N: int N; double ratio = (Rmax - Rmin) / Rmin; @@ -1232,10 +1232,10 @@ void MagfieldCoils::Magsource2CentralCoil(int i, double z0, double Rocen) if (N > 1000) N = 1000; // Initialization of Bcen1[n] and Pp[n]: - for (int n = 0; n <= nmax; n++) - Bcen1[n] = 0.; - Pp[0] = 0.; - Pp[1] = 1.; + for (int n = 0; n <= fnmax; n++) + fBcen1[n] = 0.; + fPp[0] = 0.; + fPp[1] = 1.; // C: double C = -0.5 * mu0 * sigma * Rocen; // Radial integration nodes and weights: @@ -1255,12 +1255,12 @@ void MagfieldCoils::Magsource2CentralCoil(int i, double z0, double Rocen) double b = a; double d = 1. / (ros * ros * ros); double D = C * R2 * sign[iz] * d * w[integ]; - Bcen1[0] += 0.5 * mu0 * sigma * sign[iz] * us * w[integ]; - Bcen1[1] += D; - for (int n = 2; n <= nmax; n++) { - double usPp = us * Pp[n - 1]; - Pp[n] = 2. * usPp - Pp[n - 2] + (usPp - Pp[n - 2]) * C1[n]; - Bcen1[n] += D * C0[n] * b * Pp[n]; + fBcen1[0] += 0.5 * mu0 * sigma * sign[iz] * us * w[integ]; + fBcen1[1] += D; + for (int n = 2; n <= fnmax; n++) { + double usPp = us * fPp[n - 1]; + fPp[n] = 2. * usPp - fPp[n - 2] + (usPp - fPp[n - 2]) * fC1[n]; + fBcen1[n] += D * fC0[n] * b * fPp[n]; b *= a; } } @@ -1276,7 +1276,7 @@ void MagfieldCoils::Magsource2CentralCoil(int i, double z0, double Rocen) void MagfieldCoils::MagsourceCentralCoils() // This function computes the central source constants for all coils. // The results are written into the data file dirname+objectname+magsource_central.txt: -// Ncoil: number of coils, +// fNcoil: number of coils, // nmax: number of source coefficients for a fixed coil and // source point is nmax+1 (n=0,...,nmax), // i: coil index, @@ -1287,35 +1287,35 @@ void MagfieldCoils::MagsourceCentralCoils() // rocen: central convergence radius for a fixed source point and coil. { // Output to file dirname+objectname+magsource_central.txt : - string filename = dirname + objectname + "magsource_central.txt"; + string filename = fdirname + fobjectname + "magsource_central.txt"; ofstream output; output.precision(16); output.open(filename.c_str()); // Dynamic memory allocations (Bcen1, Nsp, z0cen, rocen, Bcen): - Bcen1 = new double[nmax + 1]; - Nsp = new int[Ncoil]; - z0cen = new double*[Ncoil]; - rocen = new double*[Ncoil]; - Bcen = new double**[Ncoil]; + fBcen1 = new double[fnmax + 1]; + fNsp = new int[fNcoil]; + fz0cen = new double*[fNcoil]; + frocen = new double*[fNcoil]; + fBcen = new double**[fNcoil]; // Coil index loop: // output << "coils" << endl; - for (int i = 0; i < Ncoil; i++) { + for (int i = 0; i < fNcoil; i++) { CentralSourcepointsCoil(i); // central source point calculation - output << setw(9) << i << setw(9) << Nsp[i] << endl; + output << setw(9) << i << setw(9) << fNsp[i] << endl; // Central source constant calculation (source point loop): - Bcen[i] = new double*[Nsp[i]]; + fBcen[i] = new double*[fNsp[i]]; // Source point loop: // cout << "i, Nsp[i]= " << i << " " << Nsp[i] << endl; - for (int j = 0; j < Nsp[i]; j++) { - Magsource2CentralCoil(i, z0cen[i][j], rocen[i][j]); + for (int j = 0; j < fNsp[i]; j++) { + Magsource2CentralCoil(i, fz0cen[i][j], frocen[i][j]); // Output to file dirname+objectname+magsource_central.txt (j, z0, rocen, Bcen1[n]) - output << scientific << setw(9) << j << setw(26) << z0cen[i][j] << setw(26) << rocen[i][j] << endl; - Bcen[i][j] = new double[nmax + 1]; - for (int n = 0; n <= nmax; n++) { - if (fabs(Bcen1[n]) < 1.e-30) - Bcen1[n] = 0.; - output << scientific << setw(9) << n << setw(28) << Bcen1[n] << endl; - Bcen[i][j][n] = Bcen1[n]; + output << scientific << setw(9) << j << setw(26) << fz0cen[i][j] << setw(26) << frocen[i][j] << endl; + fBcen[i][j] = new double[fnmax + 1]; + for (int n = 0; n <= fnmax; n++) { + if (fabs(fBcen1[n]) < 1.e-30) + fBcen1[n] = 0.; + output << scientific << setw(9) << n << setw(28) << fBcen1[n] << endl; + fBcen[i][j][n] = fBcen1[n]; } } // end of source point loop } // end of coil index loop @@ -1332,7 +1332,7 @@ void MagfieldCoils::CentralSourcepointsCoil(int i) vector Z0cen; Z0cen.push_back(0.); double z0; - double roremC = rorem[i]; + double roremC = frorem[i]; double z0max = 2.5 * roremC; int nsp = 0; // number of central source points on positive z side int k = 0; @@ -1350,21 +1350,21 @@ void MagfieldCoils::CentralSourcepointsCoil(int i) nsp += 1; } while (z0 < z0max); // Number of central source points for coil i: - Nsp[i] = 2 * nsp + 1; + fNsp[i] = 2 * nsp + 1; // Central source points and conv. radii for coil i (source point loop): - z0cen[i] = new double[Nsp[i]]; - rocen[i] = new double[Nsp[i]]; + fz0cen[i] = new double[fNsp[i]]; + frocen[i] = new double[fNsp[i]]; // int j = nsp; - z0 = z0cen[i][j] = 0.; - rocen[i][j] = Funrocen(i, z0); + z0 = fz0cen[i][j] = 0.; + frocen[i][j] = Funrocen(i, z0); for (int k = 1; k <= nsp; k++) { j = nsp + k; - z0 = z0cen[i][j] = Z0cen[k]; - rocen[i][j] = Funrocen(i, z0); + z0 = fz0cen[i][j] = Z0cen[k]; + frocen[i][j] = Funrocen(i, z0); j = nsp - k; - z0 = z0cen[i][j] = -Z0cen[k]; - rocen[i][j] = Funrocen(i, z0); + z0 = fz0cen[i][j] = -Z0cen[k]; + frocen[i][j] = Funrocen(i, z0); } Z0cen.clear(); } @@ -1384,52 +1384,52 @@ void MagfieldCoils::MagsourceCentralGroups() // j: central source point index for a fixed group, { // Output to file dirname+objectname+magsource_central.txt : - string filename = dirname + objectname + "magsource_central.txt"; + string filename = fdirname + fobjectname + "magsource_central.txt"; ofstream output; output.precision(16); output.open(filename.c_str(), ios::app); // Dynamic memory allocations (NspG, z0cenG, rocenG, BcenG): - NspG = new int[Ng]; - z0cenG = new double*[Ng]; - rocenG = new double*[Ng]; - BcenG = new double**[Ng]; + fNspG = new int[fNg]; + fz0cenG = new double*[fNg]; + frocenG = new double*[fNg]; + fBcenG = new double**[fNg]; // Group index loop: // output << "groups" << endl; - for (int g = 0; g < Ng; g++) { - if (Nc[g] == 1) // no group calculation if the group has only 1 coil + for (int g = 0; g < fNg; g++) { + if (fNc[g] == 1) // no group calculation if the group has only 1 coil { - z0cenG[g] = new double[1]; - rocenG[g] = new double[1]; - NspG[g] = 1; - BcenG[g] = new double*[1]; - BcenG[g][0] = new double[nmax + 1]; + fz0cenG[g] = new double[1]; + frocenG[g] = new double[1]; + fNspG[g] = 1; + fBcenG[g] = new double*[1]; + fBcenG[g][0] = new double[fnmax + 1]; continue; } CentralSourcepointsGroup(g); // central source point calculation for group g // Number of central source points for group g: NspG[g] - output << setw(9) << g << setw(9) << NspG[g] << endl; + output << setw(9) << g << setw(9) << fNspG[g] << endl; // Central source constant calculation (source point loop): - BcenG[g] = new double*[NspG[g]]; - for (int j = 0; j < NspG[g]; j++) { + fBcenG[g] = new double*[fNspG[g]]; + for (int j = 0; j < fNspG[g]; j++) { // Output to file dirname+objectname+magsource_central.txt (j, z0, rocen) - output << scientific << setw(9) << j << setw(28) << z0cenG[g][j] << setw(28) << rocenG[g][j] << endl; + output << scientific << setw(9) << j << setw(28) << fz0cenG[g][j] << setw(28) << frocenG[g][j] << endl; // BcenG[g][j][n] initialization: - BcenG[g][j] = new double[nmax + 1]; - for (int n = 0; n <= nmax; n++) - BcenG[g][j][n] = 0.; + fBcenG[g][j] = new double[fnmax + 1]; + for (int n = 0; n <= fnmax; n++) + fBcenG[g][j][n] = 0.; // Local coil index loop: - for (int c = 0; c < Nc[g]; c++) // c: local coil index in group g + for (int c = 0; c < fNc[g]; c++) // c: local coil index in group g { - int i = Cin[g][c]; // global coil index - Magsource2CentralCoil(i, z0cenG[g][j] - Z[g][c], rocenG[g][j]); // source point relative to coil center - for (int n = 0; n <= nmax; n++) - BcenG[g][j][n] += Bcen1[n]; + int i = fCin[g][c]; // global coil index + Magsource2CentralCoil(i, fz0cenG[g][j] - fZ[g][c], frocenG[g][j]); // source point relative to coil center + for (int n = 0; n <= fnmax; n++) + fBcenG[g][j][n] += fBcen1[n]; } // Output to file dirname+objectname+magsource_central.txt of BcenG[g][j][n]: - for (int n = 0; n <= nmax; n++) { - if (fabs(BcenG[g][j][n]) < 1.e-30) - BcenG[g][j][n] = 0.; - output << scientific << setw(9) << n << setw(28) << BcenG[g][j][n] << endl; + for (int n = 0; n <= fnmax; n++) { + if (fabs(fBcenG[g][j][n]) < 1.e-30) + fBcenG[g][j][n] = 0.; + output << scientific << setw(9) << n << setw(28) << fBcenG[g][j][n] << endl; } } // end of source point j loop } // end of group index g loop @@ -1446,9 +1446,9 @@ void MagfieldCoils::CentralSourcepointsGroup(int g) { vector Z0cen; double z0; - double roremC = roremG[g]; - double z0min = z0remG[g] - 2. * roremC; - double z0max = z0remG[g] + 2. * roremC; + double roremC = froremG[g]; + double z0min = fz0remG[g] - 2. * roremC; + double z0max = fz0remG[g] + 2. * roremC; Z0cen.push_back(z0min); // j=0 int nsp = 1; // number of central source points in group g int j = 0; // --> z0min @@ -1466,14 +1466,14 @@ void MagfieldCoils::CentralSourcepointsGroup(int g) nsp += 1; } while (z0 < z0max); // Number of central source points for group g: - NspG[g] = nsp; + fNspG[g] = nsp; // Central source points and conv. radii for group g (source point loop): - z0cenG[g] = new double[nsp]; - rocenG[g] = new double[nsp]; + fz0cenG[g] = new double[nsp]; + frocenG[g] = new double[nsp]; // for (int j = 0; j < nsp; j++) { - z0 = z0cenG[g][j] = Z0cen[j]; - rocenG[g][j] = FunrocenG(g, z0); + z0 = fz0cenG[g][j] = Z0cen[j]; + frocenG[g][j] = FunrocenG(g, z0); } } @@ -1506,11 +1506,11 @@ bool MagfieldCoils::Magfield2Remote(bool bcoil, int ig, double z, double r, doub double z0, Rorem; // source point and conv. radius if (bcoil == true) { z0 = 0.; - Rorem = rorem[ig]; + Rorem = frorem[ig]; } else { - z0 = z0remG[ig]; - Rorem = roremG[ig]; + z0 = fz0remG[ig]; + Rorem = froremG[ig]; } double delz = z - z0; double ro = sqrt(r * r + delz * delz); @@ -1521,23 +1521,23 @@ bool MagfieldCoils::Magfield2Remote(bool bcoil, int ig, double z, double r, doub rc = Rorem / ro; // convergence ratio double rcn = rc * rc; // If rc>rclimit: the Legendre polynomial series is too slow, or not convergent (for rc>1) !!! - if (rc > rclimit) { + if (rc > frclimit) { return false; } // First 2 terms of Legendre polynomial P and its derivative Pp (P-primed) - P[0] = 1.; - P[1] = u; - Pp[0] = 0.; - Pp[1] = 1.; + fP[0] = 1.; + fP[1] = u; + fPp[0] = 0.; + fPp[1] = 1.; // bool even = false; // true: only the even n terms are used if (bcoil == false) // group calculation { double sum = 0., sumodd = 0.; for (int n = 2; n <= 12; n++) - sum += fabs(BremG[ig][n]); + sum += fabs(fBremG[ig][n]); for (int n = 3; n <= 12; n += 2) - sumodd += fabs(BremG[ig][n]); + sumodd += fabs(fBremG[ig][n]); if (sumodd < sum * 1.e-8) even = true; } @@ -1552,64 +1552,64 @@ bool MagfieldCoils::Magfield2Remote(bool bcoil, int ig, double z, double r, doub // Brem[ig][n]=0 for odd n) for (int n = 2; n <= 4; n++) { if (bcoil == true) - brem = Brem[ig][n]; + brem = fBrem[ig][n]; else - brem = BremG[ig][n]; + brem = fBremG[ig][n]; // P[n]=c1[n]*u*P[n-1]-c2[n]*P[n-2]; // Pp[n]=c3[n]*u*Pp[n-1]-c4[n]*Pp[n-2]; - double uPp = u * Pp[n - 1]; - double uP = u * P[n - 1]; - P[n] = 2. * uP - P[n - 2] - (uP - P[n - 2]) * C0[n]; - Pp[n] = 2. * uPp - Pp[n - 2] + (uPp - Pp[n - 2]) * C1[n]; + double uPp = u * fPp[n - 1]; + double uP = u * fP[n - 1]; + fP[n] = 2. * uP - fP[n - 2] - (uP - fP[n - 2]) * fC0[n]; + fPp[n] = 2. * uPp - fPp[n - 2] + (uPp - fPp[n - 2]) * fC1[n]; rcn = rcn * rc; - Bzplus[n] = brem * rcn * P[n]; - Brplus[n] = sr * brem * c5[n] * rcn * Pp[n]; - Bz += Bzplus[n]; - Br += Brplus[n]; + fBzplus[n] = brem * rcn * fP[n]; + fBrplus[n] = sr * brem * fc5[n] * rcn * fPp[n]; + Bz += fBzplus[n]; + Br += fBrplus[n]; } double u2 = u * u; double rc2 = rc * rc; - for (int n = 6; n <= nmax - 1; n += 2) { + for (int n = 6; n <= fnmax - 1; n += 2) { if (bcoil == true) - brem = Brem[ig][n]; + brem = fBrem[ig][n]; else - brem = BremG[ig][n]; + brem = fBremG[ig][n]; nlast = n; - P[n] = (c7[n] * u2 - c8[n]) * P[n - 2] - c9[n] * P[n - 4]; - Pp[n] = (c10[n] * u2 - c11[n]) * Pp[n - 2] - c12[n] * Pp[n - 4]; + fP[n] = (fc7[n] * u2 - fc8[n]) * fP[n - 2] - fc9[n] * fP[n - 4]; + fPp[n] = (fc10[n] * u2 - fc11[n]) * fPp[n - 2] - fc12[n] * fPp[n - 4]; rcn = rcn * rc2; - Bzplus[n] = brem * rcn * P[n]; - Brplus[n] = sr * brem * c5[n] * rcn * Pp[n]; - Bz += Bzplus[n]; - Br += Brplus[n]; + fBzplus[n] = brem * rcn * fP[n]; + fBrplus[n] = sr * brem * fc5[n] * rcn * fPp[n]; + Bz += fBzplus[n]; + Br += fBrplus[n]; double Beps = 1.e-15 * (fabs(Bz) + fabs(Br)); - double Bdelta = fabs(Bzplus[n]) + fabs(Brplus[n]) + fabs(Bzplus[n - 2]) + fabs(Brplus[n - 2]); + double Bdelta = fabs(fBzplus[n]) + fabs(fBrplus[n]) + fabs(fBzplus[n - 2]) + fabs(fBrplus[n - 2]); if (Bdelta < Beps || Bdelta < 1.e-20) break; } } else // even + odd group calculation: both even and odd terms are needed in the series { - for (int n = 2; n <= nmax - 1; n++) { - brem = BremG[ig][n]; + for (int n = 2; n <= fnmax - 1; n++) { + brem = fBremG[ig][n]; nlast = n; - P[n] = c1[n] * u * P[n - 1] - c2[n] * P[n - 2]; - Pp[n] = c3[n] * u * Pp[n - 1] - c4[n] * Pp[n - 2]; + fP[n] = fc1[n] * u * fP[n - 1] - fc2[n] * fP[n - 2]; + fPp[n] = fc3[n] * u * fPp[n - 1] - fc4[n] * fPp[n - 2]; rcn = rcn * rc; - Bzplus[n] = brem * rcn * P[n]; - Brplus[n] = sr * brem * c5[n] * rcn * Pp[n]; - Bz += Bzplus[n]; - Br += Brplus[n]; + fBzplus[n] = brem * rcn * fP[n]; + fBrplus[n] = sr * brem * fc5[n] * rcn * fPp[n]; + Bz += fBzplus[n]; + Br += fBrplus[n]; if (n > 5) { double Beps = 1.e-15 * (fabs(Bz) + fabs(Br)); - double Bdelta = fabs(Bzplus[n]) + fabs(Brplus[n]) + fabs(Bzplus[n - 1]) + fabs(Brplus[n - 1]) + - fabs(Bzplus[n - 2]) + fabs(Brplus[n - 2]) + fabs(Bzplus[n - 3]) + fabs(Brplus[n - 3]); + double Bdelta = fabs(fBzplus[n]) + fabs(fBrplus[n]) + fabs(fBzplus[n - 1]) + fabs(fBrplus[n - 1]) + + fabs(fBzplus[n - 2]) + fabs(fBrplus[n - 2]) + fabs(fBzplus[n - 3]) + fabs(fBrplus[n - 3]); if (Bdelta < Beps || Bdelta < 1.e-20) break; } } } - if (nlast >= nmax - 2) + if (nlast >= fnmax - 2) return false; else return true; @@ -1642,10 +1642,10 @@ bool MagfieldCoils::Magfield2Magcharge(int i, double z, double r, double& Bz, do const double mu0 = 4. * M_PI * 1.e-7; double Hzmax, Hrmax, Hzmin, Hrmin, Hz, Hr, Mz, rcmin, rcmax; double L, Zmin, Zmax, Rmin, Rmax, sigma; - L = coil[i][7]; - Rmin = coil[i][8]; - Rmax = coil[i][9]; // coil length, inner and outer radius - sigma = coil[i][10]; // current density + L = fcoil[i][7]; + Rmin = fcoil[i][8]; + Rmax = fcoil[i][9]; // coil length, inner and outer radius + sigma = fcoil[i][10]; // current density Zmin = -L / 2.; Zmax = L / 2.; // coil endpoints relative to coil center bool hfieldmax = Hfield(i, z - Zmax, r, Hzmax, Hrmax, rcmax); @@ -1685,7 +1685,7 @@ bool MagfieldCoils::Hfield(int i, double z, double r, double& Hz, double& Hr, do // rc>rclimit, and the return value is false // (in this case one should use another computation method !!!) double z0 = 0.; - double rorem = coil[i][9]; // Rmax + double rorem = fcoil[i][9]; // Rmax double delz = z - z0; double ro = sqrt(r * r + delz * delz); if (ro < 1.e-9) // field point is very close to the source point --> rc>1 ! @@ -1695,18 +1695,18 @@ bool MagfieldCoils::Hfield(int i, double z, double r, double& Hz, double& Hr, do rc = rorem / ro; // convergence ratio double rcn1 = rc * rc; // If rc>rclimit: the Legendre polynomial series is too slow, or not convergent (for rc>1) !!! - if (rc > rclimit) { + if (rc > frclimit) { return false; } // First 2 terms of Legendre polynomial P and its derivative Pp (P-primed) - P[0] = 1.; - P[1] = u; - Pp[0] = 0.; - Pp[1] = 1.; + fP[0] = 1.; + fP[1] = u; + fPp[0] = 0.; + fPp[1] = 1.; // for (int n = 2; n <= 3; n++) { - P[n] = c1[n] * u * P[n - 1] - c2[n] * P[n - 2]; - Pp[n] = c3[n] * u * Pp[n - 1] - c4[n] * Pp[n - 2]; + fP[n] = fc1[n] * u * fP[n - 1] - fc2[n] * fP[n - 2]; + fPp[n] = fc3[n] * u * fPp[n - 1] - fc4[n] * fPp[n - 2]; } double u2 = u * u; double rc2 = rc * rc; @@ -1717,14 +1717,14 @@ bool MagfieldCoils::Hfield(int i, double z, double r, double& Hz, double& Hr, do // Only the odd-n terms are needed in the series. double Hzplus, Hrplus, Hzplusold, Hrplusold, Heps, Hdelta; int nlast = 0; - for (int n = 1; n <= nmax - 1; n += 2) { + for (int n = 1; n <= fnmax - 1; n += 2) { nlast = n; if (n >= 4) { - P[n] = (c7[n] * u2 - c8[n]) * P[n - 2] - c9[n] * P[n - 4]; - Pp[n] = (c10[n] * u2 - c11[n]) * Pp[n - 2] - c12[n] * Pp[n - 4]; + fP[n] = (fc7[n] * u2 - fc8[n]) * fP[n - 2] - fc9[n] * fP[n - 4]; + fPp[n] = (fc10[n] * u2 - fc11[n]) * fPp[n - 2] - fc12[n] * fPp[n - 4]; } - Hzplus = Cz * Vrem[i][n - 1] * n * rcn1 * P[n]; - Hrplus = Cr * Vrem[i][n - 1] * rcn1 * Pp[n]; + Hzplus = Cz * fVrem[i][n - 1] * n * rcn1 * fP[n]; + Hrplus = Cr * fVrem[i][n - 1] * rcn1 * fPp[n]; Hz += Hzplus; Hr += Hrplus; Heps = 1.e-14 * (fabs(Hz) + fabs(Hr)); @@ -1737,7 +1737,7 @@ bool MagfieldCoils::Hfield(int i, double z, double r, double& Hz, double& Hr, do Hzplusold = Hzplus; Hrplusold = Hrplus; } - if (nlast >= nmax - 2) { + if (nlast >= fnmax - 2) { rc = 1.; return false; } @@ -1776,21 +1776,21 @@ bool MagfieldCoils::Magfield2Central(bool bcoil, int ig, int j, double z, double const double mu0 = 4. * M_PI * 1.e-7; double z0, Rocen; // source point and conv. radius if (bcoil == true) { - z0 = z0cen[ig][j]; - Rocen = rocen[ig][j]; + z0 = fz0cen[ig][j]; + Rocen = frocen[ig][j]; } else { - z0 = z0cenG[ig][j]; - Rocen = rocenG[ig][j]; + z0 = fz0cenG[ig][j]; + Rocen = frocenG[ig][j]; } double delz = z - z0; double ro = sqrt(r * r + delz * delz); // If the field point is very close to the source point: if (ro < 1.e-14) { if (bcoil == true) - Bz = Bcen[ig][j][0]; + Bz = fBcen[ig][j][0]; else - Bz = BcenG[ig][j][0]; + Bz = fBcenG[ig][j][0]; Br = 0.; rc = 0.; return true; @@ -1800,22 +1800,22 @@ bool MagfieldCoils::Magfield2Central(bool bcoil, int ig, int j, double z, double rc = ro / Rocen; // convergence ratio double rcn = rc; // If rc>=rclimit: the Legendre polynomial series is too slow, or not convergent (for rc>1) !!! - if (rc >= rclimit) { + if (rc >= frclimit) { return false; } // First 2 terms of Legendre polynomial P and its derivative Pp (P-primed) - P[0] = 1.; - P[1] = u; - Pp[0] = 0.; - Pp[1] = 1.; + fP[0] = 1.; + fP[1] = u; + fPp[0] = 0.; + fPp[1] = 1.; // First 2 terms of the series: if (bcoil == true) { - Bz = Bcen[ig][j][0] + Bcen[ig][j][1] * rc * u; - Br = -sr * Bcen[ig][j][1] / 2. * rc; + Bz = fBcen[ig][j][0] + fBcen[ig][j][1] * rc * u; + Br = -sr * fBcen[ig][j][1] / 2. * rc; } else { - Bz = BcenG[ig][j][0] + BcenG[ig][j][1] * rc * u; - Br = -sr * BcenG[ig][j][1] / 2. * rc; + Bz = fBcenG[ig][j][0] + fBcenG[ig][j][1] * rc * u; + Br = -sr * fBcenG[ig][j][1] / 2. * rc; } int nlast = 0; if (rc < 1.e-10) @@ -1823,42 +1823,42 @@ bool MagfieldCoils::Magfield2Central(bool bcoil, int ig, int j, double z, double // // We start here the central series expansion: double bcen; - for (int n = 2; n <= nmax - 1; n++) { + for (int n = 2; n <= fnmax - 1; n++) { if (bcoil == true) - bcen = Bcen[ig][j][n]; + bcen = fBcen[ig][j][n]; else - bcen = BcenG[ig][j][n]; + bcen = fBcenG[ig][j][n]; nlast = n; rcn *= rc; - double uPp = u * Pp[n - 1]; - double uP = u * P[n - 1]; - P[n] = 2. * uP - P[n - 2] - (uP - P[n - 2]) * C0[n]; - Pp[n] = 2. * uPp - Pp[n - 2] + (uPp - Pp[n - 2]) * C1[n]; + double uPp = u * fPp[n - 1]; + double uP = u * fP[n - 1]; + fP[n] = 2. * uP - fP[n - 2] - (uP - fP[n - 2]) * fC0[n]; + fPp[n] = 2. * uPp - fPp[n - 2] + (uPp - fPp[n - 2]) * fC1[n]; // P[n]=c1[n]*u*P[n-1]-c2[n]*P[n-2]; // Pp[n]=c3[n]*u*Pp[n-1]-c4[n]*Pp[n-2]; - Bzplus[n] = bcen * rcn * P[n]; - Brplus[n] = -sr * bcen * c6[n] * rcn * Pp[n]; - Bz += Bzplus[n]; - Br += Brplus[n]; + fBzplus[n] = bcen * rcn * fP[n]; + fBrplus[n] = -sr * bcen * fc6[n] * rcn * fPp[n]; + Bz += fBzplus[n]; + Br += fBrplus[n]; if (n > 5) { double Beps = 1.e-15 * (fabs(Bz) + fabs(Br)); - double Bdelta = fabs(Bzplus[n]) + fabs(Brplus[n]) + fabs(Bzplus[n - 1]) + fabs(Brplus[n - 1]) + - fabs(Bzplus[n - 2]) + fabs(Brplus[n - 2]) + fabs(Bzplus[n - 3]) + fabs(Brplus[n - 3]); + double Bdelta = fabs(fBzplus[n]) + fabs(fBrplus[n]) + fabs(fBzplus[n - 1]) + fabs(fBrplus[n - 1]) + + fabs(fBzplus[n - 2]) + fabs(fBrplus[n - 2]) + fabs(fBzplus[n - 3]) + fabs(fBrplus[n - 3]); if (Bdelta < Beps || Bdelta < 1.e-20) break; } } - if (nlast >= nmax - 1) { + if (nlast >= fnmax - 1) { rc = 1.; return false; } label:; // Effective central correction for Bz: if (bcoil == true) { - double Lhalf = coil[ig][7] * 0.5; + double Lhalf = fcoil[ig][7] * 0.5; if (z > -Lhalf && z < Lhalf) { - double Rmin = coil[ig][8], Rmax = coil[ig][9]; - double sigma = coil[ig][10]; // current density + double Rmin = fcoil[ig][8], Rmax = fcoil[ig][9]; + double sigma = fcoil[ig][10]; // current density if (r > Rmax) Bz += -mu0 * sigma * (Rmax - Rmin); else if (r > Rmin && r < Rmax) @@ -1894,10 +1894,10 @@ bool MagfieldCoils::MagfieldCoil(int i, const double* P, double* B) bool magfieldcoil = true; double C[3], u[3]; for (int k = 0; k < 3; k++) - C[k] = coil[i][1 + k]; // coil center - u[0] = coil[i][11]; // coil direction unit vector component ux - u[1] = coil[i][12]; // coil direction unit vector component uy - u[2] = coil[i][13]; // coil direction unit vector component uz + C[k] = fcoil[i][1 + k]; // coil center + u[0] = fcoil[i][11]; // coil direction unit vector component ux + u[1] = fcoil[i][12]; // coil direction unit vector component uy + u[2] = fcoil[i][13]; // coil direction unit vector component uz // Local cylindrical z and r coordinates of the field point P (relative to coil center and direction): double Ploc[3], Pr[3], w[3]; for (int k = 0; k < 3; k++) @@ -1917,7 +1917,7 @@ bool MagfieldCoils::MagfieldCoil(int i, const double* P, double* B) // we try first the remote series expansion // (with remote source point at the coil center) double ro2 = r2 + z * z; - double rorem2 = rorem[i] * rorem[i]; + double rorem2 = frorem[i] * frorem[i]; double Bz, Br, rc; if (ro2 > rorem2 * 2.) { bool magfield = Magfield2Remote(true, i, z, r, Bz, Br, rc); @@ -1926,8 +1926,8 @@ bool MagfieldCoils::MagfieldCoil(int i, const double* P, double* B) } // Step 2. We try the central zonal method in the beginning of tracking, // by searching all central source points. - if (jlast[i] == -2) { - jlast[i] = -1; + if (fjlast[i] == -2) { + fjlast[i] = -1; sps = SourcePointSearching(true, i, z, r, 0, jbest, rcbest); if (sps == true) { bool magfield = Magfield2Central(true, i, jbest, z, r, Bz, Br, rc); @@ -1943,9 +1943,9 @@ bool MagfieldCoils::MagfieldCoil(int i, const double* P, double* B) // has been already used (jlast[i]>-1), we search the central source // point j close to the old source point jlast[i]. sps = false; - if (jlast[i] > -1) + if (fjlast[i] > -1) sps = SourcePointSearching(true, i, z, r, 1, jbest, rcbest); - if (jlast[i] > -1 && sps == false) + if (fjlast[i] > -1 && sps == false) sps = SourcePointSearching(true, i, z, r, 2, jbest, rcbest); if (sps == true) { bool magfield = Magfield2Central(true, i, jbest, z, r, Bz, Br, rc); @@ -1964,8 +1964,8 @@ bool MagfieldCoils::MagfieldCoil(int i, const double* P, double* B) // --------------------------- // Step 5. We try now the magnetic charge method { - double L = coil[i][7]; - double Rmax = coil[i][9]; // coil length and outer radius + double L = fcoil[i][7]; + double Rmax = fcoil[i][9]; // coil length and outer radius double Zmin = -L / 2.; double Zmax = L / 2.; // coil endpoints relative to coil center double romax2 = r2 + (z - Zmax) * (z - Zmax); @@ -2023,8 +2023,8 @@ bool MagfieldCoils::MagfieldGroup(int g, const double* P, double* B) // The return value is true, if only zonal harmonic expansions are used for the magnetic field calc.; // if elliptic integral calc. is also used: the return value is false. // If the group has only 1 coil --> coil calculation: - if (Nc[g] == 1) - return MagfieldCoil(Cin[g][0], P, B); + if (fNc[g] == 1) + return MagfieldCoil(fCin[g][0], P, B); // // We start now the magnetic field calculation. // First we define the local coordinate system of the group, and the @@ -2034,9 +2034,9 @@ bool MagfieldCoils::MagfieldGroup(int g, const double* P, double* B) bool magfieldgroup = true; double C[3], u[3], w[3]; for (int k = 0; k < 3; k++) - C[k] = Line[g][k]; // group origo + C[k] = fLine[g][k]; // group origo for (int k = 0; k < 3; k++) - u[k] = Line[g][3 + k]; // group direction unit vector + u[k] = fLine[g][3 + k]; // group direction unit vector // Local cylindrical z and r coordinates of the field point P (relative to group origo and direction): double Ploc[3], Pr[3]; for (int k = 0; k < 3; k++) @@ -2054,8 +2054,8 @@ bool MagfieldCoils::MagfieldGroup(int g, const double* P, double* B) // --------------------------- // Step 1. We try first the remote series expansion // (with remote source point at the group center) - double ro2 = r2 + (z - z0remG[g]) * (z - z0remG[g]); - double rorem2 = roremG[g] * roremG[g]; + double ro2 = r2 + (z - fz0remG[g]) * (z - fz0remG[g]); + double rorem2 = froremG[g] * froremG[g]; double Bz, Br, rc; if (ro2 > rorem2 * 2.) { bool magfield = Magfield2Remote(false, g, z, r, Bz, Br, rc); @@ -2065,8 +2065,8 @@ bool MagfieldCoils::MagfieldGroup(int g, const double* P, double* B) // --------------------------- // Step 2. We try the central zonal method in the beginning of tracking, // by searching all central source points. - if (jlastG[g] == -2) { - jlastG[g] = -1; + if (fjlastG[g] == -2) { + fjlastG[g] = -1; sps = SourcePointSearching(false, g, z, r, 0, jbest, rcbest); if (sps == true) { bool magfield = Magfield2Central(false, g, jbest, z, r, Bz, Br, rc); @@ -2080,9 +2080,9 @@ bool MagfieldCoils::MagfieldGroup(int g, const double* P, double* B) // has been already used (jlastG[g]>-1), we search the central source // point j close to the old source point jlastG[g]. sps = false; - if (jlastG[g] > -1) + if (fjlastG[g] > -1) sps = SourcePointSearching(false, g, z, r, 1, jbest, rcbest); - if (jlastG[g] > -1 && sps == false) + if (fjlastG[g] > -1 && sps == false) sps = SourcePointSearching(false, g, z, r, 2, jbest, rcbest); if (sps == true) { bool magfield = Magfield2Central(false, g, jbest, z, r, Bz, Br, rc); @@ -2112,8 +2112,8 @@ bool MagfieldCoils::MagfieldGroup(int g, const double* P, double* B) B[0] = B[1] = B[2] = 0.; double Bcoil[3]; bool magfieldcoil; - for (int c = 0; c < Nc[g]; c++) { - int i = Cin[g][c]; + for (int c = 0; c < fNc[g]; c++) { + int i = fCin[g][c]; magfieldcoil = MagfieldCoil(i, P, Bcoil); for (int k = 0; k < 3; k++) B[k] += Bcoil[k]; @@ -2162,12 +2162,12 @@ bool MagfieldCoils::SourcePointSearching(bool coil, int ig, double z, double r, int i = ig; if (type == 0) // all central source points of the coil are searched { - rcmin2 = rclimit * rclimit; + rcmin2 = frclimit * frclimit; jbest = -1; - for (int j = 0; j < Nsp[i]; j++) { - delz = z - z0cen[i][j]; + for (int j = 0; j < fNsp[i]; j++) { + delz = z - fz0cen[i][j]; ro2 = r2 + delz * delz; - rocen2 = rocen[i][j] * rocen[i][j]; + rocen2 = frocen[i][j] * frocen[i][j]; if (ro2 < rocen2 * rcmin2) { rcmin2 = ro2 / rocen2; jbest = j; @@ -2175,7 +2175,7 @@ bool MagfieldCoils::SourcePointSearching(bool coil, int ig, double z, double r, } if (jbest > -1) { rcbest = sqrt(rcmin2); - jlast[i] = jbest; + fjlast[i] = jbest; return true; } else @@ -2183,25 +2183,25 @@ bool MagfieldCoils::SourcePointSearching(bool coil, int ig, double z, double r, } // end of coil=true, type=0 case else // only the 3 or the 11 source points near the last source point jlast are searched { - if (jlast[i] < 0) + if (fjlast[i] < 0) return false; int delj; if (type == 1) delj = 1; else delj = 5; - int jmin = jlast[i] - delj; + int jmin = fjlast[i] - delj; if (jmin < 0) jmin = 0; - int jmax = jlast[i] + delj; - if (jmax > Nsp[i] - 1) - jmax = Nsp[i] - 1; - rcmin2 = rclimit * rclimit; + int jmax = fjlast[i] + delj; + if (jmax > fNsp[i] - 1) + jmax = fNsp[i] - 1; + rcmin2 = frclimit * frclimit; jbest = -1; for (int j = jmin; j <= jmax; j++) { - delz = z - z0cen[i][j]; + delz = z - fz0cen[i][j]; ro2 = r2 + delz * delz; - rocen2 = rocen[i][j] * rocen[i][j]; + rocen2 = frocen[i][j] * frocen[i][j]; if (ro2 < rocen2 * rcmin2) { rcmin2 = ro2 / rocen2; jbest = j; @@ -2209,7 +2209,7 @@ bool MagfieldCoils::SourcePointSearching(bool coil, int ig, double z, double r, } if (jbest > -1) { rcbest = sqrt(rcmin2); - jlast[i] = jbest; + fjlast[i] = jbest; return true; } else @@ -2221,12 +2221,12 @@ bool MagfieldCoils::SourcePointSearching(bool coil, int ig, double z, double r, int g = ig; if (type == 0) // all central source points of the group are searched { - rcmin2 = rclimit * rclimit; + rcmin2 = frclimit * frclimit; jbest = -1; - for (int j = 0; j < NspG[g]; j++) { - delz = z - z0cenG[g][j]; + for (int j = 0; j < fNspG[g]; j++) { + delz = z - fz0cenG[g][j]; ro2 = r2 + delz * delz; - rocen2 = rocenG[g][j] * rocenG[g][j]; + rocen2 = frocenG[g][j] * frocenG[g][j]; if (ro2 < rocen2 * rcmin2) { rcmin2 = ro2 / rocen2; jbest = j; @@ -2234,7 +2234,7 @@ bool MagfieldCoils::SourcePointSearching(bool coil, int ig, double z, double r, } if (jbest > -1) { rcbest = sqrt(rcmin2); - jlastG[g] = jbest; + fjlastG[g] = jbest; return true; } else @@ -2242,25 +2242,25 @@ bool MagfieldCoils::SourcePointSearching(bool coil, int ig, double z, double r, } // end of coil=false, type=0 case else // only the 3 or the 11 source points near the last source point jlast are searched { - if (jlastG[g] < 0) + if (fjlastG[g] < 0) return false; int delj; if (type == 1) delj = 1; else delj = 5; - int jmin = jlastG[g] - delj; + int jmin = fjlastG[g] - delj; if (jmin < 0) jmin = 0; - int jmax = jlastG[g] + delj; - if (jmax > NspG[g] - 1) - jmax = NspG[g] - 1; - rcmin2 = rclimit * rclimit; + int jmax = fjlastG[g] + delj; + if (jmax > fNspG[g] - 1) + jmax = fNspG[g] - 1; + rcmin2 = frclimit * frclimit; jbest = -1; for (int j = jmin; j <= jmax; j++) { - delz = z - z0cenG[g][j]; + delz = z - fz0cenG[g][j]; ro2 = r2 + delz * delz; - rocen2 = rocenG[g][j] * rocenG[g][j]; + rocen2 = frocenG[g][j] * frocenG[g][j]; if (ro2 < rocen2 * rcmin2) { rcmin2 = ro2 / rocen2; jbest = j; @@ -2268,7 +2268,7 @@ bool MagfieldCoils::SourcePointSearching(bool coil, int ig, double z, double r, } if (jbest > -1) { rcbest = sqrt(rcmin2); - jlastG[g] = jbest; + fjlastG[g] = jbest; return true; } else @@ -2659,7 +2659,7 @@ void MagfieldCoils::GaussLegendreIntegration(int& N, double a, double b, double* The C++ code has 2 constructors: A, The first constructor has 6 input parameters; the user has to employ this one in the beginning of the calculation. This reads the coil parameters (by function CoilRead) from the file defined by the string coilfilename; these parameters are stored - in the array coil[Ncoil][14], where Ncoil denotes the number of coils. Function CoilRead + in the array coil[fNcoil][14], where fNcoil denotes the number of coils. Function CoilRead contains explanations about the coil input file structure. Then the coil parameters are tested, and symmetry groups are computed; the coil and symmetry group parameters are written into the From eee1f7965a64e40ab5b5f17921b9431599ac0dbe Mon Sep 17 00:00:00 2001 From: 2xB <31772910+2xB@users.noreply.github.com> Date: Wed, 8 May 2024 14:58:15 +0200 Subject: [PATCH 08/15] KEMField: Throw error in case of invalid nmax Not allowing nmax <= 1 and recommending nmax=500 was suggested by Ferenc. Fix previous commit Fix previous commits --- KEMField/Source/ExternalFields/MagfieldCoils/CMakeLists.txt | 1 + .../Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx | 4 ++++ 2 files changed, 5 insertions(+) diff --git a/KEMField/Source/ExternalFields/MagfieldCoils/CMakeLists.txt b/KEMField/Source/ExternalFields/MagfieldCoils/CMakeLists.txt index a492cd90e..22b8cee53 100644 --- a/KEMField/Source/ExternalFields/MagfieldCoils/CMakeLists.txt +++ b/KEMField/Source/ExternalFields/MagfieldCoils/CMakeLists.txt @@ -48,6 +48,7 @@ target_link_libraries (KEMMagfieldCoils # KEMCore # KEMMath # KEMFieldSolverCore + KEMFieldExceptions ) kasper_install_headers (${MAGFIELDCOILS_HEADERFILES}) diff --git a/KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx b/KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx index 304c76aea..e9052b21c 100755 --- a/KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx +++ b/KEMField/Source/ExternalFields/MagfieldCoils/src/MagfieldCoils.cxx @@ -1,4 +1,5 @@ #include "MagfieldCoils.h" +#include "KEMSimpleException.hh" MagfieldCoils::MagfieldCoils(string inputdirname, string inputobjectname, string inputcoilfilename, int inputNelliptic, @@ -14,6 +15,9 @@ MagfieldCoils::MagfieldCoils(string inputdirname, string inputobjectname, string fcoilfilename = inputcoilfilename; // name of coil input file fNelliptic = inputNelliptic; // radial numerical integration parameter for magnetic field calc. with elliptic integrals + if(inputnmax <= 1) { + throw KEMField::KEMSimpleException("MagfieldCoils: nmax can't be <= 1. A recommended value for nmax is 500.\nThis error was added in May 2024, before results in this case could not be trusted."); + } fnmax = inputnmax; // number of source coefficients for fixed source point: nmax+1 fepstol = inputepstol; // small tolerance parameter needed for the symmetry group definitions // Recommended values: inputNelliptic=32, inputnmax=500, inputepstol=1.e-8. From 805765f69aafb5227470b05448d980138df9ece5 Mon Sep 17 00:00:00 2001 From: 2xB <31772910+2xB@users.noreply.github.com> Date: Sun, 10 Dec 2023 14:45:29 +0100 Subject: [PATCH 09/15] Kassiopeia: Remove KSRoot::fToolbox Minor improvement - toolbox doesn't need to be a member of KSRoot --- Kassiopeia/Simulation/Include/KSRoot.h | 2 - Kassiopeia/Simulation/Source/KSRoot.cxx | 84 ++++++++++++------------- 2 files changed, 42 insertions(+), 44 deletions(-) diff --git a/Kassiopeia/Simulation/Include/KSRoot.h b/Kassiopeia/Simulation/Include/KSRoot.h index 3d417f15b..e6aaa5e58 100644 --- a/Kassiopeia/Simulation/Include/KSRoot.h +++ b/Kassiopeia/Simulation/Include/KSRoot.h @@ -56,8 +56,6 @@ class KSRoot : public KSComponentTemplate static void SignalHandler(int aSignal); private: - katrin::KToolbox& fToolbox; - KSSimulation* fSimulation; KSRun* fRun; KSEvent* fEvent; diff --git a/Kassiopeia/Simulation/Source/KSRoot.cxx b/Kassiopeia/Simulation/Source/KSRoot.cxx index 98167b894..e88eb1cef 100644 --- a/Kassiopeia/Simulation/Source/KSRoot.cxx +++ b/Kassiopeia/Simulation/Source/KSRoot.cxx @@ -50,7 +50,6 @@ bool KSRoot::fStopEventSignal = false; bool KSRoot::fStopTrackSignal = false; KSRoot::KSRoot() : - fToolbox(KToolbox::GetInstance()), fSimulation(nullptr), fRun(nullptr), fEvent(nullptr), @@ -77,126 +76,128 @@ KSRoot::KSRoot() : fStepIndex(0), fTotalExecTime(0) { + KToolbox& toolbox = KToolbox::GetInstance(); + fOnce = false; fRestartNavigation = true; this->SetName("root"); - if (fToolbox.Get("run") != nullptr) { + if (toolbox.Get("run") != nullptr) { mainmsg(eWarning) << "New Kassiopeia instance will re-use already existing root objects." << eom; - fRun = fToolbox.Get("run"); - fEvent = fToolbox.Get("event"); - fTrack = fToolbox.Get("track"); - fStep = fToolbox.Get("step"); - - fRootMagneticField = fToolbox.Get("root_magnetic_field"); - fRootElectricField = fToolbox.Get("root_electric_field"); - fRootSpace = fToolbox.Get("root_space"); - fRootGenerator = fToolbox.Get("root_generator"); - fRootTrajectory = fToolbox.Get("root_trajectory"); - fRootSpaceInteraction = fToolbox.Get("root_space_interaction"); - fRootSpaceNavigator = fToolbox.Get("root_space_navigator"); - fRootSurfaceInteraction = fToolbox.Get("root_surface_interaction"); - fRootSurfaceNavigator = fToolbox.Get("root_surface_navigator"); - fRootTerminator = fToolbox.Get("root_terminator"); - fRootWriter = fToolbox.Get("root_writer"); - fRootStepModifier = fToolbox.Get("root_step_modifier"); - fRootTrackModifier = fToolbox.Get("root_track_modifier"); - fRootEventModifier = fToolbox.Get("root_event_modifier"); - fRootRunModifier = fToolbox.Get("root_run_modifier"); + fRun = toolbox.Get("run"); + fEvent = toolbox.Get("event"); + fTrack = toolbox.Get("track"); + fStep = toolbox.Get("step"); + + fRootMagneticField = toolbox.Get("root_magnetic_field"); + fRootElectricField = toolbox.Get("root_electric_field"); + fRootSpace = toolbox.Get("root_space"); + fRootGenerator = toolbox.Get("root_generator"); + fRootTrajectory = toolbox.Get("root_trajectory"); + fRootSpaceInteraction = toolbox.Get("root_space_interaction"); + fRootSpaceNavigator = toolbox.Get("root_space_navigator"); + fRootSurfaceInteraction = toolbox.Get("root_surface_interaction"); + fRootSurfaceNavigator = toolbox.Get("root_surface_navigator"); + fRootTerminator = toolbox.Get("root_terminator"); + fRootWriter = toolbox.Get("root_writer"); + fRootStepModifier = toolbox.Get("root_step_modifier"); + fRootTrackModifier = toolbox.Get("root_track_modifier"); + fRootEventModifier = toolbox.Get("root_event_modifier"); + fRootRunModifier = toolbox.Get("root_run_modifier"); return; } fRun = new KSRun(); fRun->SetName("run"); - fToolbox.Add(fRun, "run"); + toolbox.Add(fRun, "run"); fEvent = new KSEvent(); fEvent->SetName("event"); - fToolbox.Add(fEvent); + toolbox.Add(fEvent); fTrack = new KSTrack(); fTrack->SetName("track"); - fToolbox.Add(fTrack); + toolbox.Add(fTrack); fStep = new KSStep(); fStep->SetName("step"); - fToolbox.Add(fStep); + toolbox.Add(fStep); fRootMagneticField = new KSRootMagneticField(); fRootMagneticField->SetName("root_magnetic_field"); - fToolbox.Add(fRootMagneticField); + toolbox.Add(fRootMagneticField); fRootElectricField = new KSRootElectricField(); fRootElectricField->SetName("root_electric_field"); - fToolbox.Add(fRootElectricField); + toolbox.Add(fRootElectricField); fRootSpace = new KSRootSpace(); fRootSpace->SetName("root_space"); - fToolbox.Add(fRootSpace); + toolbox.Add(fRootSpace); fRootGenerator = new KSRootGenerator(); fRootGenerator->SetName("root_generator"); fRootGenerator->SetEvent(fEvent); - fToolbox.Add(fRootGenerator); + toolbox.Add(fRootGenerator); fRootTrajectory = new KSRootTrajectory(); fRootTrajectory->SetName("root_trajectory"); fRootTrajectory->SetStep(fStep); - fToolbox.Add(fRootTrajectory); + toolbox.Add(fRootTrajectory); fRootSpaceInteraction = new KSRootSpaceInteraction(); fRootSpaceInteraction->SetName("root_space_interaction"); fRootSpaceInteraction->SetStep(fStep); fRootSpaceInteraction->SetTrajectory(fRootTrajectory); - fToolbox.Add(fRootSpaceInteraction); + toolbox.Add(fRootSpaceInteraction); fRootSpaceNavigator = new KSRootSpaceNavigator(); fRootSpaceNavigator->SetName("root_space_navigator"); fRootSpaceNavigator->SetStep(fStep); fRootSpaceNavigator->SetTrajectory(fRootTrajectory); - fToolbox.Add(fRootSpaceNavigator); + toolbox.Add(fRootSpaceNavigator); fRootSurfaceInteraction = new KSRootSurfaceInteraction(); fRootSurfaceInteraction->SetName("root_surface_interaction"); fRootSurfaceInteraction->SetStep(fStep); - fToolbox.Add(fRootSurfaceInteraction); + toolbox.Add(fRootSurfaceInteraction); fRootSurfaceNavigator = new KSRootSurfaceNavigator(); fRootSurfaceNavigator->SetName("root_surface_navigator"); fRootSurfaceNavigator->SetStep(fStep); - fToolbox.Add(fRootSurfaceNavigator); + toolbox.Add(fRootSurfaceNavigator); fRootTerminator = new KSRootTerminator(); fRootTerminator->SetName("root_terminator"); fRootTerminator->SetStep(fStep); - fToolbox.Add(fRootTerminator); + toolbox.Add(fRootTerminator); fRootWriter = new KSRootWriter(); fRootWriter->SetName("root_writer"); - fToolbox.Add(fRootWriter); + toolbox.Add(fRootWriter); fRootStepModifier = new KSRootStepModifier(); fRootStepModifier->SetName("root_step_modifier"); fRootStepModifier->SetStep(fStep); - fToolbox.Add(fRootStepModifier); + toolbox.Add(fRootStepModifier); fRootTrackModifier = new KSRootTrackModifier(); fRootTrackModifier->SetName("root_track_modifier"); fRootTrackModifier->SetTrack(fTrack); - fToolbox.Add(fRootTrackModifier); + toolbox.Add(fRootTrackModifier); fRootEventModifier = new KSRootEventModifier(); fRootEventModifier->SetName("root_event_modifier"); fRootEventModifier->SetEvent(fEvent); - fToolbox.Add(fRootEventModifier); + toolbox.Add(fRootEventModifier); fRootRunModifier = new KSRootRunModifier(); fRootRunModifier->SetName("root_run_modifier"); fRootRunModifier->SetRun(fRun); - fToolbox.Add(fRootRunModifier); + toolbox.Add(fRootRunModifier); // convert KEMField objects to Kassiopeia components for (auto& name : KToolbox::GetInstance().FindAll()) { @@ -214,7 +215,6 @@ KSRoot::KSRoot() : } KSRoot::KSRoot(const KSRoot& aCopy) : KSComponent(aCopy), - fToolbox(KToolbox::GetInstance()), fSimulation(nullptr), fRun(aCopy.fRun), fEvent(aCopy.fEvent), From 6b5f21542e7eeb39fe414b9c8088f5689f6ad5b8 Mon Sep 17 00:00:00 2001 From: 2xB <2xB@users.noreply.github.com> Date: Mon, 10 Jun 2024 20:56:56 +0200 Subject: [PATCH 10/15] Fix compiler warning about uninitialized member gcc does not think fHolder was initialized. This fixes that by using a unique_ptr that is initialized automatically. --- Kassiopeia/Objects/Include/KSObject.h | 9 +++------ Kassiopeia/Objects/Source/KSObject.cxx | 4 ++-- 2 files changed, 5 insertions(+), 8 deletions(-) diff --git a/Kassiopeia/Objects/Include/KSObject.h b/Kassiopeia/Objects/Include/KSObject.h index a4ba1a3ac..5f02dd2d0 100644 --- a/Kassiopeia/Objects/Include/KSObject.h +++ b/Kassiopeia/Objects/Include/KSObject.h @@ -3,6 +3,7 @@ #include "KSObjectsMessage.h" #include "KTagged.h" +#include namespace Kassiopeia { @@ -53,7 +54,7 @@ class KSObject : public katrin::KTagged XType* fObject; }; - mutable KSHolder* fHolder; + mutable std::unique_ptr fHolder; }; inline KSObject::KSHolder::KSHolder() = default; @@ -141,12 +142,8 @@ template<> inline const KSObject* KSObject::As() const template inline void KSObject::Set(XType* anObject) { - if (fHolder != nullptr) { - delete fHolder; - fHolder = nullptr; - } auto* tHolder = new KSHolderTemplate(anObject); - fHolder = tHolder; + fHolder = std::unique_ptr(tHolder); return; } diff --git a/Kassiopeia/Objects/Source/KSObject.cxx b/Kassiopeia/Objects/Source/KSObject.cxx index e33b0af8f..8660d414a 100644 --- a/Kassiopeia/Objects/Source/KSObject.cxx +++ b/Kassiopeia/Objects/Source/KSObject.cxx @@ -3,8 +3,8 @@ namespace Kassiopeia { -KSObject::KSObject() : KTagged(), fHolder(nullptr) {} -KSObject::KSObject(const KSObject& aCopy) : KTagged(aCopy), fHolder(nullptr) {} +KSObject::KSObject() : KTagged() {}; +KSObject::KSObject(const KSObject& aCopy) : KTagged(aCopy) {} KSObject::~KSObject() = default; } // namespace Kassiopeia From 81d617309a80ebebc5471019eb88c53366db0c11 Mon Sep 17 00:00:00 2001 From: 2xB <31772910+2xB@users.noreply.github.com> Date: Mon, 10 Jun 2024 21:32:34 +0200 Subject: [PATCH 11/15] CI: Fedora 37 => Fedora 40 --- .github/workflows/test.yml | 12 +++--------- Documentation/gh-pages/source/setup_manual.rst | 2 +- 2 files changed, 4 insertions(+), 10 deletions(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 5c1d7183b..9560ffe2f 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -8,7 +8,7 @@ on: pull_request: branches: [ "main" ] -# Keep the following in sync with Documentation/gh-pages/source/compiling.rst ! +# Keep the following in sync with Documentation/gh-pages/source/setup_manual.rst ! jobs: ubuntu_20_04: strategy: @@ -51,20 +51,14 @@ jobs: source root/bin/thisroot.sh && source install/bin/kasperenv.sh && UnitTestKasper shell: bash - fedora_37: + fedora_40: strategy: matrix: use_clang: [false] # FIXME add "true" after solving https://github.com/KATRIN-Experiment/Kassiopeia/issues/87 runs-on: ubuntu-latest - container: fedora:37 + container: fedora:40 steps: - uses: actions/checkout@v3 - - name: Run workaround for log4cxx bug (no longer needed with Fedora 39) - run: | - dnf update -y - dnf install -y --setopt=install_weak_deps=False dnf-plugins-core - dnf clean all - dnf copr enable thofmann/log4xx-1.x -y - name: Install dependencies run: | dnf install -y \ diff --git a/Documentation/gh-pages/source/setup_manual.rst b/Documentation/gh-pages/source/setup_manual.rst index b61169dad..b671c70d9 100644 --- a/Documentation/gh-pages/source/setup_manual.rst +++ b/Documentation/gh-pages/source/setup_manual.rst @@ -100,7 +100,7 @@ On a RedHat/Fedora Linux system, the packages can be installed through the comma boost-devel fftw-devel gsl-devel hdf5-devel libomp-devel liburing-devel libxml2-devel log4cxx-devel \ ocl-icd-devel openmpi-devel openssl-devel sqlite-devel vtk-devel zlib-devel -Tested on Fedora Linux 37. +Tested on Fedora Linux 40. Required dependencies ---------------------- From 8cd6b78a5ccc0652dc18420b05fa3a067b706a48 Mon Sep 17 00:00:00 2001 From: 2xB <31772910+2xB@users.noreply.github.com> Date: Mon, 10 Jun 2024 23:53:47 +0200 Subject: [PATCH 12/15] kasperenv.sh: Fix PATH pollution Previously, code that was intended to remove old KASPER entries from PATH variables was broken, this is fixed. The code didn't work since it was introduced in January 2020 due to using single quotes and therefore avoiding variable replacement. Since "/" is also a letter that is frequently used in paths, to avoid future problems the letter was replaced with ";", which `thisroot.sh` also uses for this cause. --- kasperenv.sh.in | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/kasperenv.sh.in b/kasperenv.sh.in index b77d87467..7ea995341 100755 --- a/kasperenv.sh.in +++ b/kasperenv.sh.in @@ -67,10 +67,20 @@ fi export KEMFIELD_CACHE mkdir -p ${KEMFIELD_CACHE} -export PATH=${KASPER_INSTALL_BIN}:$(echo $PATH | sed 's/${OLD_PATH}//g') -export LD_LIBRARY_PATH=${KASPER_INSTALL_LIB}:$(echo $LD_LIBRARY_PATH | sed 's/${OLD_LDLIBPATH}//g') -export PKG_CONFIG_PATH=${KASPER_INSTALL_LIB}/pkgconfig:$(echo $PKG_CONFIG_PATH | sed 's/${OLD_PKGCFGPATH}//g') -export PYTHONPATH=${KASPER_INSTALL_LIB}/python:$(echo $PYTHONPATH | sed 's/${OLD_PYTHONPATH}//g') -export CMAKE_PREFIX_PATH=${KASPER_INSTALL}:$(echo $CMAKE_PREFIX_PATH | sed 's/${OLD_CMAKE_PREF}//g') + +[[ -n $OLD_PATH ]] && export PATH=$(echo $PATH | sed "s;${OLD_PATH};;g") +export PATH=${KASPER_INSTALL_BIN}:$PATH + +[[ -n $OLD_LDLIBPATH ]] && export LD_LIBRARY_PATH=$(echo $LD_LIBRARY_PATH | sed "s;${OLD_LDLIBPATH};;g") +export LD_LIBRARY_PATH=${KASPER_INSTALL_LIB}:$LD_LIBRARY_PATH + +[[ -n $OLD_PKGCFGPATH ]] && export PKG_CONFIG_PATH=$(echo $PKG_CONFIG_PATH | sed "s;${OLD_PKGCFGPATH};;g") +export PKG_CONFIG_PATH=${KASPER_INSTALL_LIB}/pkgconfig:$PKG_CONFIG_PATH + +[[ -n $OLD_PYTHONPATH ]] && export PYTHONPATH=$(echo $PYTHONPATH | sed "s;${OLD_PYTHONPATH};;g") +export PYTHONPATH=${KASPER_INSTALL_LIB}/python:$PYTHONPATH + +[[ -n $OLD_CMAKE_PREF ]] && export CMAKE_PREFIX_PATH=$(echo $CMAKE_PREFIX_PATH | sed "s;${OLD_CMAKE_PREF};;g") +export CMAKE_PREFIX_PATH=${KASPER_INSTALL}:$CMAKE_PREFIX_PATH return 0 From 332bac557a7d6cead169f3fe7a4b9f7aa1add3b5 Mon Sep 17 00:00:00 2001 From: 2xB <31772910+2xB@users.noreply.github.com> Date: Mon, 10 Jun 2024 23:54:41 +0200 Subject: [PATCH 13/15] kasperenv.sh: Fix "Error in cling::[...]" Adding KASPER include path to `ROOT_INCLUDE_PATH` to avoid errors of type `Error in cling::AutoLoadingVisitor::InsertIntoAutoLoadingState:`. --- kasperenv.sh.in | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/kasperenv.sh.in b/kasperenv.sh.in index 7ea995341..903ce438d 100755 --- a/kasperenv.sh.in +++ b/kasperenv.sh.in @@ -14,6 +14,7 @@ if [ -n "${KASPERSYS}" ] OLD_PKGCFGPATH=${OLD_KASPERSYS}/lib/pkgconfig: OLD_PYTHONPATH=${OLD_KASPERSYS}/lib/python: OLD_CMAKE_PREF=${OLD_KASPERSYS}: + OLD_ROOT_INCLUDE_PATH=${OLD_KASPERSYS}/include: fi if [ -n "${KASPER_INSTALL}" ] @@ -24,6 +25,7 @@ if [ -n "${KASPER_INSTALL}" ] OLD_PKGCFGPATH=${OLD_KASPER_INSTALL}/lib/pkgconfig: OLD_PYTHONPATH=${OLD_KASPER_INSTALL}/lib/python: OLD_CMAKE_PREF=${OLD_KASPER_INSTALL}: + OLD_ROOT_INCLUDE_PATH=${OLD_KASPER_INSTALL}/include: fi if [ -n "${KASPER_SOURCE}" ] @@ -36,6 +38,7 @@ export KASPER_INSTALL=@KASPER_INSTALL_DIR@ KASPER_INSTALL_BIN=@BIN_INSTALL_DIR@ KASPER_INSTALL_LIB=@LIB_INSTALL_DIR@ KASPER_INSTALL_CACHE=@CACHE_INSTALL_DIR@ +KASPER_INSTALL_INCLUDE=@INCLUDE_INSTALL_DIR@ printf "\033[32;1m** KASPER source directory set to ${KASPER_SOURCE}\033[0m\n" printf "\033[32;1m** KASPER install directory set to ${KASPER_INSTALL}\033[0m\n" @@ -83,4 +86,7 @@ export PYTHONPATH=${KASPER_INSTALL_LIB}/python:$PYTHONPATH [[ -n $OLD_CMAKE_PREF ]] && export CMAKE_PREFIX_PATH=$(echo $CMAKE_PREFIX_PATH | sed "s;${OLD_CMAKE_PREF};;g") export CMAKE_PREFIX_PATH=${KASPER_INSTALL}:$CMAKE_PREFIX_PATH +[[ -n $OLD_ROOT_INCLUDE_PATH ]] && export ROOT_INCLUDE_PATH=$(echo $ROOT_INCLUDE_PATH | sed "s;${OLD_ROOT_INCLUDE_PATH};;g") +export ROOT_INCLUDE_PATH=${KASPER_INSTALL_INCLUDE}:$ROOT_INCLUDE_PATH + return 0 From 54be51f4f3d0e35cb56236731ce29cd29089057f Mon Sep 17 00:00:00 2001 From: 2xB <31772910+2xB@users.noreply.github.com> Date: Mon, 10 Jun 2024 23:58:10 +0200 Subject: [PATCH 14/15] Docker: Added diffutils & paraview to full image --- Docker/packages.full | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Docker/packages.full b/Docker/packages.full index 07a1543c6..b3f4884c3 100644 --- a/Docker/packages.full +++ b/Docker/packages.full @@ -3,6 +3,8 @@ vim zsh zsh-syntax-highlighting htop +diffutils +paraview root root-notebook From 350dbb5769ea262113fb125c35dba7ad0c0685ff Mon Sep 17 00:00:00 2001 From: 2xB <31772910+2xB@users.noreply.github.com> Date: Mon, 10 Jun 2024 23:59:22 +0200 Subject: [PATCH 15/15] Docker: Install uproot in full image --- Dockerfile | 1 + 1 file changed, 1 insertion(+) diff --git a/Dockerfile b/Dockerfile index 283dbbcc9..ee551fd03 100644 --- a/Dockerfile +++ b/Dockerfile @@ -136,6 +136,7 @@ RUN pip3 install --no-cache-dir jupyterlab \ && pip3 install --no-cache-dir jupyter-server-proxy \ && pip3 install --no-cache-dir jupyterhub \ && pip3 install --no-cache-dir ipympl \ + && pip3 install --no-cache-dir uproot \ && pip3 install --no-cache-dir iminuit # Ensure if LDAP is used on a JupyterHub, user names are correctly resolved