Garfield++ v2r0
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

Collection of numerical routines. More...

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)
 Numerical integration using 15-point Gauss-Kronrod algorithm.
 
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 (const std::vector< std::vector< std::vector< double > > > &value, const std::vector< double > &xAxis, const std::vector< double > &yAxis, const std::vector< double > &zAxis, const int nx, const int ny, const int nz, const double xx, const double yy, const double zz, double &f, const int iOrder)
 

Detailed Description

Collection of numerical routines.

Function Documentation

◆ BesselI0S()

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

Modified Bessel functions. Series expansions from Abramowitz and Stegun.

Definition at line 36 of file Numerics.hh.

36 {
37 return 1. + 3.5156229 * pow(xx / 3.75, 2) + 3.0899424 * pow(xx / 3.75, 4) +
38 1.2067492 * pow(xx / 3.75, 6) + 0.2659732 * pow(xx / 3.75, 8) +
39 0.0360768 * pow(xx / 3.75, 10) + 0.0045813 * pow(xx / 3.75, 12);
40}

Referenced by BesselK0S().

◆ BesselI1S()

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

Definition at line 42 of file Numerics.hh.

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

Referenced by BesselK1S().

◆ BesselK0L()

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

Definition at line 57 of file Numerics.hh.

57 {
58 return (exp(-xx) / sqrt(xx)) *
59 (1.25331414 - 0.07832358 * (2. / xx) + 0.02189568 * pow(2. / xx, 2) -
60 0.01062446 * pow(2. / xx, 3) + 0.00587872 * pow(2. / xx, 4) -
61 0.00251540 * pow(2. / xx, 5) + 0.00053208 * pow(2. / xx, 6));
62}

◆ BesselK0S()

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

Definition at line 50 of file Numerics.hh.

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

◆ BesselK1L()

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

Definition at line 72 of file Numerics.hh.

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

◆ BesselK1S()

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

Definition at line 64 of file Numerics.hh.

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

◆ Boxin3()

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

Definition at line 770 of file Numerics.cc.

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

Referenced by Garfield::Medium::CloneTable(), Garfield::Medium::CloneTensor(), Garfield::Medium::ElectronAttachment(), Garfield::Medium::ElectronDiffusion(), Garfield::Medium::ElectronLorentzAngle(), 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 do {
723 const int isub = ix + l;
724 if ((1 > isub) || (isub > n)) {
725 // Skip point.
726 npts = mplus;
727 } else {
728 // Insert point.
729 ip++;
730 t[ip - 1] = a[isub - 1];
731 d[ip - 1] = f[isub - 1];
732 }
733 if (ip < npts) {
734 l = -l;
735 if (l >= 0) {
736 l++;
737 }
738 }
739 } while (ip < npts);
740
741 bool extra = npts != mplus;
742 // Replace d by the leading diagonal of a divided-difference table,
743 // supplemented by an extra line if EXTRA is True.
744 for (l = 1; l <= m; l++) {
745 if (extra) {
746 const int isub = mplus - l;
747 d[m + 1] = (d[m + 1] - d[m - 1]) / (t[m + 1] - t[isub - 1]);
748 }
749 int i = mplus;
750 for (int j = l; j <= m; j++) {
751 const int isub = i - l;
752 d[i - 1] = (d[i - 1] - d[i - 1 - 1]) / (t[i - 1] - t[isub - 1]);
753 i--;
754 }
755 }
756 // Evaluate the Newton interpolation formula at X, averaging two values
757 // of last difference if EXTRA is True.
758 double sum = d[mplus - 1];
759 if (extra) {
760 sum = 0.5 * (sum + d[m + 1]);
761 }
762 int j = m;
763 for (l = 1; l <= m; l++) {
764 sum = d[j - 1] + (x - t[j - 1]) * sum;
765 j--;
766 }
767 return sum;
768}

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

◆ GaussKronrod15()

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

Numerical integration using 15-point Gauss-Kronrod algorithm.

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}