Garfield++ v1r0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Functions
Garfield::Numerics Namespace Reference

Functions

void Dfact (const int n, std::vector< std::vector< double > > &a, std::vector< int > &ir, int &ifail, double &det, int &jfail)
 
void Dfeqn (const int n, std::vector< std::vector< double > > &a, std::vector< int > &ir, std::vector< double > &b)
 
void Dfinv (const int n, std::vector< std::vector< double > > &a, std::vector< int > &ir)
 
void Deqinv (const int n, std::vector< std::vector< double > > &a, int &ifail, std::vector< double > &b)
 
void Cfact (const int n, std::vector< std::vector< std::complex< double > > > &a, std::vector< int > &ir, int &ifail, std::complex< double > &det, int &jfail)
 
void Cfinv (const int n, std::vector< std::vector< std::complex< double > > > &a, std::vector< int > &ir)
 
void Cinv (const int n, std::vector< std::vector< std::complex< double > > > &a, int &ifail)
 
double GaussKronrod15 (double(*f)(const double), const double a, const double b)
 
double BesselI0S (const double xx)
 
double BesselI1S (const double xx)
 
double BesselK0S (const double xx)
 
double BesselK0L (const double xx)
 
double BesselK1S (const double xx)
 
double BesselK1L (const double xx)
 
double Divdif (const std::vector< double > &f, const std::vector< double > &a, int nn, double x, int mm)
 
bool Boxin3 (std::vector< std::vector< std::vector< double > > > &value, std::vector< double > &xAxis, std::vector< double > &yAxis, std::vector< double > &zAxis, int nx, int ny, int nz, double xx, double yy, double zz, double &f, int iOrder)
 

Function Documentation

◆ BesselI0S()

double Garfield::Numerics::BesselI0S ( const double  xx)
inline

Definition at line 37 of file Numerics.hh.

37 {
38 return 1. + 3.5156229 * pow(xx / 3.75, 2) + 3.0899424 * pow(xx / 3.75, 4) +
39 1.2067492 * pow(xx / 3.75, 6) + 0.2659732 * pow(xx / 3.75, 8) +
40 0.0360768 * pow(xx / 3.75, 10) + 0.0045813 * pow(xx / 3.75, 12);
41}
pow
DoubleAc pow(const DoubleAc &f, double p)
Definition: DoubleAc.cpp:336

Referenced by BesselK0S().

◆ BesselI1S()

double Garfield::Numerics::BesselI1S ( const double  xx)
inline

Definition at line 43 of file Numerics.hh.

43 {
44 return xx *
45 (0.5 + 0.87890594 * pow(xx / 3.75, 2) +
46 0.51498869 * pow(xx / 3.75, 4) + 0.15084934 * pow(xx / 3.75, 6) +
47 0.02658733 * pow(xx / 3.75, 8) + 0.00301532 * pow(xx / 3.75, 10) +
48 0.00032411 * pow(xx / 3.75, 12));
49}

Referenced by BesselK1S().

◆ BesselK0L()

double Garfield::Numerics::BesselK0L ( const double  xx)
inline

Definition at line 58 of file Numerics.hh.

58 {
59 return (exp(-xx) / sqrt(xx)) *
60 (1.25331414 - 0.07832358 * (2. / xx) + 0.02189568 * pow(2. / xx, 2) -
61 0.01062446 * pow(2. / xx, 3) + 0.00587872 * pow(2. / xx, 4) -
62 0.00251540 * pow(2. / xx, 5) + 0.00053208 * pow(2. / xx, 6));
63}
sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:313
exp
DoubleAc exp(const DoubleAc &f)
Definition: DoubleAc.cpp:376

◆ BesselK0S()

double Garfield::Numerics::BesselK0S ( const double  xx)
inline

Definition at line 51 of file Numerics.hh.

51 {
52 return -log(xx / 2.) * BesselI0S(xx) - 0.57721566 +
53 0.42278420 * pow(xx / 2., 2) + 0.23069756 * pow(xx / 2., 4) +
54 0.03488590 * pow(xx / 2., 6) + 0.00262698 * pow(xx / 2., 8) +
55 0.00010750 * pow(xx / 2., 10) + 0.00000740 * pow(xx / 2., 12);
56}
Garfield::Numerics::BesselI0S
double BesselI0S(const double xx)
Definition: Numerics.hh:37

◆ BesselK1L()

double Garfield::Numerics::BesselK1L ( const double  xx)
inline

Definition at line 73 of file Numerics.hh.

73 {
74 return (exp(-xx) / sqrt(xx)) *
75 (1.25331414 + 0.23498619 * (2. / xx) - 0.03655620 * pow(2. / xx, 2) +
76 0.01504268 * pow(2. / xx, 3) - 0.00780353 * pow(2. / xx, 4) +
77 0.00325614 * pow(2. / xx, 5) - 0.00068245 * pow(2. / xx, 6));
78}

◆ BesselK1S()

double Garfield::Numerics::BesselK1S ( const double  xx)
inline

Definition at line 65 of file Numerics.hh.

65 {
66 return log(xx / 2.) * BesselI1S(xx) +
67 (1. / xx) *
68 (1. + 0.15443144 * pow(xx / 2., 2) - 0.67278579 * pow(xx / 2., 4) -
69 0.18156897 * pow(xx / 2., 6) - 0.01919402 * pow(xx / 2., 8) -
70 0.00110404 * pow(xx / 2., 10) - 0.00004686 * pow(xx / 2., 12));
71}
Garfield::Numerics::BesselI1S
double BesselI1S(const double xx)
Definition: Numerics.hh:43

◆ Boxin3()

bool Garfield::Numerics::Boxin3 ( std::vector< std::vector< std::vector< double > > > &  value,
std::vector< double > &  xAxis,
std::vector< double > &  yAxis,
std::vector< double > &  zAxis,
int  nx,
int  ny,
int  nz,
double  xx,
double  yy,
double  zz,
double &  f,
int  iOrder 
)

Definition at line 771 of file Numerics.cc.

774 {
775
776 // std::cout << nx << ", " << ny << ", " << nz << "\n";
777 //-----------------------------------------------------------------------
778 // BOXIN3 - interpolation of order 1 and 2 in an irregular rectangular
779 // 3-dimensional grid.
780 // (Last changed on 13/ 2/00.)
781 //-----------------------------------------------------------------------
782
783 int iX0 = 0, iX1 = 0;
784 int iY0 = 0, iY1 = 0;
785 int iZ0 = 0, iZ1 = 0;
786 double fX[4], fY[4], fZ[4];
787
788 // Ensure we are in the grid.
789 const double x = std::min(std::max(xx, std::min(xAxis[0], xAxis[nx - 1])),
790 std::max(xAxis[0], xAxis[nx - 1]));
791 const double y = std::min(std::max(yy, std::min(yAxis[0], yAxis[ny - 1])),
792 std::max(yAxis[0], yAxis[ny - 1]));
793 const double z = std::min(std::max(zz, std::min(zAxis[0], zAxis[nz - 1])),
794 std::max(zAxis[0], zAxis[nz - 1]));
795
796 // Make sure we have enough points.
797 if (iOrder < 0 || iOrder > 2 || nx < 1 || ny < 1 || nz < 1) {
798 std::cerr << "Boxin3:\n";
799 std::cerr << " Incorrect order or number of points.\n";
800 std::cerr << " No interpolation.\n";
801 f = 0.;
802 return false;
803 }
804
805 if (iOrder == 0 || nx == 1) {
806 // Zeroth order interpolation in x.
807 // Find the nearest node.
808 double dist = fabs(x - xAxis[0]);
809 int iNode = 0;
810 for (int i = 1; i < nx; i++) {
811 if (fabs(x - xAxis[i]) < dist) {
812 dist = fabs(x - xAxis[i]);
813 iNode = i;
814 }
815 }
816 // Set the summing range.
817 iX0 = iNode;
818 iX1 = iNode;
819 // Establish the shape functions.
820 fX[0] = 0.;
821 fX[1] = 0.;
822 fX[2] = 0.;
823 fX[3] = 0.;
824 } else if (iOrder == 1 || nx == 2) {
825 // First order interpolation in x.
826 // Find the grid segment containing this point.
827 int iGrid = 0;
828 for (int i = 1; i < nx; i++) {
829 if ((xAxis[i - 1] - x) * (x - xAxis[i]) >= 0.) {
830 iGrid = i;
831 }
832 }
833 // Ensure there won't be divisions by zero.
834 if (xAxis[iGrid] == xAxis[iGrid - 1]) {
835 std::cerr << "Boxin3:\n";
836 std::cerr << " Incorrect grid; no interpolation.\n";
837 f = 0.;
838 return false;
839 }
840 // Compute local coordinates.
841 const double xLocal =
842 (x - xAxis[iGrid - 1]) / (xAxis[iGrid] - xAxis[iGrid - 1]);
843 // Set the summing range.
844 iX0 = iGrid - 1;
845 iX1 = iGrid;
846 // Set the shape functions.
847 fX[0] = 1. - xLocal;
848 fX[1] = xLocal;
849 fX[2] = 0.;
850 fX[3] = 0.;
851 } else if (iOrder == 2) {
852 // Second order interpolation in x.
853 // Find the grid segment containing this point.
854 int iGrid = 0;
855 for (int i = 1; i < nx; i++) {
856 if ((xAxis[i - 1] - x) * (x - xAxis[i]) >= 0.) {
857 iGrid = i;
858 }
859 }
860 // Compute the local coordinate for this grid segment.
861 const double xLocal =
862 (x - xAxis[iGrid - 1]) / (xAxis[iGrid] - xAxis[iGrid - 1]);
863 // Set the summing range and shape functions.
864 if (iGrid == 1) {
865 iX0 = iGrid - 1;
866 iX1 = iGrid + 1;
867 if (xAxis[iX0] == xAxis[iX0 + 1] || xAxis[iX0] == xAxis[iX0 + 2] ||
868 xAxis[iX0 + 1] == xAxis[iX0 + 2]) {
869 std::cerr << "Boxin3:\n";
870 std::cerr << " One or more grid points in x coincide.\n";
871 std::cerr << " No interpolation.\n";
872 f = 0.;
873 return false;
874 }
875 fX[0] = (x - xAxis[iX0 + 1]) * (x - xAxis[iX0 + 2]) /
876 ((xAxis[iX0] - xAxis[iX0 + 1]) * (xAxis[iX0] - xAxis[iX0 + 2]));
877 fX[1] =
878 (x - xAxis[iX0]) * (x - xAxis[iX0 + 2]) /
879 ((xAxis[iX0 + 1] - xAxis[iX0]) * (xAxis[iX0 + 1] - xAxis[iX0 + 2]));
880 fX[2] =
881 (x - xAxis[iX0]) * (x - xAxis[iX0 + 1]) /
882 ((xAxis[iX0 + 2] - xAxis[iX0]) * (xAxis[iX0 + 2] - xAxis[iX0 + 1]));
883 } else if (iGrid == nx - 1) {
884 iX0 = iGrid - 2;
885 iX1 = iGrid;
886 if (xAxis[iX0] == xAxis[iX0 + 1] || xAxis[iX0] == xAxis[iX0 + 2] ||
887 xAxis[iX0 + 1] == xAxis[iX0 + 2]) {
888 std::cerr << "Boxin3:\n";
889 std::cerr << " One or more grid points in x coincide.\n";
890 std::cerr << " No interpolation.\n";
891 f = 0.;
892 return false;
893 }
894 fX[0] = (x - xAxis[iX0 + 1]) * (x - xAxis[iX0 + 2]) /
895 ((xAxis[iX0] - xAxis[iX0 + 1]) * (xAxis[iX0] - xAxis[iX0 + 2]));
896 fX[1] =
897 (x - xAxis[iX0]) * (x - xAxis[iX0 + 2]) /
898 ((xAxis[iX0 + 1] - xAxis[iX0]) * (xAxis[iX0 + 1] - xAxis[iX0 + 2]));
899 fX[2] =
900 (x - xAxis[iX0]) * (x - xAxis[iX0 + 1]) /
901 ((xAxis[iX0 + 2] - xAxis[iX0]) * (xAxis[iX0 + 2] - xAxis[iX0 + 1]));
902 } else {
903 iX0 = iGrid - 2;
904 iX1 = iGrid + 1;
905 if (xAxis[iX0] == xAxis[iX0 + 1] || xAxis[iX0] == xAxis[iX0 + 2] ||
906 xAxis[iX0] == xAxis[iX0 + 3] || xAxis[iX0 + 1] == xAxis[iX0 + 2] ||
907 xAxis[iX0 + 1] == xAxis[iX0 + 3] ||
908 xAxis[iX0 + 2] == xAxis[iX0 + 3]) {
909 std::cerr << "Boxin3:\n";
910 std::cerr << " One or more grid points in x coincide.\n";
911 std::cerr << " No interpolation.\n";
912 f = 0.;
913 return false;
914 }
915 fX[0] = (x - xAxis[iX0 + 1]) * (x - xAxis[iX0 + 2]) /
916 ((xAxis[iX0] - xAxis[iX0 + 1]) * (xAxis[iX0] - xAxis[iX0 + 2]));
917 fX[1] =
918 (x - xAxis[iX0]) * (x - xAxis[iX0 + 2]) /
919 ((xAxis[iX0 + 1] - xAxis[iX0]) * (xAxis[iX0 + 1] - xAxis[iX0 + 2]));
920 fX[2] =
921 (x - xAxis[iX0]) * (x - xAxis[iX0 + 1]) /
922 ((xAxis[iX0 + 2] - xAxis[iX0]) * (xAxis[iX0 + 2] - xAxis[iX0 + 1]));
923 fX[0] *= (1. - xLocal);
924 fX[1] = fX[1] * (1. - xLocal) + xLocal * (x - xAxis[iX0 + 2]) *
925 (x - xAxis[iX0 + 3]) /
926 ((xAxis[iX0 + 1] - xAxis[iX0 + 2]) *
927 (xAxis[iX0 + 1] - xAxis[iX0 + 3]));
928 fX[2] = fX[2] * (1. - xLocal) + xLocal * (x - xAxis[iX0 + 1]) *
929 (x - xAxis[iX0 + 3]) /
930 ((xAxis[iX0 + 2] - xAxis[iX0 + 1]) *
931 (xAxis[iX0 + 2] - xAxis[iX0 + 3]));
932 fX[3] = xLocal * (x - xAxis[iX0 + 1]) * (x - xAxis[iX0 + 2]) /
933 ((xAxis[iX0 + 3] - xAxis[iX0 + 1]) *
934 (xAxis[iX0 + 3] - xAxis[iX0 + 2]));
935 }
936 }
937
938 if (iOrder == 0 || ny == 1) {
939 // Zeroth order interpolation in y.
940 // Find the nearest node.
941 double dist = fabs(y - yAxis[0]);
942 int iNode = 0;
943 for (int i = 1; i < ny; i++) {
944 if (fabs(y - yAxis[i]) < dist) {
945 dist = fabs(y - yAxis[i]);
946 iNode = i;
947 }
948 }
949 // Set the summing range.
950 iY0 = iNode;
951 iY1 = iNode;
952 // Establish the shape functions.
953 fY[0] = 1.;
954 fY[1] = 0.;
955 fY[2] = 0.;
956 } else if (iOrder == 1 || ny == 2) {
957 // First order interpolation in y.
958 // Find the grid segment containing this point.
959 int iGrid = 0;
960 for (int i = 1; i < ny; i++) {
961 if ((yAxis[i - 1] - y) * (y - yAxis[i]) >= 0.) {
962 iGrid = i;
963 }
964 }
965 // Ensure there won't be divisions by zero.
966 if (yAxis[iGrid] == yAxis[iGrid - 1]) {
967 std::cerr << "Boxin3:\n";
968 std::cerr << " Incorrect grid; no interpolation.\n";
969 f = 0.;
970 return false;
971 }
972 // Compute local coordinates.
973 const double yLocal =
974 (y - yAxis[iGrid - 1]) / (yAxis[iGrid] - yAxis[iGrid - 1]);
975 // Set the summing range.
976 iY0 = iGrid - 1;
977 iY1 = iGrid;
978 // Set the shape functions.
979 fY[0] = 1. - yLocal;
980 fY[1] = yLocal;
981 fY[2] = 0.;
982 } else if (iOrder == 2) {
983 // Second order interpolation in y.
984 // Find the grid segment containing this point.
985 int iGrid = 0;
986 for (int i = 1; i < ny; i++) {
987 if ((yAxis[i - 1] - y) * (y - yAxis[i]) >= 0.) {
988 iGrid = i;
989 }
990 }
991 // Compute the local coordinate for this grid segment.
992 const double yLocal =
993 (y - yAxis[iGrid - 1]) / (yAxis[iGrid] - yAxis[iGrid - 1]);
994 // Set the summing range and shape functions.
995 // These assignments are shared by all of the following conditions,
996 // so it's easier to take them out.
997 fY[0] = (y - yAxis[iY0 + 1]) * (y - yAxis[iY0 + 2]) /
998 ((yAxis[iY0] - yAxis[iY0 + 1]) * (yAxis[iY0] - yAxis[iY0 + 2]));
999 fY[1] = (y - yAxis[iY0]) * (y - yAxis[iY0 + 2]) /
1000 ((yAxis[iY0 + 1] - yAxis[iY0]) * (yAxis[iY0 + 1] - yAxis[iY0 + 2]));
1001 fY[2] = (y - yAxis[iY0]) * (y - yAxis[iY0 + 1]) /
1002 ((yAxis[iY0 + 2] - yAxis[iY0]) * (yAxis[iY0 + 2] - yAxis[iY0 + 1]));
1003
1004 if (iGrid == 1) {
1005 iY0 = iGrid - 1;
1006 iY1 = iGrid + 1;
1007 if (yAxis[iY0] == yAxis[iY0 + 1] || yAxis[iY0] == yAxis[iY0 + 2] ||
1008 yAxis[iY0 + 1] == yAxis[iY0 + 2]) {
1009 std::cerr << "Boxin3:\n";
1010 std::cerr << " One or more grid points in y coincide.\n";
1011 std::cerr << " No interpolation.\n";
1012 f = 0.;
1013 return false;
1014 }
1015 fY[0] = (y - yAxis[iY0 + 1]) * (y - yAxis[iY0 + 2]) /
1016 ((yAxis[iY0] - yAxis[iY0 + 1]) * (yAxis[iY0] - yAxis[iY0 + 2]));
1017 fY[1] =
1018 (y - yAxis[iY0]) * (y - yAxis[iY0 + 2]) /
1019 ((yAxis[iY0 + 1] - yAxis[iY0]) * (yAxis[iY0 + 1] - yAxis[iY0 + 2]));
1020 fY[2] =
1021 (y - yAxis[iY0]) * (y - yAxis[iY0 + 1]) /
1022 ((yAxis[iY0 + 2] - yAxis[iY0]) * (yAxis[iY0 + 2] - yAxis[iY0 + 1]));
1023 } else if (iGrid == ny - 1) {
1024 iY0 = iGrid - 2;
1025 iY1 = iGrid;
1026 if (yAxis[iY0] == yAxis[iY0 + 1] || yAxis[iY0] == yAxis[iY0 + 2] ||
1027 yAxis[iY0 + 1] == yAxis[iY0 + 2]) {
1028 std::cerr << "Boxin3:\n";
1029 std::cerr << " One or more grid points in y coincide.\n";
1030 std::cerr << " No interpolation.\n";
1031 f = 0.;
1032 return false;
1033 }
1034 fY[0] = (y - yAxis[iY0 + 1]) * (y - yAxis[iY0 + 2]) /
1035 ((yAxis[iY0] - yAxis[iY0 + 1]) * (yAxis[iY0] - yAxis[iY0 + 2]));
1036 fY[1] =
1037 (y - yAxis[iY0]) * (y - yAxis[iY0 + 2]) /
1038 ((yAxis[iY0 + 1] - yAxis[iY0]) * (yAxis[iY0 + 1] - yAxis[iY0 + 2]));
1039 fY[2] =
1040 (y - yAxis[iY0]) * (y - yAxis[iY0 + 1]) /
1041 ((yAxis[iY0 + 2] - yAxis[iY0]) * (yAxis[iY0 + 2] - yAxis[iY0 + 1]));
1042 } else {
1043 iY0 = iGrid - 2;
1044 iY1 = iGrid + 1;
1045 if (yAxis[iY0] == yAxis[iY0 + 1] || yAxis[iY0] == yAxis[iY0 + 2] ||
1046 yAxis[iY0] == yAxis[iY0 + 3] || yAxis[iY0 + 1] == yAxis[iY0 + 2] ||
1047 yAxis[iY0 + 1] == yAxis[iY0 + 3] ||
1048 yAxis[iY0 + 2] == yAxis[iY0 + 3]) {
1049 std::cerr << "Boxin3:\n";
1050 std::cerr << " One or more grid points in y coincide.\n";
1051 std::cerr << " No interpolation.\n";
1052 f = 0.;
1053 return false;
1054 }
1055 fY[0] = (y - yAxis[iY0 + 1]) * (y - yAxis[iY0 + 2]) /
1056 ((yAxis[iY0] - yAxis[iY0 + 1]) * (yAxis[iY0] - yAxis[iY0 + 2]));
1057 fY[1] =
1058 (y - yAxis[iY0]) * (y - yAxis[iY0 + 2]) /
1059 ((yAxis[iY0 + 1] - yAxis[iY0]) * (yAxis[iY0 + 1] - yAxis[iY0 + 2]));
1060 fY[2] =
1061 (y - yAxis[iY0]) * (y - yAxis[iY0 + 1]) /
1062 ((yAxis[iY0 + 2] - yAxis[iY0]) * (yAxis[iY0 + 2] - yAxis[iY0 + 1]));
1063
1064 fY[0] *= (1. - yLocal);
1065 fY[1] = fY[1] * (1. - yLocal) + yLocal * (y - yAxis[iY0 + 2]) *
1066 (y - yAxis[iY0 + 3]) /
1067 ((yAxis[iY0 + 1] - yAxis[iY0 + 2]) *
1068 (yAxis[iY0 + 1] - yAxis[iY0 + 3]));
1069 fY[2] = fY[2] * (1. - yLocal) + yLocal * (y - yAxis[iY0 + 1]) *
1070 (y - yAxis[iY0 + 3]) /
1071 ((yAxis[iY0 + 2] - yAxis[iY0 + 1]) *
1072 (yAxis[iY0 + 2] - yAxis[iY0 + 3]));
1073 fY[3] = yLocal * (y - yAxis[iY0 + 1]) * (y - yAxis[iY0 + 2]) /
1074 ((yAxis[iY0 + 3] - yAxis[iY0 + 1]) *
1075 (yAxis[iY0 + 3] - yAxis[iY0 + 2]));
1076 }
1077 }
1078
1079 if (iOrder == 0 || nz == 1) {
1080 // Zeroth order interpolation in z.
1081 // Find the nearest node.
1082 double dist = fabs(z - zAxis[0]);
1083 int iNode = 0;
1084 for (int i = 1; i < nz; i++) {
1085 if (fabs(z - zAxis[i]) < dist) {
1086 dist = fabs(z - zAxis[i]);
1087 iNode = i;
1088 }
1089 }
1090 // Set the summing range.
1091 iZ0 = iNode;
1092 iZ1 = iNode;
1093 // Establish the shape functions.
1094 fZ[0] = 1.;
1095 fZ[1] = 0.;
1096 fZ[2] = 0.;
1097 } else if (iOrder == 1 || nz == 2) {
1098 // First order interpolation in z.
1099 // Find the grid segment containing this point.
1100 int iGrid = 0;
1101 for (int i = 1; i < nz; i++) {
1102 if ((zAxis[i - 1] - z) * (z - zAxis[i]) >= 0.) {
1103 iGrid = i;
1104 }
1105 }
1106 // Ensure there won't be divisions by zero.
1107 if (zAxis[iGrid] == zAxis[iGrid - 1]) {
1108 std::cerr << "Boxin3:\n";
1109 std::cerr << " Incorrect grid; no interpolation.\n";
1110 f = 0.;
1111 return false;
1112 }
1113 // Compute local coordinates.
1114 const double zLocal =
1115 (z - zAxis[iGrid - 1]) / (zAxis[iGrid] - zAxis[iGrid - 1]);
1116 // Set the summing range.
1117 iZ0 = iGrid - 1;
1118 iZ1 = iGrid;
1119 // Set the shape functions.
1120 fZ[0] = 1. - zLocal;
1121 fZ[1] = zLocal;
1122 fZ[2] = 0.;
1123 } else if (iOrder == 2) {
1124 // Second order interpolation in z.
1125 // Find the grid segment containing this point.
1126 int iGrid = 0;
1127 for (int i = 1; i < nz; i++) {
1128 if ((zAxis[i - 1] - z) * (z - zAxis[i]) >= 0.) {
1129 iGrid = i;
1130 }
1131 }
1132 // Compute the local coordinate for this grid segment.
1133 const double zLocal =
1134 (z - zAxis[iGrid - 1]) / (zAxis[iGrid] - zAxis[iGrid - 1]);
1135 // Set the summing range and shape functions.
1136 // These assignments are shared by all of the following conditions,
1137 // so it's easier to take them out.
1138 fZ[0] = (z - zAxis[iZ0 + 1]) * (z - zAxis[iZ0 + 2]) /
1139 ((zAxis[iZ0] - zAxis[iZ0 + 1]) * (zAxis[iZ0] - zAxis[iZ0 + 2]));
1140 fZ[1] = (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 2]) /
1141 ((zAxis[iZ0 + 1] - zAxis[iZ0]) * (zAxis[iZ0 + 1] - zAxis[iZ0 + 2]));
1142 fZ[2] = (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 1]) /
1143 ((zAxis[iZ0 + 2] - zAxis[iZ0]) * (zAxis[iZ0 + 2] - zAxis[iZ0 + 1]));
1144
1145 if (iGrid == 1) {
1146 iZ0 = iGrid - 1;
1147 iZ1 = iGrid + 1;
1148 if (zAxis[iZ0] == zAxis[iZ0 + 1] || zAxis[iZ0] == zAxis[iZ0 + 2] ||
1149 zAxis[iZ0 + 1] == zAxis[iZ0 + 2]) {
1150 std::cerr << "Boxin3:\n";
1151 std::cerr << " One or more grid points in z coincide.\n";
1152 std::cerr << " No interpolation.\n";
1153 f = 0.;
1154 return false;
1155 }
1156 fZ[0] = (z - zAxis[iZ0 + 1]) * (z - zAxis[iZ0 + 2]) /
1157 ((zAxis[iZ0] - zAxis[iZ0 + 1]) * (zAxis[iZ0] - zAxis[iZ0 + 2]));
1158 fZ[1] =
1159 (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 2]) /
1160 ((zAxis[iZ0 + 1] - zAxis[iZ0]) * (zAxis[iZ0 + 1] - zAxis[iZ0 + 2]));
1161 fZ[2] =
1162 (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 1]) /
1163 ((zAxis[iZ0 + 2] - zAxis[iZ0]) * (zAxis[iZ0 + 2] - zAxis[iZ0 + 1]));
1164 } else if (iGrid == nz - 1) {
1165 iZ0 = iGrid - 2;
1166 iZ1 = iGrid;
1167 if (zAxis[iZ0] == zAxis[iZ0 + 1] || zAxis[iZ0] == zAxis[iZ0 + 2] ||
1168 zAxis[iZ0 + 1] == zAxis[iZ0 + 2]) {
1169 std::cerr << "Boxin3:\n";
1170 std::cerr << " One or more grid points in z coincide.\n";
1171 std::cerr << " No interpolation.\n";
1172 f = 0.;
1173 return false;
1174 }
1175 fZ[0] = (z - zAxis[iZ0 + 1]) * (z - zAxis[iZ0 + 2]) /
1176 ((zAxis[iZ0] - zAxis[iZ0 + 1]) * (zAxis[iZ0] - zAxis[iZ0 + 2]));
1177 fZ[1] =
1178 (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 2]) /
1179 ((zAxis[iZ0 + 1] - zAxis[iZ0]) * (zAxis[iZ0 + 1] - zAxis[iZ0 + 2]));
1180 fZ[2] =
1181 (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 1]) /
1182 ((zAxis[iZ0 + 2] - zAxis[iZ0]) * (zAxis[iZ0 + 2] - zAxis[iZ0 + 1]));
1183 } else {
1184 iZ0 = iGrid - 2;
1185 iZ1 = iGrid + 1;
1186
1187 if (zAxis[iZ0] == zAxis[iZ0 + 1] || zAxis[iZ0] == zAxis[iZ0 + 2] ||
1188 zAxis[iZ0] == zAxis[iZ0 + 3] || zAxis[iZ0 + 1] == zAxis[iZ0 + 2] ||
1189 zAxis[iZ0 + 1] == zAxis[iZ0 + 3] ||
1190 zAxis[iZ0 + 2] == zAxis[iZ0 + 3]) {
1191 std::cerr << "Boxin3:\n";
1192 std::cerr << " One or more grid points in z coincide.\n";
1193 std::cerr << " No interpolation.\n";
1194 f = 0.;
1195 return false;
1196 }
1197
1198 fZ[0] = (z - zAxis[iZ0 + 1]) * (z - zAxis[iZ0 + 2]) /
1199 ((zAxis[iZ0] - zAxis[iZ0 + 1]) * (zAxis[iZ0] - zAxis[iZ0 + 2]));
1200 fZ[1] =
1201 (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 2]) /
1202 ((zAxis[iZ0 + 1] - zAxis[iZ0]) * (zAxis[iZ0 + 1] - zAxis[iZ0 + 2]));
1203 fZ[2] =
1204 (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 1]) /
1205 ((zAxis[iZ0 + 2] - zAxis[iZ0]) * (zAxis[iZ0 + 2] - zAxis[iZ0 + 1]));
1206
1207 fZ[0] *= (1. - zLocal);
1208 fZ[1] = fZ[1] * (1. - zLocal) + zLocal * (z - zAxis[iZ0 + 2]) *
1209 (z - zAxis[iZ0 + 3]) /
1210 ((zAxis[iZ0 + 1] - zAxis[iZ0 + 2]) *
1211 (zAxis[iZ0 + 1] - zAxis[iZ0 + 3]));
1212 fZ[2] = fZ[2] * (1. - zLocal) + zLocal * (z - zAxis[iZ0 + 1]) *
1213 (z - zAxis[iZ0 + 3]) /
1214 ((zAxis[iZ0 + 2] - zAxis[iZ0 + 1]) *
1215 (zAxis[iZ0 + 2] - zAxis[iZ0 + 3]));
1216 fZ[3] = zLocal * (z - zAxis[iZ0 + 1]) * (z - zAxis[iZ0 + 2]) /
1217 ((zAxis[iZ0 + 3] - zAxis[iZ0 + 1]) *
1218 (zAxis[iZ0 + 3] - zAxis[iZ0 + 2]));
1219 }
1220 }
1221
1222 f = 0.;
1223 for (int i = iX0; i <= iX1; ++i) {
1224 for (int j = iY0; j <= iY1; ++j) {
1225 for (int k = iZ0; k <= iZ1; ++k) {
1226 // std::cout << "i = " << i << ", j = " << j << ", k = " << k << "\n";
1227 // std::cout << "value: " << value[i][j][k] << "\n";
1228 // std::cout << "fX = " << fX[i - iX0] << ", fY = " << fY[j - iY0] << ",
1229 // fZ = " << fZ[k - iZ0] << "\n";
1230 f += value[i][j][k] * fX[i - iX0] * fY[j - iY0] * fZ[k - iZ0];
1231 }
1232 }
1233 }
1234 // std::cout << f << std::endl;
1235 return true;
1236}
fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:616

Referenced by Garfield::Medium::CloneTable(), Garfield::Medium::CloneTensor(), Garfield::Medium::ElectronAttachment(), Garfield::Medium::ElectronDiffusion(), Garfield::Medium::ElectronTownsend(), Garfield::Medium::ElectronVelocity(), Garfield::Medium::HoleAttachment(), Garfield::Medium::HoleDiffusion(), Garfield::Medium::HoleTownsend(), Garfield::Medium::HoleVelocity(), Garfield::Medium::IonDiffusion(), Garfield::Medium::IonDissociation(), and Garfield::Medium::IonVelocity().

◆ Cfact()

void Garfield::Numerics::Cfact ( const int  n,
std::vector< std::vector< std::complex< double > > > &  a,
std::vector< int > &  ir,
int &  ifail,
std::complex< double > &  det,
int &  jfail 
)

Definition at line 333 of file Numerics.cc.

335 {
336
337 const double g1 = 1.e-19;
338 const double g2 = 1.e-19;
339
340 std::complex<double> tf;
341 double p, q, t;
342 std::complex<double> s11, s12;
343 int k;
344
345 ifail = jfail = 0;
346
347 int nxch = 0;
348 det = std::complex<double>(1., 0.);
349
350 for (int j = 1; j <= n; ++j) {
351 k = j;
352 p = std::max(fabs(real(a[j - 1][j - 1])), fabs(imag(a[j - 1][j - 1])));
353 if (j == n) {
354 if (p <= 0.) {
355 det = std::complex<double>(0., 0.);
356 ifail = -1;
357 jfail = 0;
358 return;
359 }
360 det *= a[j - 1][j - 1];
361 a[j - 1][j - 1] = std::complex<double>(1., 0.) / a[j - 1][j - 1];
362 t = std::max(fabs(real(det)), fabs(imag(det)));
363 if (t < g1) {
364 det = std::complex<double>(0., 0.);
365 if (jfail == 0) jfail = -1;
366 } else if (t > g2) {
367 det = std::complex<double>(1., 0.);
368 if (jfail == 0) jfail = +1;
369 }
370 continue;
371 }
372 for (int i = j + 1; i <= n; ++i) {
373 q = std::max(fabs(real(a[i - 1][j - 1])), fabs(imag(a[i - 1][j - 1])));
374 if (q <= p) continue;
375 k = i;
376 p = q;
377 }
378 if (k != j) {
379 for (int l = 1; l <= n; ++l) {
380 tf = a[j - 1][l - 1];
381 a[j - 1][l - 1] = a[k - 1][l - 1];
382 a[k - 1][l - 1] = tf;
383 }
384 ++nxch;
385 ir[nxch - 1] = j * 4096 + k;
386 } else if (p <= 0.) {
387 det = std::complex<double>(0., 0.);
388 ifail = -1;
389 jfail = 0;
390 return;
391 }
392 det *= a[j - 1][j - 1];
393 a[j - 1][j - 1] = 1. / a[j - 1][j - 1];
394 t = std::max(fabs(real(det)), fabs(imag(det)));
395 if (t < g1) {
396 det = std::complex<double>(0., 0.);
397 if (jfail == 0) jfail = -1;
398 } else if (t > g2) {
399 det = std::complex<double>(1., 0.);
400 if (jfail == 0) jfail = +1;
401 }
402 for (k = j + 1; k <= n; ++k) {
403 s11 = -a[j - 1][k - 1];
404 s12 = -a[k - 1][j];
405 if (j == 1) {
406 a[j - 1][k - 1] = -s11 * a[j - 1][j - 1];
407 a[k - 1][j] = -(s12 + a[j - 1][j] * a[k - 1][j - 1]);
408 continue;
409 }
410 for (int i = 1; i <= j - 1; ++i) {
411 s11 += a[i - 1][k - 1] * a[j - 1][i - 1];
412 s12 += a[i - 1][j] * a[k - 1][i - 1];
413 }
414 a[j - 1][k - 1] = -s11 * a[j - 1][j - 1];
415 a[k - 1][j] = -a[j - 1][j] * a[k - 1][j - 1] - s12;
416 }
417 }
418
419 if (nxch % 2 != 0) det = -det;
420 if (jfail != 0) det = std::complex<double>(0., 0.);
421 ir[n - 1] = nxch;
422}

Referenced by Cinv().

◆ Cfinv()

void Garfield::Numerics::Cfinv ( const int  n,
std::vector< std::vector< std::complex< double > > > &  a,
std::vector< int > &  ir 
)

Definition at line 424 of file Numerics.cc.

425 {
426
427 std::complex<double> ti;
428 std::complex<double> s31, s32, s33, s34;
429
430 if (n <= 1) return;
431 a[1][0] = -a[1][1] * a[0][0] * a[1][0];
432 a[0][1] = -a[0][1];
433 if (n > 2) {
434 for (int i = 3; i <= n; ++i) {
435 for (int j = 1; j <= i - 2; ++j) {
436 s31 = std::complex<double>(0., 0.);
437 s32 = a[j - 1][i - 1];
438 for (int k = j; k <= i - 2; ++k) {
439 s31 += a[k - 1][j - 1] * a[i - 1][k - 1];
440 s32 += a[j - 1][k] * a[k][i - 1];
441 }
442 a[i - 1][j - 1] =
443 -a[i - 1][i - 1] * (s31 + a[i - 2][j - 1] * a[i - 1][i - 2]);
444 a[j - 1][i - 1] = -s32;
445 }
446 a[i - 1][i - 2] = -a[i - 1][i - 1] * a[i - 2][i - 2] * a[i - 1][i - 2];
447 a[i - 2][i - 1] = -a[i - 2][i - 1];
448 }
449 }
450
451 for (int i = 1; i <= n - 1; ++i) {
452 for (int j = 1; j <= i; ++j) {
453 s33 = a[i - 1][j - 1];
454 for (int k = 1; k <= n - i; ++k) {
455 s33 += a[i + k - 1][j - 1] * a[i - 1][i + k - 1];
456 }
457 a[i - 1][j - 1] = s33;
458 }
459 for (int j = 1; j <= n - i; ++j) {
460 s34 = std::complex<double>(0., 0.);
461 for (int k = j; k <= n - i; ++k) {
462 s34 += a[i + k - 1][i + j - 1] * a[i - 1][i + k - 1];
463 }
464 a[i - 1][i + j - 1] = s34;
465 }
466 }
467
468 int nxch = ir[n - 1];
469 if (nxch == 0) return;
470
471 for (int m = 1; m <= nxch; ++m) {
472 int k = nxch - m + 1;
473 int ij = ir[k - 1];
474 int i = ij / 4096;
475 int j = ij % 4096;
476 for (k = 1; k <= n; ++k) {
477 ti = a[k - 1][i - 1];
478 a[k - 1][i - 1] = a[k - 1][j - 1];
479 a[k - 1][j - 1] = ti;
480 }
481 }
482}

Referenced by Cinv().

◆ Cinv()

void Garfield::Numerics::Cinv ( const int  n,
std::vector< std::vector< std::complex< double > > > &  a,
int &  ifail 
)

Definition at line 494 of file Numerics.cc.

495 {
496
497 double t1, t2, t3;
498 std::complex<double> det, temp, s;
499 std::complex<double> c11, c12, c13, c21, c22, c23, c31, c32, c33;
500
501 std::vector<int> ir;
502 ir.clear();
503 ir.resize(n);
504 for (int i = 0; i < n; ++i) ir[i] = 0;
505
506 // TEST FOR PARAMETER ERRORS.
507 if (n < 1) {
508 ifail = 1;
509 return;
510 }
511
512 ifail = 0;
513 int jfail = 0;
514
515 if (n > 3) {
516 // n > 3 CASES.
517 // FACTORIZE MATRIX AND INVERT.
518 Cfact(n, a, ir, ifail, det, jfail);
519 if (ifail != 0) return;
520 Cfinv(n, a, ir);
521 } else if (n == 3) {
522 // n = 3 CASE.
523 // COMPUTE COFACTORS.
524 c11 = a[1][1] * a[2][2] - a[1][2] * a[2][1];
525 c12 = a[1][2] * a[2][0] - a[1][0] * a[2][2];
526 c13 = a[1][0] * a[2][1] - a[1][1] * a[2][0];
527 c21 = a[2][1] * a[0][2] - a[2][2] * a[0][1];
528 c22 = a[2][2] * a[0][0] - a[2][0] * a[0][2];
529 c23 = a[2][0] * a[0][1] - a[2][1] * a[0][0];
530 c31 = a[0][1] * a[1][2] - a[0][2] * a[1][1];
531 c32 = a[0][2] * a[1][0] - a[0][0] * a[1][2];
532 c33 = a[0][0] * a[1][1] - a[0][1] * a[1][0];
533 t1 = fabs(real(a[0][0])) + fabs(imag(a[0][0]));
534 t2 = fabs(real(a[1][0])) + fabs(imag(a[1][0]));
535 t3 = fabs(real(a[2][0])) + fabs(imag(a[2][0]));
536
537 // SET temp = PIVOT AND det = PIVOT * det.
538 if (t1 >= t2) {
539 if (t3 >= t1) {
540 // PIVOT IS A31
541 temp = a[2][0];
542 det = c23 * c12 - c22 * c13;
543 } else {
544 // PIVOT IS A11
545 temp = a[0][0];
546 det = c22 * c33 - c23 * c32;
547 }
548 } else {
549 if (t3 >= t2) {
550 // PIVOT IS A31
551 temp = a[2][0];
552 det = c23 * c12 - c22 * c13;
553 } else {
554 // PIVOT IS A21
555 temp = a[1][0];
556 det = c13 * c32 - c12 * c33;
557 }
558 }
559 // SET ELEMENTS OF INVERSE IN A.
560 if (real(det) == 0. && imag(det) == 0.) {
561 ifail = -1;
562 return;
563 }
564 s = temp / det;
565 a[0][0] = s * c11;
566 a[0][1] = s * c21;
567 a[0][2] = s * c31;
568 a[1][0] = s * c12;
569 a[1][1] = s * c22;
570 a[1][2] = s * c32;
571 a[2][0] = s * c13;
572 a[2][1] = s * c23;
573 a[2][2] = s * c33;
574 } else if (n == 2) {
575 // n=2 CASE BY CRAMERS RULE.
576 det = a[0][0] * a[1][1] - a[0][1] * a[1][0];
577 if (real(det) == 0. && imag(det) == 0.) {
578 ifail = -1;
579 return;
580 }
581 s = std::complex<double>(1., 0.) / det;
582 c11 = s * a[1][1];
583 a[0][1] = -s * a[0][1];
584 a[1][0] = -s * a[1][0];
585 a[1][1] = s * a[0][0];
586 a[0][0] = c11;
587 } else {
588 // n = 1 CASE.
589 if (real(a[0][0]) == 0. && imag(a[0][0]) == 0.) {
590 ifail = -1;
591 return;
592 }
593 a[0][0] = std::complex<double>(1., 0.) / a[0][0];
594 }
595}

◆ Deqinv()

void Garfield::Numerics::Deqinv ( const int  n,
std::vector< std::vector< double > > &  a,
int &  ifail,
std::vector< double > &  b 
)

Definition at line 216 of file Numerics.cc.

217 {
218
219 double t1, t2, t3;
220 double det, temp, s;
221 double b1, b2, c11, c12, c13, c21, c22, c23, c31, c32, c33;
222
223 std::vector<int> ir;
224 ir.clear();
225 ir.resize(n);
226 for (int i = 0; i < n; ++i) ir[i] = 0;
227
228 // TEST FOR PARAMETER ERRORS.
229 if (n < 1) {
230 ifail = +1;
231 return;
232 }
233
234 ifail = 0;
235 int jfail = 0;
236
237 if (n > 3) {
238 // n > 3 CASES.
239 // FACTORIZE MATRIX, INVERT AND SOLVE SYSTEM.
240 Dfact(n, a, ir, ifail, det, jfail);
241 if (ifail != 0) return;
242 Dfeqn(n, a, ir, b);
243 Dfinv(n, a, ir);
244 } else if (n == 3) {
245 // n = 3 CASE.
246 // COMPUTE COFACTORS.
247 c11 = a[1][1] * a[2][2] - a[1][2] * a[2][1];
248 c12 = a[1][2] * a[2][0] - a[1][0] * a[2][2];
249 c13 = a[1][0] * a[2][1] - a[1][1] * a[2][0];
250 c21 = a[2][1] * a[0][2] - a[2][2] * a[0][1];
251 c22 = a[2][2] * a[0][0] - a[2][0] * a[0][2];
252 c23 = a[2][0] * a[0][1] - a[2][1] * a[0][0];
253 c31 = a[0][1] * a[1][2] - a[0][2] * a[1][1];
254 c32 = a[0][2] * a[1][0] - a[0][0] * a[1][2];
255 c33 = a[0][0] * a[1][1] - a[0][1] * a[1][0];
256 t1 = fabs(a[0][0]);
257 t2 = fabs(a[1][0]);
258 t3 = fabs(a[2][0]);
259
260 // SET temp = PIVOT AND det = PIVOT * det.
261 if (t1 >= t2) {
262 if (t3 >= t1) {
263 // PIVOT IS A31
264 temp = a[2][0];
265 det = c23 * c12 - c22 * c13;
266 } else {
267 // PIVOT IS A11
268 temp = a[0][0];
269 det = c22 * c33 - c23 * c32;
270 }
271 } else {
272 if (t3 >= t2) {
273 // PIVOT IS A31
274 temp = a[2][0];
275 det = c23 * c12 - c22 * c13;
276 } else {
277 // PIVOT IS A21
278 temp = a[1][0];
279 det = c13 * c32 - c12 * c33;
280 }
281 }
282
283 // SET ELEMENTS OF INVERSE IN A.
284 if (det == 0.) {
285 ifail = -1;
286 return;
287 }
288 s = temp / det;
289 a[0][0] = s * c11;
290 a[0][1] = s * c21;
291 a[0][2] = s * c31;
292 a[1][0] = s * c12;
293 a[1][1] = s * c22;
294 a[1][2] = s * c32;
295 a[2][0] = s * c13;
296 a[2][1] = s * c23;
297 a[2][2] = s * c33;
298
299 // REPLACE B BY AINV*B.
300 b1 = b[0];
301 b2 = b[1];
302 b[0] = a[0][0] * b1 + a[0][1] * b2 + a[0][2] * b[2];
303 b[1] = a[1][0] * b1 + a[1][1] * b2 + a[1][2] * b[2];
304 b[2] = a[2][0] * b1 + a[2][1] * b2 + a[2][2] * b[2];
305 } else if (n == 2) {
306 // n = 2 CASE BY CRAMERS RULE.
307 det = a[0][0] * a[1][1] - a[0][1] * a[1][0];
308 if (det == 0.) {
309 ifail = -1;
310 return;
311 }
312 s = 1. / det;
313 c11 = s * a[1][1];
314 a[0][1] = -s * a[0][1];
315 a[1][0] = -s * a[1][0];
316 a[1][1] = s * a[0][0];
317 a[0][0] = c11;
318
319 b1 = b[0];
320 b[0] = c11 * b1 + a[0][1] * b[1];
321 b[1] = a[1][0] * b1 + a[1][1] * b[1];
322 } else {
323 // n = 1 CASE.
324 if (a[0][0] == 0.) {
325 ifail = -1;
326 return;
327 }
328 a[0][0] = 1. / a[0][0];
329 b[0] = a[0][0] * b[0];
330 }
331}

◆ Dfact()

void Garfield::Numerics::Dfact ( const int  n,
std::vector< std::vector< double > > &  a,
std::vector< int > &  ir,
int &  ifail,
double &  det,
int &  jfail 
)

Definition at line 10 of file Numerics.cc.

11 {
12
13 const double g1 = 1.e-19;
14 const double g2 = 1.e-19;
15
16 double tf, p, q, t, s11, s12;
17 int k;
18
19 ifail = jfail = 0;
20
21 int nxch = 0;
22 det = 1.;
23
24 for (int j = 1; j <= n; ++j) {
25 k = j;
26 p = fabs(a[j - 1][j - 1]);
27 if (j == n) {
28 if (p <= 0.) {
29 det = 0.;
30 ifail = -1;
31 jfail = 0;
32 return;
33 }
34 det *= a[j - 1][j - 1];
35 a[j - 1][j - 1] = 1. / a[j - 1][j - 1];
36 t = fabs(det);
37 if (t < g1) {
38 det = 0.;
39 if (jfail == 0) jfail = -1;
40 } else if (t > g2) {
41 det = 1.;
42 if (jfail == 0) jfail = +1;
43 }
44 continue;
45 }
46 for (int i = j + 1; i <= n; ++i) {
47 q = fabs(a[i - 1][j - 1]);
48 if (q <= p) continue;
49 k = i;
50 p = q;
51 }
52 if (k != j) {
53 for (int l = 1; l <= n; ++l) {
54 tf = a[j - 1][l - 1];
55 a[j - 1][l - 1] = a[k - 1][l - 1];
56 a[k - 1][l - 1] = tf;
57 }
58 ++nxch;
59 ir[nxch - 1] = j * 4096 + k;
60 } else if (p <= 0.) {
61 det = 0.;
62 ifail = -1;
63 jfail = 0;
64 return;
65 }
66 det *= a[j - 1][j - 1];
67 a[j - 1][j - 1] = 1. / a[j - 1][j - 1];
68 t = fabs(det);
69 if (t < g1) {
70 det = 0.;
71 if (jfail == 0) jfail = -1;
72 } else if (t > g2) {
73 det = 1.;
74 if (jfail == 0) jfail = +1;
75 }
76 for (k = j + 1; k <= n; ++k) {
77 s11 = -a[j - 1][k - 1];
78 s12 = -a[k - 1][j];
79 if (j == 1) {
80 a[j - 1][k - 1] = -s11 * a[j - 1][j - 1];
81 a[k - 1][j] = -(s12 + a[j - 1][j] * a[k - 1][j - 1]);
82 continue;
83 }
84 for (int i = 1; i <= j - 1; ++i) {
85 s11 += a[i - 1][k - 1] * a[j - 1][i - 1];
86 s12 += a[i - 1][j] * a[k - 1][i - 1];
87 }
88 a[j - 1][k - 1] = -s11 * a[j - 1][j - 1];
89 a[k - 1][j] = -a[j - 1][j] * a[k - 1][j - 1] - s12;
90 }
91 }
92
93 if (nxch % 2 != 0) det = -det;
94 if (jfail != 0) det = 0.;
95 ir[n - 1] = nxch;
96}

Referenced by Deqinv().

◆ Dfeqn()

void Garfield::Numerics::Dfeqn ( const int  n,
std::vector< std::vector< double > > &  a,
std::vector< int > &  ir,
std::vector< double > &  b 
)

Definition at line 98 of file Numerics.cc.

99 {
100
101 double te;
102 double s21, s22;
103
104 if (n <= 0) return;
105
106 int nxch = ir[n - 1];
107 if (nxch != 0) {
108 for (int m = 1; m <= nxch; ++m) {
109 int ij = ir[m - 1];
110 int i = ij / 4096;
111 int j = ij % 4096;
112 te = b[i - 1];
113 b[i - 1] = b[j - 1];
114 b[j - 1] = te;
115 }
116 }
117
118 b[0] *= a[0][0];
119 if (n == 1) return;
120
121 for (int i = 2; i <= n; ++i) {
122 s21 = -b[i - 1];
123 for (int j = 1; j <= i - 1; ++j) {
124 s21 += a[i - 1][j - 1] * b[j - 1];
125 }
126 b[i - 1] = -a[i - 1][i - 1] * s21;
127 }
128
129 for (int i = 1; i <= n - 1; ++i) {
130 s22 = -b[n - i - 1];
131 for (int j = 1; j <= i; ++j) {
132 s22 += a[n - i - 1][n - j] * b[n - j];
133 }
134 b[n - i - 1] = -s22;
135 }
136}

Referenced by Deqinv().

◆ Dfinv()

void Garfield::Numerics::Dfinv ( const int  n,
std::vector< std::vector< double > > &  a,
std::vector< int > &  ir 
)

Definition at line 138 of file Numerics.cc.

139 {
140
141 double ti;
142 double s31, s32, s33, s34;
143
144 if (n <= 1) return;
145 a[1][0] = -a[1][1] * a[0][0] * a[1][0];
146 a[0][1] = -a[0][1];
147 if (n > 2) {
148 for (int i = 3; i <= n; ++i) {
149 for (int j = 1; j <= i - 2; ++j) {
150 s31 = 0.;
151 s32 = a[j - 1][i - 1];
152 for (int k = j; k <= i - 2; ++k) {
153 s31 += a[k - 1][j - 1] * a[i - 1][k - 1];
154 s32 += a[j - 1][k] * a[k][i - 1];
155 }
156 a[i - 1][j - 1] =
157 -a[i - 1][i - 1] * (s31 + a[i - 2][j - 1] * a[i - 1][i - 2]);
158 a[j - 1][i - 1] = -s32;
159 }
160 a[i - 1][i - 2] = -a[i - 1][i - 1] * a[i - 2][i - 2] * a[i - 1][i - 2];
161 a[i - 2][i - 1] = -a[i - 2][i - 1];
162 }
163 }
164
165 for (int i = 1; i <= n - 1; ++i) {
166 for (int j = 1; j <= i; ++j) {
167 s33 = a[i - 1][j - 1];
168 for (int k = 1; k <= n - i; ++k) {
169 s33 += a[i + k - 1][j - 1] * a[i - 1][i + k - 1];
170 }
171 a[i - 1][j - 1] = s33;
172 }
173 for (int j = 1; j <= n - i; ++j) {
174 s34 = 0.;
175 for (int k = j; k <= n - i; ++k) {
176 s34 += a[i + k - 1][i + j - 1] * a[i - 1][i + k - 1];
177 }
178 a[i - 1][i + j - 1] = s34;
179 }
180 }
181
182 int nxch = ir[n - 1];
183 if (nxch == 0) return;
184
185 for (int m = 1; m <= nxch; ++m) {
186 int k = nxch - m + 1;
187 int ij = ir[k - 1];
188 int i = ij / 4096;
189 int j = ij % 4096;
190 for (k = 1; k <= n; ++k) {
191 ti = a[k - 1][i - 1];
192 a[k - 1][i - 1] = a[k - 1][j - 1];
193 a[k - 1][j - 1] = ti;
194 }
195 }
196}

Referenced by Deqinv().

◆ Divdif()

double Garfield::Numerics::Divdif ( const std::vector< double > &  f,
const std::vector< double > &  a,
int  nn,
double  x,
int  mm 
)

Definition at line 650 of file Numerics.cc.

651 {
652
653 // C++ version of DIVDIF (CERN program library E105) which performs
654 // tabular interpolation using symmetrically placed argument points.
655
656 double t[20], d[20];
657
658 const int mmax = 10;
659
660 // Check the arguments.
661 if (nn < 2) {
662 std::cerr << "Divdif:\n";
663 std::cerr << " Array length < 2.\n";
664 return 0.;
665 }
666 if (mm < 1) {
667 std::cerr << "Divdif:\n";
668 std::cerr << " Interpolation order < 1.\n";
669 return 0.;
670 }
671
672 // Deal with the case that X is located at A(1) or A(N).
673 if (fabs(x - a[0]) < 1.e-6 * (fabs(a[0]) + fabs(a[nn - 1]))) {
674 return f[0];
675 }
676 if (fabs(x - a[nn - 1]) < 1.e-6 * (fabs(a[0]) + fabs(a[nn - 1]))) {
677 return f[nn - 1];
678 }
679
680 // Find subscript IX of X in array A.
681 int n = nn;
682 int m;
683 if (mm <= mmax && mm <= n - 1) {
684 m = mm;
685 } else {
686 if (mmax <= n - 1) {
687 m = mmax;
688 } else {
689 m = n - 1;
690 }
691 }
692 int mplus = m + 1;
693 int ix = 0;
694 int iy = n + 1;
695 int mid;
696 if (a[0] > a[n - 1]) {
697 // Search decreasing arguments.
698 do {
699 mid = (ix + iy) / 2;
700 if (x > a[mid - 1]) {
701 iy = mid;
702 } else {
703 ix = mid;
704 }
705 } while (iy - ix > 1);
706 } else {
707 // Search increasing arguments.
708 do {
709 mid = (ix + iy) / 2;
710 if (x < a[mid - 1]) {
711 iy = mid;
712 } else {
713 ix = mid;
714 }
715 } while (iy - ix > 1);
716 }
717 // Copy reordered interpolation points into (T[I],D[I]), setting
718 // EXTRA to True if M+2 points to be used.
719 int npts = m + 2 - (m % 2);
720 int ip = 0;
721 int l = 0;
722 int isub;
723 do {
724 isub = ix + l;
725 if ((1 > isub) || (isub > n)) {
726 // Skip point.
727 npts = mplus;
728 } else {
729 // Insert point.
730 ip++;
731 t[ip - 1] = a[isub - 1];
732 d[ip - 1] = f[isub - 1];
733 }
734 if (ip < npts) {
735 l = -l;
736 if (l >= 0) {
737 l++;
738 }
739 }
740 } while (ip < npts);
741
742 bool extra = npts != mplus;
743 // Replace d by the leading diagonal of a divided-difference table,
744 // supplemented by an extra line if EXTRA is True.
745 for (int l = 1; l <= m; l++) {
746 if (extra) {
747 isub = mplus - l;
748 d[m + 1] = (d[m + 1] - d[m - 1]) / (t[m + 1] - t[isub - 1]);
749 }
750 int i = mplus;
751 for (int j = l; j <= m; j++) {
752 isub = i - l;
753 d[i - 1] = (d[i - 1] - d[i - 1 - 1]) / (t[i - 1] - t[isub - 1]);
754 i--;
755 }
756 }
757 // Evaluate the Newton interpolation formula at X, averaging two values
758 // of last difference if EXTRA is True.
759 double sum = d[mplus - 1];
760 if (extra) {
761 sum = 0.5 * (sum + d[m + 1]);
762 }
763 int j = m;
764 for (int l = 1; l <= m; l++) {
765 sum = d[j - 1] + (x - t[j - 1]) * sum;
766 j--;
767 }
768 return sum;
769}

Referenced by Garfield::Sensor::ComputeThresholdCrossings(), and Garfield::Medium::Interpolate1D().

◆ GaussKronrod15()

double Garfield::Numerics::GaussKronrod15 ( double(*)(const double)  f,
const double  a,
const double  b 
)

Definition at line 599 of file Numerics.cc.

600 {
601
602 // Abscissae of the 15-point Kronrod rule
603 // xGK[1], xGK[3], ... abscissae of the 7-point Gauss rule
604 // xGK[0], xGK[2], ... abscissae which are optimally added
605 // to the 7-point Gauss rule
606 const double xGK[8] = {9.914553711208126e-01, 9.491079123427585e-01,
607 8.648644233597691e-01, 7.415311855993944e-01,
608 5.860872354676911e-01, 4.058451513773972e-01,
609 2.077849550078985e-01, 0.0e+00};
610 // Weights of the 15-point Kronrod rule
611 const double wGK[8] = {2.293532201052922e-02, 6.309209262997855e-02,
612 1.047900103222502e-01, 1.406532597155259e-01,
613 1.690047266392679e-01, 1.903505780647854e-01,
614 2.044329400752989e-01, 2.094821410847278e-01};
615 // Weights of the 7-point Gauss rule
616 const double wG[4] = {1.294849661688697e-01, 2.797053914892767e-01,
617 3.818300505051189e-01, 4.179591836734694e-01};
618
619 // Mid-point of the interval
620 const double center = 0.5 * (a + b);
621 // Half-length of the interval
622 const double halfLength = 0.5 * (b - a);
623
624 double fC = f(center);
625 // Result of the 7-point Gauss formula
626 double resG = fC * wG[3];
627 // Result of the 15-point Kronrod formula
628 double resK = fC * wGK[7];
629
630 for (int j = 0; j < 3; ++j) {
631 const int i = j * 2 + 1;
632 // Abscissa
633 const double x = halfLength * xGK[i];
634 // Function value
635 const double fSum = f(center - x) + f(center + x);
636 resG += wG[j] * fSum;
637 resK += wGK[i] * fSum;
638 }
639
640 for (int j = 0; j < 4; ++j) {
641 const int i = j * 2;
642 const double x = halfLength * xGK[i];
643 const double fSum = f(center - x) + f(center + x);
644 resK += wGK[i] * fSum;
645 }
646
647 return resK * halfLength;
648}