15 m_chargeCheck =
false;
24 std::cerr <<
m_className <<
"::GetVoltageRange:\n Unable to return "
25 <<
"voltage range (could not set up the cell).\n";
36 double& x1,
double& y1,
44 if (!m_cellset)
return false;
56 const double x1,
const double y1,
58 double& xc,
double& yc,
double& zc) {
64 if (m_w.empty())
return false;
66 const double dx = x1 - x0;
67 const double dy = y1 - y0;
68 const double d2 = dx * dx + dy * dy;
70 if (d2 < Small)
return false;
73 if ((m_perx && fabs(dx) >= m_sx) || (m_pery && fabs(dy) >= m_sy)) {
75 <<
" Particle crossed more than one period.\n";
83 const double xm = 0.5 * (x0 + x1);
84 const double ym = 0.5 * (y0 + y1);
86 for (
unsigned int i = 0; i < m_nWires; ++i) {
87 double xw = m_w[i].x, yw = m_w[i].y;
89 xw += m_sx * int(round((xm - xw) / m_sx));
92 yw += m_sy * int(round((ym - yw) / m_sy));
95 const double xIn0 = dx * (xw - x0) + dy * (yw - y0);
97 if (xIn0 < 0.)
continue;
98 const double xIn1 = -(dx * (xw - x1) + dy * (yw - y1));
100 if (xIn1 < 0.)
continue;
102 const double xw0 = xw - x0;
103 const double xw1 = xw - x1;
104 const double yw0 = yw - y0;
105 const double yw1 = yw - y1;
106 const double dw02 = xw0 * xw0 + yw0 * yw0;
107 const double dw12 = xw1 * xw1 + yw1 * yw1;
108 if (xIn1 * xIn1 * dw02 > xIn0 * xIn0 * dw12) {
109 dMin2 = dw02 - xIn0 * xIn0 / d2;
111 dMin2 = dw12 - xIn1 * xIn1 / d2;
114 const double r2 = 0.25 * m_w[i].d * m_w[i].d;
118 const double p = -xIn0 / d2;
119 const double q = (dw02 - r2) / d2;
120 const double t1 = -p + sqrt(p * p - q);
121 const double t2 = -p - sqrt(p * p - q);
122 const double t = std::min(t1, t2);
125 zc = z0 + t * (z1 - z0);
133 const double yin,
const double zin,
134 double& xw,
double& yw,
140 int nX = 0, nY = 0, nPhi = 0;
142 nX = int(round(xin / m_sx));
145 if (m_pery && m_tube) {
146 Cartesian2Polar(xin, yin, x0, y0);
147 nPhi = int(round(DegreeToRad * y0 / m_sy));
148 y0 -= RadToDegree * m_sy * nPhi;
149 Polar2Cartesian(x0, y0, x0, y0);
151 nY = int(round(yin / m_sy));
156 if (m_perx && m_ynplan[0] && x0 <= m_coplan[0]) x0 += m_sx;
157 if (m_perx && m_ynplan[1] && x0 >= m_coplan[1]) x0 -= m_sx;
158 if (m_pery && m_ynplan[2] && y0 <= m_coplan[2]) y0 += m_sy;
159 if (m_pery && m_ynplan[3] && y0 >= m_coplan[3]) y0 -= m_sy;
161 for (
unsigned int i = 0; i < m_nWires; ++i) {
163 if (qin * m_w[i].e > 0.)
continue;
164 const double dxw0 = m_w[i].x - x0;
165 const double dyw0 = m_w[i].y - y0;
166 const double r2 = dxw0 * dxw0 + dyw0 * dyw0;
167 const double rTrap = 0.5 * m_w[i].d * m_w[i].nTrap;
168 if (r2 < rTrap * rTrap) {
172 if (m_perx && m_ynplan[0] && x0 <= m_coplan[0]) x0 -= m_sx;
173 if (m_perx && m_ynplan[1] && x0 >= m_coplan[1]) x0 += m_sx;
174 if (m_pery && m_ynplan[2] && y0 <= m_coplan[2]) y0 -= m_sy;
175 if (m_pery && m_ynplan[3] && y0 >= m_coplan[3]) y0 += m_sy;
176 if (m_pery && m_tube) {
178 Cartesian2Polar(xw, yw, rhow, phiw);
179 phiw += RadToDegree * m_sy * nPhi;
180 Polar2Cartesian(rhow, phiw, xw, yw);
184 if (m_perx) xw += m_sx * nX;
187 std::cout <<
" (" << xin <<
", " << yin <<
", " << zin <<
")"
188 <<
" within trap radius of wire " << i <<
".\n";
198 const double diameter,
199 const double voltage,
200 const std::string& label,
201 const double length,
const double tension,
202 double rho,
const int ntrap) {
205 if (diameter <= 0.) {
207 std::cerr <<
" Unphysical wire diameter.\n";
213 std::cerr <<
" Unphysical wire tension.\n";
219 std::cerr <<
" Unphysical wire density.\n";
225 std::cerr <<
" Unphysical wire length.\n";
231 std::cerr <<
" Number of trap radii must be > 0.\n";
238 newWire.d = diameter;
241 newWire.type = label;
244 newWire.nTrap = ntrap;
246 m_w.push_back(newWire);
256 const std::string& label) {
261 <<
" Unphysical tube dimension.\n";
265 if (nEdges < 3 && nEdges != 0) {
267 <<
" Unphysical number of tube edges (" << nEdges <<
")\n";
274 <<
" Warning: Existing tube settings will be overwritten.\n";
283 m_cotube2 = radius * radius;
288 planes[4].type = label;
297 const std::string& lab) {
299 if (m_ynplan[0] && m_ynplan[1]) {
301 std::cerr <<
" There are already two x planes defined.\n";
309 planes[1].type = lab;
315 planes[0].type = lab;
325 const std::string& lab) {
327 if (m_ynplan[2] && m_ynplan[3]) {
329 std::cerr <<
" There are already two y planes defined.\n";
337 planes[3].type = lab;
343 planes[2].type = lab;
353 const double x,
const double smin,
355 const std::string& label,
358 if (!m_ynplan[0] && !m_ynplan[1]) {
359 std::cerr <<
m_className <<
"::AddStripOnPlaneX:\n";
360 std::cerr <<
" There are no planes at constant x defined.\n";
364 if (direction !=
'y' && direction !=
'Y' && direction !=
'z' &&
366 std::cerr <<
m_className <<
"::AddStripOnPlaneX:\n";
367 std::cerr <<
" Invalid direction (" << direction <<
").\n";
368 std::cerr <<
" Only strips in y or z direction are possible.\n";
372 if (fabs(smax - smin) < Small) {
373 std::cerr <<
m_className <<
"::AddStripOnPlaneX:\n";
374 std::cerr <<
" Strip width must be greater than zero.\n";
379 newStrip.type = label;
381 newStrip.smin = std::min(smin, smax);
382 newStrip.smax = std::max(smin, smax);
391 const double d0 = fabs(m_coplan[0] - x);
392 const double d1 = fabs(m_coplan[1] - x);
393 if (d1 < d0) iplane = 1;
396 if (direction ==
'y' || direction ==
'Y') {
397 planes[iplane].strips1.push_back(newStrip);
399 planes[iplane].strips2.push_back(newStrip);
404 const double y,
const double smin,
406 const std::string& label,
409 if (!m_ynplan[2] && !m_ynplan[3]) {
410 std::cerr <<
m_className <<
"::AddStripOnPlaneY:\n";
411 std::cerr <<
" There are no planes at constant y defined.\n";
415 if (direction !=
'x' && direction !=
'X' && direction !=
'z' &&
417 std::cerr <<
m_className <<
"::AddStripOnPlaneY:\n";
418 std::cerr <<
" Invalid direction (" << direction <<
").\n";
419 std::cerr <<
" Only strips in x or z direction are possible.\n";
423 if (fabs(smax - smin) < Small) {
424 std::cerr <<
m_className <<
"::AddStripOnPlaneY:\n";
425 std::cerr <<
" Strip width must be greater than zero.\n";
430 newStrip.type = label;
432 newStrip.smin = std::min(smin, smax);
433 newStrip.smax = std::max(smin, smax);
442 const double d2 = fabs(m_coplan[2] - y);
443 const double d3 = fabs(m_coplan[3] - y);
444 if (d3 < d2) iplane = 3;
447 if (direction ==
'x' || direction ==
'X') {
448 planes[iplane].strips1.push_back(newStrip);
450 planes[iplane].strips2.push_back(newStrip);
455 const double x,
const double ymin,
const double ymax,
const double zmin,
456 const double zmax,
const std::string& label,
const double gap) {
458 if (!m_ynplan[0] && !m_ynplan[1]) {
459 std::cerr <<
m_className <<
"::AddPixelOnPlaneX:\n";
460 std::cerr <<
" There are no planes at constant x defined.\n";
464 if (fabs(ymax - ymin) < Small || fabs(zmax - zmin) < Small) {
465 std::cerr <<
m_className <<
"::AddSPixelOnPlaneX:\n";
466 std::cerr <<
" Pixel width must be greater than zero.\n";
471 newPixel.type = label;
473 newPixel.smin = std::min(ymin, ymax);
474 newPixel.smax = std::max(ymin, ymax);
475 newPixel.zmin = std::min(zmin, zmax);
476 newPixel.zmax = std::max(zmin, zmax);
485 const double d0 = fabs(m_coplan[0] - x);
486 const double d1 = fabs(m_coplan[1] - x);
487 if (d1 < d0) iplane = 1;
490 planes[iplane].pixels.push_back(newPixel);
494 const double y,
const double xmin,
const double xmax,
const double zmin,
495 const double zmax,
const std::string& label,
const double gap) {
497 if (!m_ynplan[2] && !m_ynplan[3]) {
498 std::cerr <<
m_className <<
"::AddPixelOnPlaneY:\n";
499 std::cerr <<
" There are no planes at constant y defined.\n";
503 if (fabs(xmax - xmin) < Small || fabs(zmax - zmin) < Small) {
504 std::cerr <<
m_className <<
"::AddPixelOnPlaneY:\n";
505 std::cerr <<
" Pixel width must be greater than zero.\n";
510 newPixel.type = label;
512 newPixel.smin = std::min(xmin, xmax);
513 newPixel.smax = std::max(xmin, xmax);
514 newPixel.zmin = std::min(zmin, zmax);
515 newPixel.zmax = std::max(zmin, zmax);
524 const double d0 = fabs(m_coplan[2] - y);
525 const double d1 = fabs(m_coplan[3] - y);
526 if (d1 < d0) iplane = 3;
529 planes[iplane].pixels.push_back(newPixel);
535 std::cerr <<
m_className <<
"::SetPeriodicityX:\n";
536 std::cerr <<
" Periodic length must be greater than zero.\n";
548 std::cerr <<
m_className <<
"::SetPeriodicityY:\n";
549 std::cerr <<
" Periodic length must be greater than zero.\n";
580void ComponentAnalyticField::UpdatePeriodicity() {
589 std::cerr <<
m_className <<
"::UpdatePeriodicity:\n";
590 std::cerr <<
" Periodicity in x direction was enabled"
591 <<
" but periodic length is not set.\n";
605 std::cerr <<
m_className <<
"::UpdatePeriodicity:\n";
606 std::cerr <<
" Periodicity in y direction was enabled"
607 <<
" but periodic length is not set.\n";
617 std::cerr <<
m_className <<
"::UpdatePeriodicity:\n";
618 std::cerr <<
" Periodicity in z is not possible.\n";
622 std::cerr <<
m_className <<
"::UpdatePeriodicity:\n";
623 std::cerr <<
" Mirror periodicity is not possible.\n";
627 std::cerr <<
m_className <<
"::UpdatePeriodicity:\n";
628 std::cerr <<
" Axial periodicity is not possible.\n";
632 std::cerr <<
m_className <<
"::UpdatePeriodicity:\n";
633 std::cerr <<
" Rotation symmetry is not possible.\n";
638 const double z,
const double q) {
645 newCharge.e = q / FourPiEpsilon0;
646 m_ch3d.push_back(newCharge);
659 if (m_ch3d.empty()) {
660 std::cout <<
" No charges present.\n";
663 std::cout <<
" x [cm] y [cm] z [cm] charge [fC]\n";
664 const unsigned int n3d = m_ch3d.size();
665 for (
unsigned int i = 0; i < n3d; ++i) {
666 std::cout <<
" " << std::setw(9) << m_ch3d[i].x <<
" " << std::setw(9)
667 << m_ch3d[i].y <<
" " << std::setw(9) << m_ch3d[i].z <<
" "
668 << std::setw(11) << m_ch3d[i].e * FourPiEpsilon0 <<
"\n";
674 if (m_ynplan[0] && m_ynplan[1]) {
676 }
else if (m_ynplan[0] || m_ynplan[1]) {
684 if (m_ynplan[2] && m_ynplan[3]) {
686 }
else if (m_ynplan[2] || m_ynplan[3]) {
693 double& diameter,
double& voltage,
694 std::string& label,
double& length,
695 double& charge,
int& ntrap)
const {
699 std::cerr <<
" Wire index is out of range.\n";
710 ntrap = m_w[i].nTrap;
715 double& x,
double& voltage,
716 std::string& label)
const {
718 if (i >= 2 || (i == 1 && !m_ynplan[1])) {
720 std::cerr <<
" Plane index is out of range.\n";
725 voltage = m_vtplan[i];
726 label = planes[i].type;
731 double& y,
double& voltage,
732 std::string& label)
const {
734 if (i >= 2 || (i == 1 && !m_ynplan[3])) {
736 std::cerr <<
" Plane index is out of range.\n";
741 voltage = m_vtplan[i + 2];
742 label = planes[i + 2].type;
747 std::string& label)
const {
749 if (!m_tube)
return false;
753 label = planes[4].type;
757int ComponentAnalyticField::Field(
const double xin,
const double yin,
758 const double zin,
double& ex,
double& ey,
759 double& ez,
double& volt,
const bool opt) {
780 ex = ey = ez = volt = 0.;
784 if (!Prepare())
return -11;
787 double xpos = xin, ypos = yin;
791 xpos -= m_sx * int(round(xin / m_sx));
794 if (m_pery && m_tube) {
795 Cartesian2Polar(xin, yin, xpos, ypos);
796 arot = RadToDegree * m_sy * int(round(DegreeToRad * ypos / m_sy));
798 Polar2Cartesian(xpos, ypos, xpos, ypos);
800 ypos -= m_sy * int(round(yin / m_sy));
804 if (m_perx && m_ynplan[0] && xpos <= m_coplan[0]) xpos += m_sx;
805 if (m_perx && m_ynplan[1] && xpos >= m_coplan[1]) xpos -= m_sx;
806 if (m_pery && m_ynplan[2] && ypos <= m_coplan[2]) ypos += m_sy;
807 if (m_pery && m_ynplan[3] && ypos >= m_coplan[3]) ypos -= m_sy;
811 if (!InTube(xpos, ypos, m_cotube, m_ntube)) {
816 if (m_ynplan[0] && xpos < m_coplan[0]) {
820 if (m_ynplan[1] && xpos > m_coplan[1]) {
824 if (m_ynplan[2] && ypos < m_coplan[2]) {
828 if (m_ynplan[3] && ypos > m_coplan[3]) {
835 for (
int i = m_nWires; i--;) {
836 double dx = xpos - m_w[i].x;
837 double dy = ypos - m_w[i].y;
839 if (m_perx) dx -= m_sx * int(round(dx / m_sx));
840 if (m_pery) dy -= m_sy * int(round(dy / m_sy));
842 if (dx * dx + dy * dy < 0.25 * m_w[i].d * m_w[i].d) {
849 switch (m_cellType) {
851 FieldA00(xpos, ypos, ex, ey, volt, opt);
854 FieldB1X(xpos, ypos, ex, ey, volt, opt);
857 FieldB1Y(xpos, ypos, ex, ey, volt, opt);
860 FieldB2X(xpos, ypos, ex, ey, volt, opt);
863 FieldB2Y(xpos, ypos, ex, ey, volt, opt);
866 FieldC10(xpos, ypos, ex, ey, volt, opt);
869 FieldC2X(xpos, ypos, ex, ey, volt, opt);
872 FieldC2Y(xpos, ypos, ex, ey, volt, opt);
875 FieldC30(xpos, ypos, ex, ey, volt, opt);
878 FieldD10(xpos, ypos, ex, ey, volt, opt);
881 FieldD20(xpos, ypos, ex, ey, volt, opt);
884 FieldD30(xpos, ypos, ex, ey, volt, opt);
889 std::cerr <<
" Unknown cell type (id " << m_cellType <<
")\n";
896 double exd = 0., eyd = 0., voltd = 0.;
897 switch (m_cellType) {
922 if (m_pery && m_tube) {
924 Cartesian2Polar(ex, ey, xaux, yaux);
926 Polar2Cartesian(xaux, yaux, ex, ey);
932 volt += m_corvta * xpos + m_corvtb * ypos + m_corvtc;
935 if (!m_ch3d.empty()) {
936 double ex3d = 0., ey3d = 0., ez3d = 0., volt3d = 0.;
937 switch (m_cellType) {
941 Field3dA00(xin, yin, zin, ex3d, ey3d, ez3d, volt3d);
944 Field3dB2X(xin, yin, zin, ex3d, ey3d, ez3d, volt3d);
947 Field3dB2Y(xin, yin, zin, ex3d, ey3d, ez3d, volt3d);
950 Field3dD10(xin, yin, zin, ex3d, ey3d, ez3d, volt3d);
953 Field3dA00(xin, yin, zin, ex3d, ey3d, ez3d, volt3d);
965void ComponentAnalyticField::CellInit() {
978 m_xmin = m_xmax = 0.;
979 m_ymin = m_ymax = 0.;
980 m_zmin = m_zmax = 0.;
984 m_perx = m_pery =
false;
990 m_scellTypeFourier =
"A ";
991 fperx = fpery =
false;
992 mxmin = mxmax = mymin = mymax = 0;
1016 m_zmult = std::complex<double>(0., 0.);
1017 m_p1 = m_p2 = m_c1 = 0.;
1026 m_corvta = m_corvtb = m_corvtc = 0.;
1031 for (
int i = 0; i < 4; ++i) {
1032 m_ynplan[i] =
false;
1037 m_ynplax = m_ynplay =
false;
1038 m_coplax = m_coplay = 1.;
1040 for (
int i = 0; i < 5; ++i) {
1041 planes[i].type =
'?';
1043 planes[i].ewxcor = 0.;
1044 planes[i].ewycor = 0.;
1045 planes[i].strips1.clear();
1046 planes[i].strips2.clear();
1047 planes[i].pixels.clear();
1069 down[0] = down[1] = 0.;
1073bool ComponentAnalyticField::Prepare() {
1078 <<
" The cell does not meet the requirements.\n";
1086 std::cerr <<
" Type identification of the cell failed.\n";
1091 std::cout <<
" Cell is of type " << m_scellType <<
".\n";
1097 std::cerr <<
" Calculation of charges failed.\n";
1102 std::cout <<
" Calculation of charges was successful.\n";
1106 if (!PrepareStrips()) {
1108 std::cerr <<
" Strip/pixel preparation failed.\n";
1116bool ComponentAnalyticField::CellCheck() {
1131 double conew1 = m_coplan[0] - m_sx * int(round(m_coplan[0] / m_sx));
1132 double conew2 = m_coplan[1] - m_sx * int(round(m_coplan[1] / m_sx));
1134 if (m_ynplan[0] && m_ynplan[1] && conew1 == conew2) {
1141 if ((conew1 != m_coplan[0] && m_ynplan[0]) ||
1142 (conew2 != m_coplan[1] && m_ynplan[1])) {
1144 std::cout <<
" The planes in x or r are moved to the basic period.\n";
1145 std::cout <<
" This should not affect the results.";
1147 m_coplan[0] = conew1;
1148 m_coplan[1] = conew2;
1151 if (m_ynplan[0] && m_ynplan[1] &&
fabs(m_coplan[1] - m_coplan[0]) != m_sx) {
1153 std::cerr <<
" The separation of the x or r planes"
1154 <<
" does not match the period.\b";
1155 std::cerr <<
" The periodicity is cancelled.\n";
1159 if (m_ynplan[0] && m_ynplan[1] && m_vtplan[0] != m_vtplan[1]) {
1161 std::cerr <<
" The voltages of the two x (or r) planes differ.\n";
1162 std::cerr <<
" The periodicity is cancelled.\n";
1169 double conew3 = m_coplan[2] - m_sy * int(round(m_coplan[2] / m_sy));
1170 double conew4 = m_coplan[3] - m_sy * int(round(m_coplan[3] / m_sy));
1172 if (m_ynplan[2] && m_ynplan[3] && conew3 == conew4) {
1179 if ((conew3 != m_coplan[2] && m_ynplan[2]) ||
1180 (conew4 != m_coplan[3] && m_ynplan[3])) {
1182 std::cout <<
" The planes in y are moved to the basic period.\n";
1183 std::cout <<
" This should not affect the results.";
1185 m_coplan[2] = conew3;
1186 m_coplan[3] = conew4;
1189 if (m_ynplan[2] && m_ynplan[3] &&
fabs(m_coplan[3] - m_coplan[2]) != m_sy) {
1191 std::cerr <<
" The separation of the two y planes"
1192 <<
" does not match the period.\b";
1193 std::cerr <<
" The periodicity is cancelled.\n";
1197 if (m_ynplan[2] && m_ynplan[3] && m_vtplan[2] != m_vtplan[3]) {
1199 std::cerr <<
" The voltages of the two y planes differ.\n";
1200 std::cerr <<
" The periodicity is cancelled.\n";
1206 for (
int i = 0; i < 2; ++i) {
1207 for (
int j = 2; j < 3; ++j) {
1208 if (m_ynplan[i] && m_ynplan[j] && m_vtplan[i] != m_vtplan[j]) {
1210 std::cerr <<
" Conflicting potential of 2 crossing planes.\n";
1211 std::cerr <<
" One y (or phi) plane is removed.\n";
1212 m_ynplan[j] =
false;
1218 for (
int i = 0; i < 3; i += 2) {
1219 if (m_ynplan[i] && m_ynplan[i + 1]) {
1220 if (m_coplan[i] == m_coplan[i + 1]) {
1222 std::cerr <<
" Two planes are on top of each other.\n";
1223 std::cerr <<
" One of them is removed.\n";
1224 m_ynplan[i + 1] =
false;
1226 if (m_coplan[i] > m_coplan[i + 1]) {
1229 std::cout <<
" Planes " << i <<
" and " << i + 1
1230 <<
" are interchanged.\n";
1233 const double cohlp = m_coplan[i];
1234 m_coplan[i] = m_coplan[i + 1];
1235 m_coplan[i + 1] = cohlp;
1237 const double vthlp = m_vtplan[i];
1238 m_vtplan[i] = m_vtplan[i + 1];
1239 m_vtplan[i + 1] = vthlp;
1241 plane plahlp = planes[i];
1242 planes[i] = planes[i + 1];
1243 planes[i + 1] = plahlp;
1250 for (
unsigned int i = 0; i < m_nWires; ++i) {
1251 const double xnew = m_w[i].x - m_sx * int(round(m_w[i].x / m_sx));
1252 if (
int(round(m_w[i].x / m_sx)) != 0) {
1253 double xprt = m_w[i].x;
1254 double yprt = m_w[i].y;
1255 if (m_polar) RTheta2RhoPhi(m_w[i].x, m_w[i].y, xprt, yprt);
1256 std::cout <<
m_className <<
"::CellCheck:\n The " << m_w[i].type
1257 <<
"-wire at (" << xprt <<
", " << yprt
1258 <<
") is moved to the basic x (or r) period.\n"
1259 <<
" This should not affect the results.\n";
1266 if (m_tube && m_pery) {
1267 for (
unsigned int i = 0; i < m_nWires; ++i) {
1268 double xnew = m_w[i].x;
1269 double ynew = m_w[i].y;
1270 Cartesian2Polar(xnew, ynew, xnew, ynew);
1271 if (
int(round(DegreeToRad * ynew / m_sy)) != 0) {
1273 std::cout <<
" The " << m_w[i].type <<
"-wire at (" << m_w[i].x
1275 <<
") is moved to the basic phi period.\n";
1276 std::cout <<
" This should not affect the results.\n";
1277 ynew -= RadToDegree * m_sy * int(round(DegreeToRad * ynew / m_sy));
1278 Polar2Cartesian(xnew, ynew, m_w[i].x, m_w[i].y);
1281 }
else if (m_pery) {
1282 for (
unsigned int i = 0; i < m_nWires; ++i) {
1283 double ynew = m_w[i].y - m_sy * int(round(m_w[i].y / m_sy));
1284 if (
int(round(m_w[i].y / m_sy)) != 0) {
1285 double xprt = m_w[i].x;
1286 double yprt = m_w[i].y;
1287 if (m_polar) RTheta2RhoPhi(m_w[i].x, m_w[i].y, xprt, yprt);
1288 std::cout <<
m_className <<
"::CellCheck:\n The " << m_w[i].type
1289 <<
"-wire at (" << xprt <<
", " << yprt
1290 <<
") is moved to the basic y period.\n"
1291 <<
" This should not affect the results.\n";
1298 int iplan1 = 0, iplan2 = 0, iplan3 = 0, iplan4 = 0;
1299 for (
unsigned int i = 0; i < m_nWires; ++i) {
1300 if (m_ynplan[0] && m_w[i].x <= m_coplan[0]) ++iplan1;
1301 if (m_ynplan[1] && m_w[i].x <= m_coplan[1]) ++iplan2;
1302 if (m_ynplan[2] && m_w[i].y <= m_coplan[2]) ++iplan3;
1303 if (m_ynplan[3] && m_w[i].y <= m_coplan[3]) ++iplan4;
1307 if (m_ynplan[0] && m_ynplan[1]) {
1308 if (iplan1 >
int(m_nWires) / 2) {
1309 m_ynplan[1] =
false;
1314 if (iplan2 <
int(m_nWires) / 2) {
1315 m_ynplan[0] =
false;
1321 if (m_ynplan[0] && !m_ynplan[1]) {
1322 if (iplan1 >
int(m_nWires) / 2)
1327 if (m_ynplan[1] && !m_ynplan[0]) {
1328 if (iplan2 <
int(m_nWires) / 2)
1334 if (m_ynplan[2] && m_ynplan[3]) {
1335 if (iplan3 >
int(m_nWires) / 2) {
1336 m_ynplan[3] =
false;
1341 if (iplan4 <
int(m_nWires) / 2) {
1342 m_ynplan[2] =
false;
1348 if (m_ynplan[2] && !m_ynplan[3]) {
1349 if (iplan3 >
int(m_nWires) / 2)
1354 if (m_ynplan[3] && !m_ynplan[2]) {
1355 if (iplan4 <
int(m_nWires) / 2)
1363 m_ynplan[0] =
false;
1365 m_coplan[1] = m_coplan[0];
1366 m_vtplan[1] = m_vtplan[0];
1367 planes[1] = planes[0];
1371 m_ynplan[1] =
false;
1373 m_coplan[0] = m_coplan[1];
1374 m_vtplan[0] = m_vtplan[1];
1375 planes[0] = planes[1];
1379 m_ynplan[2] =
false;
1381 m_coplan[3] = m_coplan[2];
1382 m_vtplan[3] = m_vtplan[2];
1383 planes[3] = planes[2];
1387 m_ynplan[3] =
false;
1389 m_coplan[2] = m_coplan[3];
1390 m_vtplan[2] = m_vtplan[3];
1391 planes[2] = planes[3];
1394 std::vector<bool> wrong(m_nWires,
false);
1396 for (
unsigned int i = 0; i < m_nWires; ++i) {
1397 const double rw = 0.5 * m_w[i].d;
1398 if (m_ynplan[0] && m_w[i].x - rw <= m_coplan[0]) wrong[i] =
true;
1399 if (m_ynplan[1] && m_w[i].x + rw >= m_coplan[1]) wrong[i] =
true;
1400 if (m_ynplan[2] && m_w[i].y - rw <= m_coplan[2]) wrong[i] =
true;
1401 if (m_ynplan[3] && m_w[i].y + rw >= m_coplan[3]) wrong[i] =
true;
1403 if (!InTube(m_w[i].x, m_w[i].y, m_cotube, m_ntube)) {
1405 std::cerr <<
" The " << m_w[i].type <<
"-wire at (" << m_w[i].x
1406 <<
", " << m_w[i].y <<
") is located outside the tube.\n";
1407 std::cerr <<
" This wire is removed.\n";
1410 }
else if (wrong[i]) {
1411 double xprt = m_w[i].x;
1412 double yprt = m_w[i].y;
1413 if (m_polar) RTheta2RhoPhi(m_w[i].x, m_w[i].y, xprt, yprt);
1414 std::cerr <<
m_className <<
"::CellCheck:\n The " << m_w[i].type
1415 <<
"-wire at (" << xprt <<
", " << yprt <<
") is located "
1416 <<
"outside the planes.\n This wire is removed.\n";
1417 }
else if ((m_perx && m_w[i].d >= m_sx) || (m_pery && m_w[i].d >= m_sy)) {
1418 double xprt = m_w[i].x;
1419 double yprt = m_w[i].y;
1420 if (m_polar) RTheta2RhoPhi(m_w[i].x, m_w[i].y, xprt, yprt);
1421 std::cerr <<
m_className <<
"::CellCheck:\n The diameter of the "
1422 << m_w[i].type <<
"-wire at (" << xprt <<
", " << yprt
1423 <<
") exceeds 1 period.\n This wire is removed.\n";
1429 for (
unsigned int i = 0; i < m_nWires; ++i) {
1430 if (wrong[i])
continue;
1431 for (
unsigned int j = i + 1; j < m_nWires; ++j) {
1432 if (wrong[j])
continue;
1437 double xaux1, xaux2, yaux1, yaux2;
1438 Cartesian2Polar(m_w[i].x, m_w[i].y, xaux1, yaux1);
1439 Cartesian2Polar(m_w[j].x, m_w[j].y, xaux2, yaux2);
1440 yaux1 -= m_sy * int(round(yaux1 / m_sy));
1441 yaux2 -= m_sy * int(round(yaux2 / m_sy));
1442 Polar2Cartesian(xaux1, yaux1, xaux1, yaux1);
1443 Polar2Cartesian(xaux2, yaux2, xaux2, yaux2);
1444 xsepar = xaux1 - xaux2;
1445 ysepar = yaux1 - yaux2;
1447 xsepar = m_w[i].x - m_w[j].x;
1448 ysepar = m_w[i].y - m_w[j].y;
1451 xsepar =
fabs(m_w[i].x - m_w[j].x);
1452 if (m_perx) xsepar -= m_sx * int(round(xsepar / m_sx));
1453 ysepar =
fabs(m_w[i].y - m_w[j].y);
1454 if (m_pery) ysepar -= m_sy * int(round(ysepar / m_sy));
1456 const double dij = m_w[i].d + m_w[j].d;
1457 if (xsepar * xsepar + ysepar * ysepar > 0.25 * dij * dij)
continue;
1458 double xprti = m_w[i].x;
1459 double yprti = m_w[i].y;
1460 double xprtj = m_w[j].x;
1461 double yprtj = m_w[j].y;
1463 RTheta2RhoPhi(xprti, yprti, xprti, yprti);
1464 RTheta2RhoPhi(xprtj, yprtj, xprtj, yprtj);
1466 std::cerr <<
m_className <<
"::CellCheck:\n Wires "
1467 << m_w[i].type <<
" at (" << xprti <<
", " << yprti <<
") and "
1468 << m_w[j].type <<
" at (" << xprtj <<
", " << yprtj
1469 <<
") overlap at least partially.\n"
1470 <<
" The latter wire is removed.\n";
1476 const int iWires = m_nWires;
1478 for (
int i = 0; i < iWires; ++i) {
1480 m_w[m_nWires].x = m_w[i].x;
1481 m_w[m_nWires].y = m_w[i].y;
1482 m_w[m_nWires].d = m_w[i].d;
1483 m_w[m_nWires].v = m_w[i].v;
1484 m_w[m_nWires].type = m_w[i].type;
1485 m_w[m_nWires].u = m_w[i].u;
1486 m_w[m_nWires].e = m_w[i].e;
1487 m_w[m_nWires].ind = m_w[i].ind;
1488 m_w[m_nWires].nTrap = m_w[i].nTrap;
1494 int nElements = m_nWires;
1495 if (m_ynplan[0]) ++nElements;
1496 if (m_ynplan[1]) ++nElements;
1497 if (m_ynplan[2]) ++nElements;
1498 if (m_ynplan[3]) ++nElements;
1499 if (m_tube) ++nElements;
1501 if (nElements < 2) {
1503 std::cerr <<
" At least 2 elements are necessary.\n";
1504 std::cerr <<
" Cell rejected.\n";
1514 m_xmin = m_xmax = 0.;
1515 m_ymin = m_ymax = 0.;
1516 m_zmin = m_zmax = 0.;
1520 for (
int i = m_nWires; i--;) {
1521 const double rw = 0.5 * m_w[i].d;
1523 m_xmin = std::min(m_xmin, m_w[i].x - rw);
1524 m_xmax = std::max(m_xmax, m_w[i].x + rw);
1526 m_xmin = m_w[i].x - rw;
1527 m_xmax = m_w[i].x + rw;
1531 m_ymin = std::min(m_ymin, m_w[i].y - rw);
1532 m_ymax = std::max(m_ymax, m_w[i].y + rw);
1534 m_ymin = m_w[i].y - rw;
1535 m_ymax = m_w[i].y + rw;
1539 m_zmin = std::min(m_zmin, -0.5 * m_w[i].u);
1540 m_zmax = std::max(m_zmax, +0.5 * m_w[i].u);
1542 m_zmin = -0.5 * m_w[i].u;
1543 m_zmax = +0.5 * m_w[i].u;
1547 vmin = std::min(vmin, m_w[i].v);
1548 vmax = std::max(vmax, m_w[i].v);
1550 vmin = vmax = m_w[i].v;
1555 for (
int i = 0; i < 4; ++i) {
1556 if (!m_ynplan[i])
continue;
1559 m_xmin = std::min(m_xmin, m_coplan[i]);
1560 m_xmax = std::max(m_xmax, m_coplan[i]);
1562 m_xmin = m_xmax = m_coplan[i];
1567 m_ymin = std::min(m_ymin, m_coplan[i]);
1568 m_ymax = std::max(m_ymax, m_coplan[i]);
1570 m_ymin = m_ymax = m_coplan[i];
1575 vmin = std::min(vmin, m_vtplan[i]);
1576 vmax = std::max(vmax, m_vtplan[i]);
1578 vmin = vmax = m_vtplan[i];
1585 m_xmin = -1.1 * m_cotube;
1586 m_xmax = +1.1 * m_cotube;
1588 m_ymin = -1.1 * m_cotube;
1589 m_ymax = +1.1 * m_cotube;
1591 vmin = std::min(vmin, m_vttube);
1592 vmax = std::max(vmax, m_vttube);
1597 if (m_perx && m_sx > (m_xmax - m_xmin)) {
1598 m_xmin = -0.5 * m_sx;
1599 m_xmax = +0.5 * m_sx;
1603 if (m_pery && m_sy > (m_ymax - m_ymin)) {
1604 m_ymin = -0.5 * m_sy;
1605 m_ymax = +0.5 * m_sy;
1609 if (m_polar && (m_ymax - m_ymin) >= TwoPi) {
1616 if (setx && m_xmin != m_xmax && (m_ymin == m_ymax || !sety)) {
1617 m_ymin -= 0.5 *
fabs(m_xmax - m_xmin);
1618 m_ymax += 0.5 *
fabs(m_xmax - m_xmin);
1621 if (sety && m_ymin != m_ymax && (m_xmin == m_xmax || !setx)) {
1622 m_xmin -= 0.5 *
fabs(m_ymax - m_ymin);
1623 m_xmax += 0.5 *
fabs(m_ymax - m_ymin);
1628 m_zmin = -0.25 * (
fabs(m_xmax - m_xmin) +
fabs(m_ymax - m_ymin));
1629 m_zmax = +0.25 * (
fabs(m_xmax - m_xmin) +
fabs(m_ymax - m_ymin));
1634 if (!(setx && sety && setz)) {
1636 std::cerr <<
" Unable to establish"
1637 <<
" default dimensions in all directions.\n";
1641 if (vmin == vmax || !setv) {
1643 std::cerr <<
" All potentials in the cell are the same.\n";
1644 std::cerr <<
" There is no point in going on.\n";
1652bool ComponentAnalyticField::CellType() {
1658 m_scellType =
"D2 ";
1661 m_scellType =
"D1 ";
1664 }
else if (m_ntube >= 3 && m_ntube <= 8) {
1666 m_scellType =
"D4 ";
1669 m_scellType =
"D3 ";
1674 <<
" Potentials for tube with " << m_ntube
1675 <<
" edges are not yet available.\n"
1676 <<
" Using a round tube instead.\n";
1677 m_scellType =
"D3 ";
1685 if (!(m_perx || m_pery) && !(m_ynplan[0] && m_ynplan[1]) &&
1686 !(m_ynplan[2] && m_ynplan[3])) {
1693 if (m_perx && !m_pery && !(m_ynplan[0] || m_ynplan[1]) &&
1694 !(m_ynplan[2] && m_ynplan[3])) {
1695 m_scellType =
"B1X";
1701 if (m_pery && !m_perx && !(m_ynplan[0] && m_ynplan[1]) &&
1702 !(m_ynplan[2] || m_ynplan[3])) {
1703 m_scellType =
"B1Y";
1709 if (m_perx && !m_pery && !(m_ynplan[2] && m_ynplan[3])) {
1710 m_scellType =
"B2X";
1715 if (!(m_perx || m_pery) && !(m_ynplan[2] && m_ynplan[3]) &&
1716 (m_ynplan[0] && m_ynplan[1])) {
1717 m_sx =
fabs(m_coplan[1] - m_coplan[0]);
1718 m_scellType =
"B2X";
1724 if (m_pery && !m_perx && !(m_ynplan[0] && m_ynplan[1])) {
1725 m_scellType =
"B2Y";
1730 if (!(m_perx || m_pery) && !(m_ynplan[0] && m_ynplan[1]) &&
1731 (m_ynplan[2] && m_ynplan[3])) {
1732 m_sy =
fabs(m_coplan[3] - m_coplan[2]);
1733 m_scellType =
"B2Y";
1739 if (!(m_ynplan[0] || m_ynplan[1] || m_ynplan[2] || m_ynplan[3]) && m_perx && m_pery) {
1740 m_scellType =
"C1 ";
1746 if (!((m_ynplan[2] && m_pery) || (m_ynplan[2] && m_ynplan[3]))) {
1747 if (m_ynplan[0] && m_ynplan[1]) {
1748 m_sx =
fabs(m_coplan[1] - m_coplan[0]);
1749 m_scellType =
"C2X";
1753 if (m_perx && m_ynplan[0]) {
1754 m_scellType =
"C2X";
1761 if (!((m_ynplan[0] && m_perx) || (m_ynplan[0] && m_ynplan[1]))) {
1762 if (m_ynplan[2] && m_ynplan[3]) {
1763 m_sy =
fabs(m_coplan[3] - m_coplan[2]);
1764 m_scellType =
"C2Y";
1768 if (m_pery && m_ynplan[2]) {
1769 m_scellType =
"C2Y";
1776 if (m_perx && m_pery) {
1777 m_scellType =
"C3 ";
1783 m_sy =
fabs(m_coplan[3] - m_coplan[2]);
1784 m_scellType =
"C3 ";
1790 m_sx =
fabs(m_coplan[1] - m_coplan[0]);
1791 m_scellType =
"C3 ";
1796 if (m_ynplan[0] && m_ynplan[1] && m_ynplan[2] && m_ynplan[3]) {
1797 m_scellType =
"C3 ";
1798 m_sx =
fabs(m_coplan[1] - m_coplan[0]);
1799 m_sy =
fabs(m_coplan[3] - m_coplan[2]);
1808bool ComponentAnalyticField::PrepareStrips() {
1815 double gapDef[4] = {0., 0., 0., 0.};
1820 gapDef[0] = m_coplan[1] - m_coplan[0];
1821 }
else if (m_nWires <= 0) {
1824 gapDef[0] = m_w[0].x - m_coplan[0];
1825 for (
int i = m_nWires; i--;) {
1826 if (m_w[i].x - m_coplan[0] < gapDef[0]) gapDef[0] = m_w[i].x - m_coplan[0];
1833 gapDef[1] = m_coplan[1] - m_coplan[0];
1834 }
else if (m_nWires <= 0) {
1837 gapDef[1] = m_coplan[1] - m_w[0].x;
1838 for (
int i = m_nWires; i--;) {
1839 if (m_coplan[1] - m_w[i].x < gapDef[1]) gapDef[1] = m_coplan[1] - m_w[i].x;
1846 gapDef[2] = m_coplan[3] - m_coplan[2];
1847 }
else if (m_nWires <= 0) {
1850 gapDef[2] = m_w[0].y - m_coplan[2];
1851 for (
int i = m_nWires; i--;) {
1852 if (m_w[i].y - m_coplan[2] < gapDef[2]) gapDef[2] = m_w[i].y - m_coplan[2];
1859 gapDef[3] = m_coplan[3] - m_coplan[2];
1860 }
else if (m_nWires <= 0) {
1863 gapDef[3] = m_coplan[3] - m_w[0].y;
1864 for (
int i = m_nWires; i--;) {
1865 if (m_coplan[3] - m_w[i].y < gapDef[3]) gapDef[3] = m_coplan[3] - m_w[i].y;
1871 for (
unsigned int i = 0; i < 4; ++i) {
1872 const unsigned int nStrips1 = planes[i].strips1.size();
1873 for (
unsigned int j = 0; j < nStrips1; ++j) {
1874 if (planes[i].strips1[j].gap < 0.) {
1875 planes[i].strips1[j].gap = gapDef[i];
1877 if (planes[i].strips1[j].gap < 0.) {
1879 std::cerr <<
" Not able to set a default anode-cathode gap\n";
1880 std::cerr <<
" for x/y-strip " << j <<
" of plane " << i <<
".\n";
1884 const unsigned int nStrips2 = planes[i].strips2.size();
1885 for (
unsigned int j = 0; j < nStrips2; ++j) {
1886 if (planes[i].strips2[j].gap < 0.) {
1887 planes[i].strips2[j].gap = gapDef[i];
1889 if (planes[i].strips2[j].gap < 0.) {
1891 std::cerr <<
" Not able to set a default anode-cathode gap\n";
1892 std::cerr <<
" for z-strip " << j <<
" of plane " << i <<
".\n";
1896 const unsigned int nPixels = planes[i].pixels.size();
1897 for (
unsigned int j = 0; j < nPixels; ++j) {
1898 if (planes[i].pixels[j].gap < 0.) {
1899 planes[i].pixels[j].gap = gapDef[i];
1901 if (planes[i].pixels[j].gap < 0.) {
1903 std::cerr <<
" Not able to set a default anode-cathode gap\n";
1904 std::cerr <<
" for pixel " << j <<
" of plane " << i <<
".\n";
1916 if (std::find(m_readout.begin(), m_readout.end(), label) != m_readout.end()) {
1918 std::cout <<
" Readout group " << label <<
" already exists.\n";
1921 m_readout.push_back(label);
1923 unsigned int nWiresFound = 0;
1924 for (
unsigned int i = 0; i < m_nWires; ++i) {
1925 if (m_w[i].type == label) ++nWiresFound;
1928 unsigned int nPlanesFound = 0;
1929 unsigned int nStripsFound = 0;
1930 unsigned int nPixelsFound = 0;
1931 for (
int i = 0; i < 5; ++i) {
1932 if (planes[i].type == label) ++nPlanesFound;
1933 const unsigned int nStrips1 = planes[i].strips1.size();
1934 for (
unsigned int j = 0; j < nStrips1; ++j) {
1935 if (planes[i].strips1[j].type == label) ++nStripsFound;
1937 const unsigned int nStrips2 = planes[i].strips2.size();
1938 for (
unsigned int j = 0; j < nStrips2; ++j) {
1939 if (planes[i].strips2[j].type == label) ++nStripsFound;
1941 const unsigned int nPixels = planes[i].pixels.size();
1942 for (
unsigned int j = 0; j < nPixels; ++j) {
1943 if (planes[i].pixels[j].type == label) ++nPixelsFound;
1947 if (nWiresFound == 0 && nPlanesFound == 0 && nStripsFound == 0 &&
1948 nPixelsFound == 0) {
1950 std::cerr <<
" At present there are no wires, planes or strips\n";
1951 std::cerr <<
" associated to readout group " << label <<
".\n";
1954 std::cout <<
" Readout group " << label <<
" comprises:\n";
1955 if (nWiresFound > 1) {
1956 std::cout <<
" " << nWiresFound <<
" wires\n";
1957 }
else if (nWiresFound == 1) {
1958 std::cout <<
" 1 wire\n";
1960 if (nPlanesFound > 1) {
1961 std::cout <<
" " << nPlanesFound <<
" planes\n";
1962 }
else if (nPlanesFound == 1) {
1963 std::cout <<
" 1 plane\n";
1965 if (nStripsFound > 1) {
1966 std::cout <<
" " << nStripsFound <<
" strips\n";
1967 }
else if (nStripsFound == 1) {
1968 std::cout <<
" 1 strip\n";
1970 if (nPixelsFound > 1) {
1971 std::cout <<
" " << nPixelsFound <<
" pixels\n";
1972 }
else if (nPixelsFound == 1) {
1973 std::cout <<
" 1 pixel\n";
1980bool ComponentAnalyticField::Setup() {
1989 m_coplax = m_coplan[0];
1991 }
else if (m_ynplan[1]) {
1992 m_coplax = m_coplan[1];
1999 m_coplay = m_coplan[2];
2001 }
else if (m_ynplan[3]) {
2002 m_coplay = m_coplan[3];
2012 m_corvtc = m_vttube;
2013 }
else if ((m_ynplan[0] && m_ynplan[1]) && !(m_ynplan[2] || m_ynplan[3])) {
2014 m_corvta = (m_vtplan[0] - m_vtplan[1]) / (m_coplan[0] - m_coplan[1]);
2016 m_corvtc = (m_vtplan[1] * m_coplan[0] - m_vtplan[0] * m_coplan[1]) /
2017 (m_coplan[0] - m_coplan[1]);
2018 }
else if ((m_ynplan[2] && m_ynplan[3]) && !(m_ynplan[0] || m_ynplan[1])) {
2020 m_corvtb = (m_vtplan[2] - m_vtplan[3]) / (m_coplan[2] - m_coplan[3]);
2021 m_corvtc = (m_vtplan[3] * m_coplan[2] - m_vtplan[2] * m_coplan[3]) /
2022 (m_coplan[2] - m_coplan[3]);
2024 m_corvta = m_corvtb = m_corvtc = 0.;
2025 if (m_ynplan[0]) m_corvtc = m_vtplan[0];
2026 if (m_ynplan[1]) m_corvtc = m_vtplan[1];
2027 if (m_ynplan[2]) m_corvtc = m_vtplan[2];
2028 if (m_ynplan[3]) m_corvtc = m_vtplan[3];
2032 if (m_nWires <= 0)
return true;
2035 m_a.assign(m_nWires, std::vector<double>(m_nWires, 0.));
2040 if (m_scellType ==
"A ") ok = SetupA00();
2041 if (m_scellType ==
"B1X") ok = SetupB1X();
2042 if (m_scellType ==
"B1Y") ok = SetupB1Y();
2043 if (m_scellType ==
"B2X") ok = SetupB2X();
2044 if (m_scellType ==
"B2Y") ok = SetupB2Y();
2045 if (m_scellType ==
"C1 ") ok = SetupC10();
2046 if (m_scellType ==
"C2X") ok = SetupC2X();
2047 if (m_scellType ==
"C2Y") ok = SetupC2Y();
2048 if (m_scellType ==
"C3 ") ok = SetupC30();
2049 if (m_scellType ==
"D1 ") ok = SetupD10();
2050 if (m_scellType ==
"D2 ") ok = SetupD20();
2051 if (m_scellType ==
"D3 ") ok = SetupD30();
2058 std::cerr <<
" Computing the dipole moments failed.\n";
2066 std::cerr <<
" Preparing the cell for field calculations"
2067 <<
" did not succeed.\n";
2073bool ComponentAnalyticField::SetupA00() {
2083 for (
unsigned int i = 0; i < m_nWires; ++i) {
2084 m_a[i][i] = 0.25 * m_w[i].d * m_w[i].d;
2086 if (m_ynplax) m_a[i][i] /= 4. *
pow(m_w[i].x - m_coplax, 2);
2087 if (m_ynplay) m_a[i][i] /= 4. *
pow(m_w[i].y - m_coplay, 2);
2089 if (m_ynplax && m_ynplay)
2091 4.0 * (
pow(m_w[i].x - m_coplax, 2) +
pow(m_w[i].y - m_coplay, 2));
2093 m_a[i][i] = -0.5 * log(m_a[i][i]);
2095 for (
unsigned int j = i + 1; j < m_nWires; ++j) {
2096 m_a[i][j] =
pow(m_w[i].x - m_w[j].x, 2) +
pow(m_w[i].y - m_w[j].y, 2);
2099 m_a[i][j] = m_a[i][j] / (
pow(m_w[i].x + m_w[j].x - 2. * m_coplax, 2) +
2100 pow(m_w[i].y - m_w[j].y, 2));
2102 m_a[i][j] = m_a[i][j] / (
pow(m_w[i].x - m_w[j].x, 2) +
2103 pow(m_w[i].y + m_w[j].y - 2. * m_coplay, 2));
2105 if (m_ynplax && m_ynplay)
2106 m_a[i][j] *=
pow(m_w[i].x + m_w[j].x - 2. * m_coplax, 2) +
2107 pow(m_w[i].y + m_w[j].y - 2. * m_coplay, 2);
2109 m_a[i][j] = -0.5 * log(m_a[i][j]);
2111 m_a[j][i] = m_a[i][j];
2118bool ComponentAnalyticField::SetupB1X() {
2130 double xx = 0., yy = 0., yymirr = 0.;
2134 for (
unsigned int i = 0; i < m_nWires; ++i) {
2135 m_a[i][i] = -log(0.5 * m_w[i].d * Pi / m_sx);
2138 yy = (Pi / m_sx) * 2. * (m_w[i].y - m_coplay);
2139 if (
fabs(yy) > 20.) m_a[i][i] +=
fabs(yy) - CLog2;
2140 if (
fabs(yy) <= 20.) m_a[i][i] += log(
fabs(sinh(yy)));
2143 for (
unsigned int j = i + 1; j < m_nWires; ++j) {
2144 xx = (Pi / m_sx) * (m_w[i].x - m_w[j].x);
2145 yy = (Pi / m_sx) * (m_w[i].y - m_w[j].y);
2146 if (
fabs(yy) > 20.) m_a[i][j] = -
fabs(yy) + CLog2;
2147 if (
fabs(yy) <= 20.) {
2148 const double sinhy = sinh(yy);
2149 const double sinx =
sin(xx);
2150 m_a[i][j] = -0.5 * log(sinhy * sinhy + sinx * sinx);
2155 yymirr = (Pi / m_sx) * (m_w[i].y + m_w[j].y - 2. * m_coplay);
2156 if (
fabs(yymirr) > 20.) r2plan =
fabs(yymirr) - CLog2;
2157 if (
fabs(yymirr) <= 20.) {
2158 const double sinhy = sinh(yymirr);
2159 const double sinx =
sin(xx);
2160 r2plan = 0.5 * log(sinhy * sinhy + sinx * sinx);
2162 m_a[i][j] += r2plan;
2165 m_a[j][i] = m_a[i][j];
2172bool ComponentAnalyticField::SetupB1Y() {
2183 double xx = 0., yy = 0., xxmirr = 0.;
2187 for (
unsigned int i = 0; i < m_nWires; ++i) {
2188 m_a[i][i] = -log(0.5 * m_w[i].d * Pi / m_sy);
2191 xx = (Pi / m_sy) * 2. * (m_w[i].x - m_coplax);
2192 if (
fabs(xx) > 20.) m_a[i][i] +=
fabs(xx) - CLog2;
2193 if (
fabs(xx) <= 20.) m_a[i][i] += log(
fabs(sinh(xx)));
2196 for (
unsigned int j = i + 1; j < m_nWires; ++j) {
2197 xx = (Pi / m_sy) * (m_w[i].x - m_w[j].x);
2198 yy = (Pi / m_sy) * (m_w[i].y - m_w[j].y);
2199 if (
fabs(xx) > 20.) m_a[i][j] = -
fabs(xx) + CLog2;
2200 if (
fabs(xx) <= 20.) {
2201 const double sinhx = sinh(xx);
2202 const double siny =
sin(yy);
2203 m_a[i][j] = -0.5 * log(sinhx * sinhx + siny * siny);
2207 xxmirr = (Pi / m_sy) * (m_w[i].x + m_w[j].x - 2. * m_coplax);
2209 if (
fabs(xxmirr) > 20.) r2plan =
fabs(xxmirr) - CLog2;
2210 if (
fabs(xxmirr) <= 20.) {
2211 const double sinhx = sinh(xxmirr);
2212 const double siny =
sin(yy);
2213 r2plan = 0.5 * log(sinhx * sinhx + siny * siny);
2215 m_a[i][j] += r2plan;
2218 m_a[j][i] = m_a[i][j];
2225bool ComponentAnalyticField::SetupB2X() {
2238 m_b2sin.resize(m_nWires);
2241 for (
unsigned int i = 0; i < m_nWires; ++i) {
2242 double xx = (Pi / m_sx) * (m_w[i].x - m_coplax);
2243 m_a[i][i] = (0.25 * m_w[i].d * Pi / m_sx) /
sin(xx);
2246 const double yymirr = (Pi / m_sx) * (m_w[i].y - m_coplay);
2247 if (
fabs(yymirr) <= 20.) {
2248 const double sinhy = sinh(yymirr);
2249 const double sinx =
sin(xx);
2250 m_a[i][i] *=
sqrt(sinhy * sinhy + sinx * sinx) / sinhy;
2254 m_a[i][i] = -log(
fabs(m_a[i][i]));
2256 for (
unsigned int j = i + 1; j < m_nWires; ++j) {
2257 xx = HalfPi * (m_w[i].x - m_w[j].x) / m_sx;
2258 const double yy = HalfPi * (m_w[i].y - m_w[j].y) / m_sx;
2259 const double xxneg = HalfPi * (m_w[i].x + m_w[j].x - 2. * m_coplax) / m_sx;
2260 if (
fabs(yy) <= 20.) {
2261 const double sinhy = sinh(yy);
2262 const double sinxx =
sin(xx);
2263 const double sinxxneg =
sin(xxneg);
2264 m_a[i][j] = (sinhy * sinhy + sinxx * sinxx) /
2265 (sinhy * sinhy + sinxxneg * sinxxneg);
2267 if (
fabs(yy) > 20.) m_a[i][j] = 1.0;
2270 const double yymirr =
2271 HalfPi * (m_w[i].y + m_w[j].y - 2. * m_coplay) / m_sx;
2272 if (
fabs(yymirr) <= 20.) {
2273 const double sinhy = sinh(yymirr);
2274 const double sinxx =
sin(xx);
2275 const double sinxxneg =
sin(xxneg);
2276 m_a[i][j] *= (sinhy * sinhy + sinxxneg * sinxxneg) /
2277 (sinhy * sinhy + sinxx * sinxx);
2281 m_a[i][j] = -0.5 * log(m_a[i][j]);
2282 m_a[j][i] = m_a[i][j];
2285 m_b2sin[i] =
sin(Pi * (m_coplax - m_w[i].x) / m_sx);
2291bool ComponentAnalyticField::SetupB2Y() {
2304 m_b2sin.resize(m_nWires);
2307 for (
unsigned int i = 0; i < m_nWires; ++i) {
2308 double yy = (Pi / m_sy) * (m_w[i].y - m_coplay);
2309 m_a[i][i] = (0.25 * m_w[i].d * Pi / m_sy) /
sin(yy);
2312 const double xxmirr = (Pi / m_sy) * (m_w[i].x - m_coplax);
2313 if (
fabs(xxmirr) <= 20.) {
2314 const double sinhx = sinh(xxmirr);
2315 const double sinyy =
sin(yy);
2316 m_a[i][i] *=
sqrt(sinhx * sinhx + sinyy * sinyy) / sinhx;
2320 m_a[i][i] = -log(
fabs(m_a[i][i]));
2322 for (
unsigned int j = i + 1; j < m_nWires; j++) {
2323 const double xx = HalfPi * (m_w[i].x - m_w[j].x) / m_sy;
2324 yy = HalfPi * (m_w[i].y - m_w[j].y) / m_sy;
2325 const double yyneg = HalfPi * (m_w[i].y + m_w[j].y - 2. * m_coplay) / m_sy;
2326 if (
fabs(xx) <= 20.) {
2327 const double sinhx = sinh(xx);
2328 const double sinyy =
sin(yy);
2329 const double sinyyneg =
sin(yyneg);
2330 m_a[i][j] = (sinhx * sinhx + sinyy * sinyy) /
2331 (sinhx * sinhx + sinyyneg * sinyyneg);
2333 if (
fabs(xx) > 20.) m_a[i][j] = 1.0;
2336 const double xxmirr =
2337 HalfPi * (m_w[i].x + m_w[j].x - 2. * m_coplax) / m_sy;
2338 if (
fabs(xxmirr) <= 20.) {
2339 const double sinhx = sinh(xxmirr);
2340 const double sinyy =
sin(yy);
2341 const double sinyyneg =
sin(yyneg);
2342 m_a[i][j] *= (sinhx * sinhx + sinyyneg * sinyyneg) /
2343 (sinhx * sinhx + sinyy * sinyy);
2347 m_a[i][j] = -0.5 * log(m_a[i][j]);
2348 m_a[j][i] = m_a[i][j];
2351 m_b2sin[i] =
sin(Pi * (m_coplay - m_w[i].y) / m_sy);
2357bool ComponentAnalyticField::SetupC10() {
2375 if (m_sy / m_sx < 8.) p =
exp(-Pi * m_sy / m_sx);
2376 m_zmult = std::complex<double>(Pi / m_sx, 0.);
2379 if (m_sx / m_sy < 8.) p =
exp(-Pi * m_sx / m_sy);
2380 m_zmult = std::complex<double>(0., Pi / m_sy);
2383 if (m_p1 > 1.e-10) m_p2 =
pow(p, 6);
2387 std::cout <<
" p, p1, p2 = " << p <<
", " << m_p1 <<
", " << m_p2
2389 std::cout <<
" zmult = " << m_zmult <<
"\n";
2390 std::cout <<
" mode = " << m_mode <<
"\n";
2394 for (
unsigned int i = 0; i < m_nWires; ++i) {
2395 for (
unsigned int j = 0; j < m_nWires; ++j) {
2396 const double xyi = m_mode == 0 ? m_w[i].x : m_w[i].y;
2397 const double xyj = m_mode == 0 ? m_w[j].x : m_w[j].y;
2398 const double temp = xyi * xyj * TwoPi / (m_sx * m_sy);
2400 m_a[i][j] = Ph2Lim(0.5 * m_w[i].d) - temp;
2402 m_a[i][j] = Ph2(m_w[i].x - m_w[j].x, m_w[i].y - m_w[j].y) - temp;
2407 if (!Charge())
return false;
2410 for (
unsigned int j = 0; j < m_nWires; ++j) {
2411 const double xyj = m_mode == 0 ? m_w[j].x : m_w[j].y;
2412 s += m_w[j].e * xyj;
2414 m_c1 = -s * 2. * Pi / (m_sx * m_sy);
2418bool ComponentAnalyticField::SetupC2X() {
2433 if (2. * m_sx <= m_sy) {
2435 if (m_sy / m_sx < 25.) p =
exp(-HalfPi * m_sy / m_sx);
2436 m_zmult = std::complex<double>(HalfPi / m_sx, 0.);
2439 if (m_sx / m_sy < 6.) p =
exp(-2. * Pi * m_sx / m_sy);
2440 m_zmult = std::complex<double>(0., Pi / m_sy);
2443 if (m_p1 > 1.e-10) m_p2 =
pow(p, 6);
2447 std::cout <<
" p, p1, p2 = " << p <<
", " << m_p1 <<
", " << m_p2
2449 std::cout <<
" zmult = " << m_zmult <<
"\n";
2450 std::cout <<
" mode = " << m_mode <<
"\n";
2454 for (
unsigned int i = 0; i < m_nWires; ++i) {
2455 const double cx = m_coplax - m_sx * int(round((m_coplax - m_w[i].x) / m_sx));
2456 for (
unsigned int j = 0; j < m_nWires; ++j) {
2459 temp = (m_w[i].x - cx) * (m_w[j].x - cx) * TwoPi / (m_sx * m_sy);
2463 Ph2Lim(0.5 * m_w[i].d) - Ph2(2. * (m_w[i].x - cx), 0.) - temp;
2465 m_a[i][j] = Ph2(m_w[i].x - m_w[j].x, m_w[i].y - m_w[j].y) -
2466 Ph2(m_w[i].x + m_w[j].x - 2. * cx, m_w[i].y - m_w[j].y) -
2472 if (!Charge())
return false;
2477 for (
unsigned int i = 0; i < m_nWires; ++i) {
2478 const double cx = m_coplax - m_sx * int(round((m_coplax - m_w[i].x) / m_sx));
2479 s += m_w[i].e * (m_w[i].x - cx);
2481 m_c1 = -s * TwoPi / (m_sx * m_sy);
2486bool ComponentAnalyticField::SetupC2Y() {
2501 if (m_sx <= 2. * m_sy) {
2503 if (m_sy / m_sx <= 6.) p =
exp(-2. * Pi * m_sy / m_sx);
2504 m_zmult = std::complex<double>(Pi / m_sx, 0.);
2507 if (m_sx / m_sy <= 25.) p =
exp(-HalfPi * m_sx / m_sy);
2508 m_zmult = std::complex<double>(0., HalfPi / m_sy);
2511 if (m_p1 > 1.e-10) m_p2 =
pow(p, 6);
2515 std::cout <<
" p, p1, p2 = " << p <<
", " << m_p1 <<
", " << m_p2
2517 std::cout <<
" zmult = " << m_zmult <<
"\n";
2518 std::cout <<
" mode = " << m_mode <<
"\n";
2522 for (
unsigned int i = 0; i < m_nWires; ++i) {
2523 const double cy = m_coplay - m_sy * int(round((m_coplay - m_w[i].y) / m_sy));
2524 for (
unsigned int j = 0; j < m_nWires; ++j) {
2527 temp = (m_w[i].y - cy) * (m_w[j].y - cy) * TwoPi / (m_sx * m_sy);
2531 Ph2Lim(0.5 * m_w[i].d) - Ph2(0., 2. * (m_w[j].y - cy)) - temp;
2533 m_a[i][j] = Ph2(m_w[i].x - m_w[j].x, m_w[i].y - m_w[j].y) -
2534 Ph2(m_w[i].x - m_w[j].x, m_w[i].y + m_w[j].y - 2. * cy) -
2540 if (!Charge())
return false;
2545 for (
unsigned int i = 0; i < m_nWires; ++i) {
2546 const double cy = m_coplay - m_sy * int(round((m_coplay - m_w[i].y) / m_sy));
2547 s += m_w[i].e * (m_w[i].y - cy);
2549 m_c1 = -s * TwoPi / (m_sx * m_sy);
2554bool ComponentAnalyticField::SetupC30() {
2571 if (m_sy / m_sx <= 13.) p =
exp(-Pi * m_sy / m_sx);
2572 m_zmult = std::complex<double>(HalfPi / m_sx, 0.);
2575 if (m_sx / m_sy <= 13.) p =
exp(-Pi * m_sx / m_sy);
2576 m_zmult = std::complex<double>(0., HalfPi / m_sy);
2579 if (m_p1 > 1.e-10) m_p2 =
pow(p, 6);
2583 std::cout <<
" p, p1, p2 = " << p <<
", " << m_p1 <<
", " << m_p2
2585 std::cout <<
" zmult = " << m_zmult <<
"\n";
2586 std::cout <<
" mode = " << m_mode <<
"\n";
2590 for (
unsigned int i = 0; i < m_nWires; ++i) {
2591 double cx = m_coplax - m_sx * int(round((m_coplax - m_w[i].x) / m_sx));
2592 double cy = m_coplay - m_sy * int(round((m_coplay - m_w[i].y) / m_sy));
2593 for (
unsigned int j = 0; j < m_nWires; ++j) {
2595 m_a[i][i] = Ph2Lim(0.5 * m_w[i].d) - Ph2(0., 2. * (m_w[i].y - cy)) -
2596 Ph2(2. * (m_w[i].x - cx), 0.) +
2597 Ph2(2. * (m_w[i].x - cx), 2. * (m_w[i].y - cy));
2600 Ph2(m_w[i].x - m_w[j].x, m_w[i].y - m_w[j].y) -
2601 Ph2(m_w[i].x - m_w[j].x, m_w[i].y + m_w[j].y - 2. * cy) -
2602 Ph2(m_w[i].x + m_w[j].x - 2. * cx, m_w[i].y - m_w[j].y) +
2603 Ph2(m_w[i].x + m_w[j].x - 2. * cx, m_w[i].y + m_w[j].y - 2. * cy);
2608 if (!Charge())
return false;
2614bool ComponentAnalyticField::SetupD10() {
2624 for (
unsigned int i = 0; i < m_nWires; ++i) {
2626 m_a[i][i] = -log(0.5 * m_w[i].d * m_cotube /
2627 (m_cotube2 - (m_w[i].x * m_w[i].x + m_w[i].y * m_w[i].y)));
2629 std::complex<double> zi(m_w[i].x, m_w[i].y);
2631 for (
unsigned int j = i + 1; j < m_nWires; ++j) {
2633 std::complex<double> zj(m_w[j].x, m_w[j].y);
2634 m_a[i][j] = -log(abs(m_cotube * (zi - zj) / (m_cotube2 - conj(zi) * zj)));
2636 m_a[j][i] = m_a[i][j];
2643bool ComponentAnalyticField::SetupD20() {
2654 for (
unsigned int i = 0; i < m_nWires; ++i) {
2656 std::complex<double> zi(m_w[i].x, m_w[i].y);
2657 if (abs(zi) < 0.5 * m_w[i].d) {
2660 for (
unsigned int j = 0; j < m_nWires; ++j) {
2664 -log(0.5 * m_w[i].d /
2666 (m_w[i].x * m_w[i].x + m_w[i].y * m_w[i].y) / m_cotube));
2669 std::complex<double> zj(m_w[j].x, m_w[j].y);
2670 m_a[j][i] = -log(abs((1. / m_cotube) * (zi - zj) /
2671 (1. - conj(zi) * zj / m_cotube2)));
2677 for (
unsigned int j = 0; j < m_nWires; ++j) {
2681 -log(abs(0.5 * m_w[i].d * m_mtube *
pow(zi, m_mtube - 1) /
2682 (
pow(m_cotube, m_mtube) *
2683 (1. -
pow((abs(zi) / m_cotube), 2 * m_mtube)))));
2686 std::complex<double> zj(m_w[j].x, m_w[j].y);
2687 m_a[j][i] = -log(abs((1 /
pow(m_cotube, m_mtube)) *
2688 (
pow(zj, m_mtube) -
pow(zi, m_mtube)) /
2689 (1. -
pow(zj * conj(zi) / m_cotube2, m_mtube))));
2698bool ComponentAnalyticField::SetupD30() {
2708 wmap.assign(m_nWires, std::complex<double>(0., 0.));
2710 std::complex<double> wd = std::complex<double>(0., 0.);
2712 InitializeCoefficientTables();
2715 m_kappa = tgamma((m_ntube + 1.) / m_ntube) *
2716 tgamma((m_ntube - 2.) / m_ntube) / tgamma((m_ntube - 1.) / m_ntube);
2718 for (
unsigned int i = 0; i < m_nWires; ++i) {
2720 ConformalMap(std::complex<double>(m_w[i].x, m_w[i].y) / m_cotube, wmap[i],
2723 m_a[i][i] = -log(abs((0.5 * m_w[i].d / m_cotube) * wd /
2724 (1. -
pow(abs(wmap[i]), 2))));
2726 for (
unsigned int j = 0; j < i; ++j) {
2728 -log(abs((wmap[i] - wmap[j]) / (1. - conj(wmap[i]) * wmap[j])));
2730 m_a[j][i] = m_a[i][j];
2737bool ComponentAnalyticField::Charge() {
2747 std::vector<double> b(m_nWires, 0.);
2748 for (
unsigned int i = 0; i < m_nWires; ++i) {
2749 b[i] = m_w[i].v - (m_corvta * m_w[i].x + m_corvtb * m_w[i].y + m_corvtc);
2756 if (!(m_ynplan[0] || m_ynplan[1] || m_ynplan[2] || m_ynplan[3] || m_tube)) {
2759 m_a.resize(m_nWires + 1);
2760 m_a[m_nWires].clear();
2761 for (
unsigned int i = 0; i < m_nWires; ++i) {
2762 m_a[i].push_back(1.);
2763 m_a[m_nWires].push_back(1.);
2765 m_a[m_nWires].push_back(0.);
2770 std::cerr <<
" Matrix inversion failed.\n";
2774 if (m_a[m_nWires][m_nWires] != 0.) {
2775 const double t = 1. / m_a[m_nWires][m_nWires];
2776 for (
unsigned int i = 0; i < m_nWires; ++i) {
2777 for (
unsigned int j = 0; j < m_nWires; ++j) {
2778 m_a[i][j] -= t * m_a[i][m_nWires] * m_a[m_nWires][j];
2783 std::cerr <<
" True inverse of the capacitance matrix"
2784 <<
" could not be calculated.\n";
2785 std::cerr <<
" Use of the FACTOR instruction should be avoided.\n";
2800 std::cerr <<
" Failure to solve the capacitance equations.\n";
2801 std::cerr <<
" No charges are available.\n";
2806 for (
unsigned int i = 0; i < m_nWires; ++i) m_w[i].e = b[i];
2811 std::cout <<
" Dump of the capacitance matrix after inversion:\n";
2812 for (
unsigned int i = 0; i < m_nWires; i += 10) {
2813 for (
unsigned int j = 0; j < m_nWires; j += 10) {
2814 std::cout <<
" (Block " << i / 10 <<
", " << j / 10 <<
")\n";
2815 for (
unsigned int ii = 0; ii < 10; ++ii) {
2816 if (i + ii >= m_nWires)
break;
2817 for (
unsigned int jj = 0; jj < 10; ++jj) {
2818 if (j + jj >= m_nWires)
break;
2819 std::cout << std::setw(6) << m_a[i + ii][j + jj] <<
" ";
2827 std::cout <<
" End of the inverted capacitance matrix.\n";
2831 if (m_chargeCheck) {
2833 std::cout <<
" Quality check of the charge calculation.\n";
2834 std::cout <<
" Wire E as obtained E reconstructed\n";
2835 for (
unsigned int i = 0; i < m_nWires; ++i) {
2837 for (
unsigned int j = 0; j < m_nWires; ++j) {
2838 b[i] += m_a[i][j] * (m_w[j].v - m_v0 -
2839 (m_corvta * m_w[j].x + m_corvtb * m_w[j].y + m_corvtc));
2841 std::cout <<
" " << i <<
" " << m_w[i].e <<
" " << b[i]
2848double ComponentAnalyticField::Ph2(
const double xpos,
const double ypos)
const {
2864 std::complex<double> zeta = m_zmult * std::complex<double>(xpos, ypos);
2865 if (
fabs(imag(zeta)) < 10.) {
2866 std::complex<double> zsin =
sin(zeta);
2867 std::complex<double> zcof = 4. * zsin * zsin - 2.;
2868 std::complex<double> zu = -m_p1 - zcof * m_p2;
2869 std::complex<double> zunew = 1. - zcof * zu - m_p2;
2870 std::complex<double> zterm = (zunew + zu) * zsin;
2871 return -log(abs(zterm));
2874 return -
fabs(imag(zeta)) + CLog2;
2877void ComponentAnalyticField::ConformalMap(
const std::complex<double>& z,
2878 std::complex<double>& ww,
2879 std::complex<double>& wd)
const {
2891 const int nterm = 15;
2898 }
else if (abs(z) < 0.75) {
2901 std::complex<double> zterm =
pow(m_kappa * z, m_ntube);
2902 std::complex<double> wdsum = 0.;
2903 std::complex<double> wsum = m_cc1[m_ntube - 3][nterm];
2904 for (
int i = nterm; i--;) {
2905 wdsum = wsum + zterm * wdsum;
2906 wsum = m_cc1[m_ntube - 3][i] + zterm * wsum;
2909 ww = m_kappa * z * wsum;
2910 wd = m_kappa * (wsum + double(m_ntube) * zterm * wdsum);
2914 double arot = -TwoPi *
2915 int(round(atan2(imag(z), real(z)) * m_ntube / TwoPi)) /
2917 const std::complex<double> zz = z * std::complex<double>(
cos(arot),
sin(arot));
2919 std::complex<double> zterm =
2920 pow(m_kappa * (1. - zz), m_ntube / (m_ntube - 2.));
2921 std::complex<double> wdsum = 0.;
2922 std::complex<double> wsum = m_cc2[m_ntube - 3][nterm];
2923 for (
int i = nterm; i--;) {
2924 wdsum = wsum + zterm * wdsum;
2925 wsum = m_cc2[m_ntube - 3][i] + zterm * wsum;
2928 ww = std::complex<double>(
cos(arot), -
sin(arot)) * (1. - zterm * wsum);
2929 wd = m_ntube * m_kappa *
pow(m_kappa * (1. - zz), 2. / (m_ntube - 2.)) *
2930 (wsum + zterm * wdsum) / (m_ntube - 2.);
2934void ComponentAnalyticField::E2Sum(
const double xpos,
const double ypos,
2935 double& ex,
double& ey)
const {
2948 const std::complex<double> icons = std::complex<double>(0., 1.);
2950 std::complex<double> wsum = 0.;
2951 std::complex<double> zsin, zcof, zu, zunew;
2952 std::complex<double> zterm1, zterm2, zeta;
2954 for (
unsigned int j = 0; j < m_nWires; ++j) {
2955 zeta = m_zmult * std::complex<double>(xpos - m_w[j].x, ypos - m_w[j].y);
2956 if (imag(zeta) > 15.) {
2957 wsum -= m_w[j].e * icons;
2958 }
else if (imag(zeta) < -15.) {
2959 wsum += m_w[j].e * icons;
2962 zcof = 4. * zsin * zsin - 2.;
2963 zu = -m_p1 - zcof * m_p2;
2964 zunew = 1. - zcof * zu - m_p2;
2965 zterm1 = (zunew + zu) * zsin;
2966 zu = -3. * m_p1 - zcof * 5. * m_p2;
2967 zunew = 1. - zcof * zu - 5. * m_p2;
2968 zterm2 = (zunew - zu) *
cos(zeta);
2969 wsum += m_w[j].e * (zterm2 / zterm1);
2972 ex = -real(-m_zmult * wsum);
2973 ey = imag(-m_zmult * wsum);
2976void ComponentAnalyticField::FieldA00(
const double xpos,
const double ypos,
2977 double& ex,
double& ey,
double& volt,
2978 const bool opt)
const {
2998 double xxmirr = 0., yymirr = 0.;
3000 for (
int i = m_nWires; i--;) {
3001 const double xx = xpos - m_w[i].x;
3002 const double yy = ypos - m_w[i].y;
3003 double r2 = xx * xx + yy * yy;
3005 double exhelp = xx / r2;
3006 double eyhelp = yy / r2;
3009 xxmirr = m_w[i].x + (xpos - 2. * m_coplax);
3010 const double r2plan = xxmirr * xxmirr + yy * yy;
3011 exhelp -= xxmirr / r2plan;
3012 eyhelp -= yy / r2plan;
3017 yymirr = m_w[i].y + (ypos - 2. * m_coplay);
3018 const double r2plan = xx * xx + yymirr * yymirr;
3019 exhelp -= xx / r2plan;
3020 eyhelp -= yymirr / r2plan;
3024 if (m_ynplax && m_ynplay) {
3025 const double r2plan = xxmirr * xxmirr + yymirr * yymirr;
3026 exhelp += xxmirr / r2plan;
3027 eyhelp += yymirr / r2plan;
3031 if (opt) volt -= 0.5 * m_w[i].e * log(r2);
3032 ex += m_w[i].e * exhelp;
3033 ey += m_w[i].e * eyhelp;
3037void ComponentAnalyticField::FieldB1X(
const double xpos,
const double ypos,
3038 double& ex,
double& ey,
double& volt,
3039 const bool opt)
const {
3049 const std::complex<double> icons = std::complex<double>(0., 1.);
3051 std::complex<double> ecompl;
3059 const double tx = Pi / m_sx;
3061 for (
unsigned int i = 0; i < m_nWires; ++i) {
3062 const double xx = tx * (xpos - m_w[i].x);
3063 const double yy = tx * (ypos - m_w[i].y);
3064 const std::complex<double> zz(xx, yy);
3066 if (yy > 20.) ecompl = -icons;
3067 if (
fabs(yy) <= 20.)
3069 icons * (
exp(2. * icons * zz) + 1.) / (
exp(2. * icons * zz) - 1.);
3070 if (yy < -20.) ecompl = icons;
3072 if (
fabs(yy) > 20.) r2 = -
fabs(yy) + CLog2;
3073 if (
fabs(yy) <= 20.) r2 = -0.5 * log(
pow(sinh(yy), 2) +
pow(
sin(xx), 2));
3077 const double yymirr = tx * (ypos + m_w[i].y - 2. * m_coplay);
3078 const std::complex<double> zzmirr(xx, yymirr);
3079 if (yymirr > 20.) ecompl += icons;
3080 if (
fabs(yymirr) <= 20.)
3081 ecompl += -icons * (
exp(2. * icons * zzmirr) + 1.) /
3082 (
exp(2. * icons * zzmirr) - 1.);
3083 if (yymirr < -20.) ecompl += -icons;
3085 if (opt &&
fabs(yymirr) > 20.) r2 +=
fabs(yymirr) - CLog2;
3086 if (opt &&
fabs(yymirr) <= 20.)
3087 r2 += 0.5 * log(
pow(sinh(yymirr), 2) +
pow(
sin(xx), 2));
3090 ex += m_w[i].e * real(ecompl);
3091 ey -= m_w[i].e * imag(ecompl);
3092 if (opt) volt += m_w[i].e * r2;
3098void ComponentAnalyticField::FieldB1Y(
const double xpos,
const double ypos,
3099 double& ex,
double& ey,
double& volt,
3100 const bool opt)
const {
3110 std::complex<double> ecompl;
3118 const double ty = Pi / m_sy;
3120 for (
int i = m_nWires; i--;) {
3121 const double xx = ty * (xpos - m_w[i].x);
3122 const double yy = ty * (ypos - m_w[i].y);
3123 const std::complex<double> zz(xx, yy);
3125 if (xx > 20.) ecompl = 1.;
3126 if (
fabs(xx) <= 20.) ecompl = (
exp(2. * zz) + 1.) / (
exp(2. * zz) - 1.);
3127 if (xx < -20.) ecompl = -1.;
3129 if (
fabs(xx) > 20.) r2 = -
fabs(xx) + CLog2;
3130 if (
fabs(xx) <= 20.) r2 = -0.5 * log(
pow(sinh(xx), 2) +
pow(
sin(yy), 2));
3134 const double xxmirr = ty * (xpos + m_w[i].x - 2. * m_coplax);
3135 const std::complex<double> zzmirr(xxmirr, yy);
3136 if (xxmirr > 20.) ecompl -= 1.;
3137 if (xxmirr < -20.) ecompl += 1.;
3138 if (
fabs(xxmirr) <= 20.)
3139 ecompl -= (
exp(2. * zzmirr) + 1.) / (
exp(2. * zzmirr) - 1.);
3140 if (opt &&
fabs(xxmirr) > 20.) r2 +=
fabs(xxmirr) - CLog2;
3141 if (opt &&
fabs(xxmirr) <= 20.)
3142 r2 += 0.5 * log(
pow(sinh(xxmirr), 2) +
pow(
sin(yy), 2));
3145 ex += m_w[i].e * real(ecompl);
3146 ey -= m_w[i].e * imag(ecompl);
3147 if (opt) volt += m_w[i].e * r2;
3153void ComponentAnalyticField::FieldB2X(
const double xpos,
const double ypos,
3154 double& ex,
double& ey,
double& volt,
3155 const bool opt)
const {
3170 const double tx = HalfPi / m_sx;
3172 for (
unsigned int i = 0; i < m_nWires; ++i) {
3173 const double xx = tx * (xpos - m_w[i].x);
3174 const double yy = tx * (ypos - m_w[i].y);
3175 const double xxneg = tx * (xpos - m_w[i].x - 2 * m_coplax);
3176 const std::complex<double> zz(xx, yy);
3177 const std::complex<double> zzneg(xxneg, yy);
3179 std::complex<double> ecompl(0., 0.);
3181 if (
fabs(yy) <= 20.) {
3182 ecompl = -m_b2sin[i] / (
sin(zz) *
sin(zzneg));
3184 const double sinhy = sinh(yy);
3185 const double sinxx =
sin(xx);
3186 const double sinxxneg =
sin(xxneg);
3187 r2 = (sinhy * sinhy + sinxx * sinxx) /
3188 (sinhy * sinhy + sinxxneg * sinxxneg);
3193 const double yymirr = tx * (ypos + m_w[i].y - 2 * m_coplay);
3194 const std::complex<double> zzmirr(xx, yymirr);
3195 const std::complex<double> zznmirr(xxneg, yymirr);
3196 if (
fabs(yymirr) <= 20.) {
3197 ecompl += m_b2sin[i] / (
sin(zzmirr) *
sin(zznmirr));
3199 const double sinhy = sinh(yymirr);
3200 const double sinxx =
sin(xx);
3201 const double sinxxneg =
sin(xxneg);
3202 const double r2plan = (sinhy * sinhy + sinxx * sinxx) /
3203 (sinhy * sinhy + sinxxneg * sinxxneg);
3209 ex += m_w[i].e * real(ecompl);
3210 ey -= m_w[i].e * imag(ecompl);
3211 if (opt) volt -= 0.5 * m_w[i].e * log(r2);
3217void ComponentAnalyticField::FieldB2Y(
const double xpos,
const double ypos,
3218 double& ex,
double& ey,
double& volt,
3219 const bool opt)
const {
3230 const std::complex<double> icons(0., 1.);
3236 const double ty = HalfPi / m_sy;
3238 for (
unsigned int i = 0; i < m_nWires; ++i) {
3239 const double xx = ty * (xpos - m_w[i].x);
3240 const double yy = ty * (ypos - m_w[i].y);
3241 const double yyneg = ty * (ypos + m_w[i].y - 2 * m_coplay);
3243 std::complex<double> ecompl(0., 0.);
3245 if (
fabs(xx) <= 20.) {
3246 const std::complex<double> zz(xx, yy);
3247 const std::complex<double> zzneg(xx, yyneg);
3248 ecompl = icons * m_b2sin[i] / (
sin(icons * zz) *
sin(icons * zzneg));
3250 const double sinhx = sinh(xx);
3251 const double sinyy =
sin(yy);
3252 const double sinyyneg =
sin(yyneg);
3253 r2 = (sinhx * sinhx + sinyy * sinyy) /
3254 (sinhx * sinhx + sinyyneg * sinyyneg);
3259 const double xxmirr = ty * (xpos + m_w[i].x - 2 * m_coplax);
3260 if (
fabs(xxmirr) <= 20.) {
3261 const std::complex<double> zzmirr(xxmirr, yy);
3262 const std::complex<double> zznmirr(xxmirr, yyneg);
3264 icons * m_b2sin[i] / (
sin(icons * zzmirr) *
sin(icons * zznmirr));
3266 const double sinhx = sinh(xxmirr);
3267 const double sinyy =
sin(yy);
3268 const double sinyyneg =
sin(yyneg);
3269 const double r2plan = (sinhx * sinhx + sinyy * sinyy) /
3270 (sinhx * sinhx + sinyyneg * sinyyneg);
3276 ex += m_w[i].e * real(ecompl);
3277 ey -= m_w[i].e * imag(ecompl);
3278 if (opt) volt -= 0.5 * m_w[i].e * log(r2);
3284void ComponentAnalyticField::FieldC10(
const double xpos,
const double ypos,
3285 double& ex,
double& ey,
double& volt,
3286 const bool opt)
const {
3296 if (m_mode == 0) volt = m_v0 + m_c1 * xpos;
3297 if (m_mode == 1) volt = m_v0 + m_c1 * ypos;
3298 for (
unsigned int i = 0; i < m_nWires; ++i) {
3299 volt += m_w[i].e * Ph2(xpos - m_w[i].x, ypos - m_w[i].y);
3304 E2Sum(xpos, ypos, ex, ey);
3305 if (m_mode == 0) ex -= m_c1;
3306 if (m_mode == 1) ey -= m_c1;
3309void ComponentAnalyticField::FieldC2X(
const double xpos,
const double ypos,
3310 double& ex,
double& ey,
double& volt,
3311 const bool opt)
const {
3319 const std::complex<double> icons(0., 1.);
3321 std::complex<double> zsin, zcof, zu, zunew;
3322 std::complex<double> zterm1, zterm2, zeta;
3325 std::complex<double> wsum1 = 0.;
3326 std::complex<double> wsum2 = 0.;
3330 for (
int i = m_nWires; i--;) {
3332 zeta = m_zmult * std::complex<double>(xpos - m_w[i].x, ypos - m_w[i].y);
3333 if (imag(zeta) > 15.) {
3334 wsum1 -= m_w[i].e * icons;
3335 if (opt) volt -= m_w[i].e * (
fabs(imag(zeta)) - CLog2);
3336 }
else if (imag(zeta) < -15.) {
3337 wsum1 += m_w[i].e * icons;
3338 if (opt) volt -= m_w[i].e * (
fabs(imag(zeta)) - CLog2);
3341 zcof = 4. * zsin * zsin - 2.;
3342 zu = -m_p1 - zcof * m_p2;
3343 zunew = 1. - zcof * zu - m_p2;
3344 zterm1 = (zunew + zu) * zsin;
3345 zu = -3. * m_p1 - zcof * 5. * m_p2;
3346 zunew = 1. - zcof * zu - 5. * m_p2;
3347 zterm2 = (zunew - zu) *
cos(zeta);
3348 wsum1 += m_w[i].e * (zterm2 / zterm1);
3349 if (opt) volt -= m_w[i].e * log(abs(zterm1));
3352 double cx = m_coplax - m_sx * int(round((m_coplax - m_w[i].x) / m_sx));
3355 std::complex<double>(2. * cx - xpos - m_w[i].x, ypos - m_w[i].y);
3356 if (imag(zeta) > 15.) {
3357 wsum2 -= m_w[i].e * icons;
3358 if (opt) volt += m_w[i].e * (
fabs(imag(zeta)) - CLog2);
3359 }
else if (imag(zeta) < -15.) {
3360 wsum2 += m_w[i].e * icons;
3361 if (opt) volt += m_w[i].e * (
fabs(imag(zeta)) - CLog2);
3364 zcof = 4. * zsin * zsin - 2.;
3365 zu = -m_p1 - zcof * m_p2;
3366 zunew = 1. - zcof * zu - m_p2;
3367 zterm1 = (zunew + zu) * zsin;
3368 zu = -3. * m_p1 - zcof * 5. * m_p2;
3369 zunew = 1. - zcof * zu - 5. * m_p2;
3370 zterm2 = (zunew - zu) *
cos(zeta);
3371 wsum2 += m_w[i].e * (zterm2 / zterm1);
3372 if (opt) volt += m_w[i].e * log(abs(zterm1));
3375 if (opt && m_mode == 0) {
3376 volt -= TwoPi * m_w[i].e * (xpos - cx) * (m_w[i].x - cx) / (m_sx * m_sy);
3380 ex = real(m_zmult * (wsum1 + wsum2));
3381 ey = -imag(m_zmult * (wsum1 - wsum2));
3383 if (m_mode == 0) ex -= m_c1;
3386void ComponentAnalyticField::FieldC2Y(
const double xpos,
const double ypos,
3387 double& ex,
double& ey,
double& volt,
3388 const bool opt)
const {
3396 const std::complex<double> icons(0., 1.);
3398 std::complex<double> zsin, zcof, zu, zunew;
3399 std::complex<double> zterm1, zterm2, zeta;
3403 std::complex<double> wsum1 = 0.;
3404 std::complex<double> wsum2 = 0.;
3407 for (
int i = m_nWires; i--;) {
3409 zeta = m_zmult * std::complex<double>(xpos - m_w[i].x, ypos - m_w[i].y);
3410 if (imag(zeta) > 15.) {
3411 wsum1 -= m_w[i].e * icons;
3412 if (opt) volt -= m_w[i].e * (
fabs(imag(zeta)) - CLog2);
3413 }
else if (imag(zeta) < -15.) {
3414 wsum1 += m_w[i].e * icons;
3415 if (opt) volt -= m_w[i].e * (
fabs(imag(zeta)) - CLog2);
3418 zcof = 4. * zsin * zsin - 2.;
3419 zu = -m_p1 - zcof * m_p2;
3420 zunew = 1. - zcof * zu - m_p2;
3421 zterm1 = (zunew + zu) * zsin;
3422 zu = -3. * m_p1 - zcof * 5. * m_p2;
3423 zunew = 1. - zcof * zu - 5. * m_p2;
3424 zterm2 = (zunew - zu) *
cos(zeta);
3425 wsum1 += m_w[i].e * (zterm2 / zterm1);
3426 if (opt) volt -= m_w[i].e * log(abs(zterm1));
3429 const double cy = m_coplay - m_sy * int(round((m_coplay - m_w[i].y) / m_sy));
3432 std::complex<double>(xpos - m_w[i].x, 2 * cy - ypos - m_w[i].y);
3433 if (imag(zeta) > 15.) {
3434 wsum2 -= m_w[i].e * icons;
3435 if (opt) volt += m_w[i].e * (
fabs(imag(zeta)) - CLog2);
3436 }
else if (imag(zeta) < -15.) {
3437 wsum2 += m_w[i].e * icons;
3438 if (opt) volt += m_w[i].e * (
fabs(imag(zeta)) - CLog2);
3441 zcof = 4. * zsin * zsin - 2.;
3442 zu = -m_p1 - zcof * m_p2;
3443 zunew = 1. - zcof * zu - m_p2;
3444 zterm1 = (zunew + zu) * zsin;
3445 zu = -3. * m_p1 - zcof * 5. * m_p2;
3446 zunew = 1. - zcof * zu - 5. * m_p2;
3447 zterm2 = (zunew - zu) *
cos(zeta);
3448 wsum2 += m_w[i].e * (zterm2 / zterm1);
3449 if (opt) volt += m_w[i].e * log(abs(zterm1));
3452 if (opt && m_mode == 1) {
3453 volt -= TwoPi * m_w[i].e * (ypos - cy) * (m_w[i].y - cy) / (m_sx * m_sy);
3457 ex = real(m_zmult * (wsum1 - wsum2));
3458 ey = -imag(m_zmult * (wsum1 + wsum2));
3460 if (m_mode == 1) ey -= m_c1;
3463void ComponentAnalyticField::FieldC30(
const double xpos,
const double ypos,
3464 double& ex,
double& ey,
double& volt,
3465 const bool opt)
const {
3473 const std::complex<double> icons(0., 1.);
3475 std::complex<double> zsin, zcof, zu, zunew;
3476 std::complex<double> zterm1, zterm2, zeta;
3479 std::complex<double> wsum1 = 0.;
3480 std::complex<double> wsum2 = 0.;
3481 std::complex<double> wsum3 = 0.;
3482 std::complex<double> wsum4 = 0.;
3486 for (
unsigned int i = 0; i < m_nWires; ++i) {
3488 zeta = m_zmult * std::complex<double>(xpos - m_w[i].x, ypos - m_w[i].y);
3489 if (imag(zeta) > 15.) {
3490 wsum1 -= m_w[i].e * icons;
3491 if (opt) volt -= m_w[i].e * (
fabs(imag(zeta)) - CLog2);
3492 }
else if (imag(zeta) < -15.) {
3493 wsum1 += m_w[i].e * icons;
3494 if (opt) volt -= m_w[i].e * (
fabs(imag(zeta)) - CLog2);
3497 zcof = 4. * zsin * zsin - 2.;
3498 zu = -m_p1 - zcof * m_p2;
3499 zunew = 1. - zcof * zu - m_p2;
3500 zterm1 = (zunew + zu) * zsin;
3501 zu = -3. * m_p1 - zcof * 5. * m_p2;
3502 zunew = 1. - zcof * zu - 5. * m_p2;
3503 zterm2 = (zunew - zu) *
cos(zeta);
3504 wsum1 += m_w[i].e * (zterm2 / zterm1);
3505 if (opt) volt -= m_w[i].e * log(abs(zterm1));
3508 const double cx = m_coplax - m_sx * int(round((m_coplax - m_w[i].x) / m_sx));
3511 std::complex<double>(2. * cx - xpos - m_w[i].x, ypos - m_w[i].y);
3512 if (imag(zeta) > 15.) {
3513 wsum2 -= m_w[i].e * icons;
3514 if (opt) volt += m_w[i].e * (
fabs(imag(zeta)) - CLog2);
3515 }
else if (imag(zeta) < -15.) {
3516 wsum2 += m_w[i].e * icons;
3517 if (opt) volt += m_w[i].e * (
fabs(imag(zeta)) - CLog2);
3520 zcof = 4. * zsin * zsin - 2.;
3521 zu = -m_p1 - zcof * m_p2;
3522 zunew = 1. - zcof * zu - m_p2;
3523 zterm1 = (zunew + zu) * zsin;
3524 zu = -3. * m_p1 - zcof * 5. * m_p2;
3525 zunew = 1. - zcof * zu - 5. * m_p2;
3526 zterm2 = (zunew - zu) *
cos(zeta);
3527 wsum2 += m_w[i].e * (zterm2 / zterm1);
3528 if (opt) volt += m_w[i].e * log(abs(zterm1));
3531 const double cy = m_coplay - m_sy * int(round((m_coplay - m_w[i].y) / m_sy));
3534 std::complex<double>(xpos - m_w[i].x, 2. * cy - ypos - m_w[i].y);
3535 if (imag(zeta) > 15.) {
3536 wsum3 -= m_w[i].e * icons;
3537 if (opt) volt += m_w[i].e * (
fabs(imag(zeta)) - CLog2);
3538 }
else if (imag(zeta) < -15.) {
3539 wsum3 += m_w[i].e * icons;
3540 if (opt) volt += m_w[i].e * (
fabs(imag(zeta)) - CLog2);
3543 zcof = 4. * zsin * zsin - 2.;
3544 zu = -m_p1 - zcof * m_p2;
3545 zunew = 1. - zcof * zu - m_p2;
3546 zterm1 = (zunew + zu) * zsin;
3547 zu = -3. * m_p1 - zcof * 5. * m_p2;
3548 zunew = 1. - zcof * zu - 5. * m_p2;
3549 zterm2 = (zunew - zu) *
cos(zeta);
3550 wsum3 += m_w[i].e * (zterm2 / zterm1);
3551 if (opt) volt += m_w[i].e * log(abs(zterm1));
3554 zeta = m_zmult * std::complex<double>(2. * cx - xpos - m_w[i].x,
3555 2. * cy - ypos - m_w[i].y);
3556 if (imag(zeta) > 15.) {
3557 wsum4 -= m_w[i].e * icons;
3558 if (opt) volt -= m_w[i].e * (
fabs(imag(zeta)) - CLog2);
3559 }
else if (imag(zeta) < -15.) {
3560 wsum4 += m_w[i].e * icons;
3561 if (opt) volt -= m_w[i].e * (
fabs(imag(zeta)) - CLog2);
3564 zcof = 4. * zsin * zsin - 2.;
3565 zu = -m_p1 - zcof * m_p2;
3566 zunew = 1. - zcof * zu - m_p2;
3567 zterm1 = (zunew + zu) * zsin;
3568 zu = -3. * m_p1 - zcof * 5. * m_p2;
3569 zunew = 1. - zcof * zu - 5. * m_p2;
3570 zterm2 = (zunew - zu) *
cos(zeta);
3571 wsum4 += m_w[i].e * (zterm2 / zterm1);
3572 if (opt) volt -= m_w[i].e * log(abs(zterm1));
3576 ex = real(m_zmult * (wsum1 + wsum2 - wsum3 - wsum4));
3577 ey = -imag(m_zmult * (wsum1 - wsum2 + wsum3 - wsum4));
3580void ComponentAnalyticField::FieldD10(
const double xpos,
const double ypos,
3581 double& ex,
double& ey,
double& volt,
3582 const bool opt)
const {
3599 const std::complex<double> zpos = std::complex<double>(xpos, ypos);
3601 for (
int i = m_nWires; i--;) {
3602 const double q = m_w[i].e;
3604 const std::complex<double> zi(m_w[i].x, m_w[i].y);
3607 volt -= q * log(abs(m_cotube * (zpos - zi) / (m_cotube2 - zpos * conj(zi))));
3610 const std::complex<double> wi =
3611 1. / conj(zpos - zi) + zi / (m_cotube2 - conj(zpos) * zi);
3617void ComponentAnalyticField::FieldD20(
const double xpos,
const double ypos,
3618 double& ex,
double& ey,
double& volt,
3619 const bool opt)
const {
3636 const std::complex<double> zpos = std::complex<double>(xpos, ypos);
3638 for (
int i = m_nWires; i--;) {
3640 const std::complex<double> zi(m_w[i].x, m_w[i].y);
3642 if (abs(zi) > 0.5 * m_w[i].d) {
3645 volt -= m_w[i].e * log(abs((1. /
pow(m_cotube, m_mtube)) *
3646 (
pow(zpos, m_mtube) -
pow(zi, m_mtube)) /
3647 (1. -
pow(zpos * conj(zi) / m_cotube2, m_mtube))));
3650 const std::complex<double> wi =
3651 double(m_mtube) *
pow(conj(zpos), m_mtube - 1) *
3652 (1. / conj(
pow(zpos, m_mtube) -
pow(zi, m_mtube)) +
3654 (
pow(m_cotube, 2 * m_mtube) -
pow(conj(zpos) * zi, m_mtube)));
3655 ex += m_w[i].e * real(wi);
3656 ey += m_w[i].e * imag(wi);
3660 volt -= m_w[i].e * log(abs((1. / m_cotube) * (zpos - zi) /
3661 (1. - zpos * conj(zi) / m_cotube2)));
3663 const std::complex<double> wi =
3664 1. / conj(zpos - zi) + zi / (m_cotube2 - conj(zpos) * zi);
3666 ex += m_w[i].e * real(wi);
3667 ey += m_w[i].e * imag(wi);
3672void ComponentAnalyticField::FieldD30(
const double xpos,
const double ypos,
3673 double& ex,
double& ey,
double& volt,
3674 const bool opt)
const {
3690 std::complex<double> whelp;
3693 std::complex<double> wpos, wdpos;
3694 ConformalMap(std::complex<double>(xpos, ypos) / m_cotube, wpos, wdpos);
3696 for (
int i = m_nWires; i--;) {
3700 m_w[i].e * log(abs((wpos - wmap[i]) / (1. - wpos * conj(wmap[i]))));
3702 whelp = wdpos * (1. -
pow(abs(wmap[i]), 2)) /
3703 ((wpos - wmap[i]) * (1. - conj(wmap[i]) * wpos));
3705 ex += m_w[i].e * real(whelp);
3706 ey -= m_w[i].e * imag(whelp);
3712bool ComponentAnalyticField::InTube(
const double x0,
const double y0,
3713 const double a,
const int n)
const {
3723 if (x0 == 0. && y0 == 0.)
return true;
3727 if (x0 * x0 + y0 * y0 > a * a)
return false;
3731 if (n < 0 || n == 1 || n == 2) {
3733 std::cerr <<
" Invalid number of edges (n = " << n <<
")\n";
3739 double phi = atan2(y0, x0);
3740 if (phi < 0.)
phi += TwoPi;
3741 phi -= TwoPi * int(0.5 * n * phi / Pi) / n;
3743 if ((x0 * x0 + y0 * y0) *
pow(
cos(Pi / n - phi), 2) >
3744 a * a *
pow(
cos(Pi / n), 2))
3750void ComponentAnalyticField::InitializeCoefficientTables() {
3752 const int nterms = 16;
3755 m_cc1.assign(6, std::vector<double>(nterms, 0.));
3756 m_cc2.assign(6, std::vector<double>(nterms, 0.));
3759 const double cc13[nterms] = {
3760 0.1000000000e+01, -.1666666865e+00, 0.3174602985e-01, -.5731921643e-02,
3761 0.1040112227e-02, -.1886279933e-03, 0.3421107249e-04, -.6204730198e-05,
3762 0.1125329618e-05, -.2040969207e-06, 0.3701631357e-07, -.6713513301e-08,
3763 0.1217605794e-08, -.2208327132e-09, 0.4005162868e-10, -.7264017512e-11};
3764 const double cc23[nterms] = {
3765 0.3333333135e+00, -.5555555597e-01, 0.1014109328e-01, -.1837154618e-02,
3766 0.3332451452e-03, -.6043842586e-04, 0.1096152027e-04, -.1988050826e-05,
3767 0.3605655365e-06, -.6539443120e-07, 0.1186035448e-07, -.2151069323e-08,
3768 0.3901317047e-09, -.7075676156e-10, 0.1283289534e-10, -.2327455936e-11};
3771 const double cc14[nterms] = {
3772 0.1000000000e+01, -.1000000238e+00, 0.8333332837e-02, -.7051283028e-03,
3773 0.5967194738e-04, -.5049648280e-05, 0.4273189802e-06, -.3616123934e-07,
3774 0.3060091514e-08, -.2589557457e-09, 0.2191374859e-10, -.1854418528e-11,
3775 0.1569274224e-12, -.1327975205e-13, 0.1123779363e-14, -.9509817570e-16};
3776 const double cc24[nterms] = {
3777 0.1000000000e+01, -.5000000000e+00, 0.3000000119e+00, -.1750000119e+00,
3778 0.1016666889e+00, -.5916666612e-01, 0.3442307562e-01, -.2002724260e-01,
3779 0.1165192947e-01, -.6779119372e-02, 0.3944106400e-02, -.2294691978e-02,
3780 0.1335057430e-02, -.7767395582e-03, 0.4519091453e-03, -.2629216760e-03};
3783 const double cc15[nterms] = {
3784 0.1000000000e+01, -.6666666269e-01, 0.1212121220e-02, -.2626262140e-03,
3785 -.3322110570e-04, -.9413293810e-05, -.2570029210e-05, -.7695705904e-06,
3786 -.2422486887e-06, -.7945993730e-07, -.2691839640e-07, -.9361642128e-08,
3787 -.3327319087e-08, -.1204430555e-08, -.4428404310e-09, -.1650302672e-09};
3788 const double cc25[nterms] = {
3789 0.1248050690e+01, -.7788147926e+00, 0.6355384588e+00, -.4899077415e+00,
3790 0.3713272810e+00, -.2838423252e+00, 0.2174729109e+00, -.1663445234e+00,
3791 0.1271933913e+00, -.9728997946e-01, 0.7442557812e-01, -.5692918226e-01,
3792 0.4354400188e-01, -.3330700099e-01, 0.2547712997e-01, -.1948769018e-01};
3795 const double cc16[nterms] = {
3796 0.1000000000e+01, -.4761904851e-01, -.1221001148e-02, -.3753788769e-03,
3797 -.9415557724e-04, -.2862767724e-04, -.9587882232e-05, -.3441659828e-05,
3798 -.1299798896e-05, -.5103651119e-06, -.2066504408e-06, -.8578405186e-07,
3799 -.3635090096e-07, -.1567239494e-07, -.6857355572e-08, -.3038770346e-08};
3800 const double cc26[nterms] = {
3801 0.1333333015e+01, -.8888888955e+00, 0.8395061493e+00, -.7242798209e+00,
3802 0.6016069055e+00, -.5107235312e+00, 0.4393203855e+00, -.3745460510e+00,
3803 0.3175755739e+00, -.2703750730e+00, 0.2308617830e+00, -.1966916919e+00,
3804 0.1672732830e+00, -.1424439549e+00, 0.1214511395e+00, -.1034612656e+00};
3807 const double cc17[nterms] = {
3808 0.1000000000e+01, -.3571428731e-01, -.2040816238e-02, -.4936389159e-03,
3809 -.1446709794e-03, -.4963850370e-04, -.1877940667e-04, -.7600909157e-05,
3810 -.3232265954e-05, -.1427365532e-05, -.6493634714e-06, -.3026190711e-06,
3811 -.1438593245e-06, -.6953911225e-07, -.3409525462e-07, -.1692310647e-07};
3812 const double cc27[nterms] = {
3813 0.1359752655e+01, -.9244638681e+00, 0.9593217969e+00, -.8771237731e+00,
3814 0.7490229011e+00, -.6677658558e+00, 0.6196745634e+00, -.5591596961e+00,
3815 0.4905325770e+00, -.4393517375e+00, 0.4029803872e+00, -.3631100059e+00,
3816 0.3199430704e+00, -.2866140604e+00, 0.2627358437e+00, -.2368256450e+00};
3819 const double cc18[nterms] = {
3820 0.1000000000e+01, -.2777777612e-01, -.2246732125e-02, -.5571441725e-03,
3821 -.1790652314e-03, -.6708275760e-04, -.2766949183e-04, -.1219387286e-04,
3822 -.5640039490e-05, -.2706697160e-05, -.1337270078e-05, -.6763995657e-06,
3823 -.3488264610e-06, -.1828456675e-06, -.9718036154e-07, -.5227070332e-07};
3824 const double cc28[nterms] = {
3825 0.1362840652e+01, -.9286670089e+00, 0.1035511017e+01, -.9800255299e+00,
3826 0.8315343261e+00, -.7592730522e+00, 0.7612683773e+00, -.7132136226e+00,
3827 0.6074471474e+00, -.5554352999e+00, 0.5699443221e+00, -.5357525349e+00,
3828 0.4329345822e+00, -.3916820884e+00, 0.4401986003e+00, -.4197303057e+00};
3830 for (
int i = 0; i < nterms; ++i) {
3831 m_cc1[0][i] = cc13[i];
3832 m_cc2[0][i] = cc23[i];
3833 m_cc1[1][i] = cc14[i];
3834 m_cc2[1][i] = cc24[i];
3835 m_cc1[2][i] = cc15[i];
3836 m_cc2[2][i] = cc25[i];
3837 m_cc1[3][i] = cc16[i];
3838 m_cc2[3][i] = cc26[i];
3839 m_cc1[4][i] = cc17[i];
3840 m_cc2[4][i] = cc27[i];
3841 m_cc1[5][i] = cc18[i];
3842 m_cc2[5][i] = cc28[i];
3846void ComponentAnalyticField::Field3dA00(
const double xpos,
const double ypos,
3847 const double zpos,
double& ex,
3848 double& ey,
double& ez,
double& volt) {
3861 double exhelp, eyhelp, ezhelp, vhelp;
3864 ex = ey = ez = volt = 0.;
3867 const unsigned int n3d = m_ch3d.size();
3868 for (
unsigned int i = 0; i < n3d; ++i) {
3870 const double dx = xpos - m_ch3d[i].x;
3871 const double dy = ypos - m_ch3d[i].y;
3872 const double dz = zpos - m_ch3d[i].z;
3873 const double r =
sqrt(dx * dx + dy * dy + dz * dz);
3874 if (
fabs(r) < Small)
continue;
3875 const double r3 =
pow(r, 3);
3881 double dxm = 0., dym = 0.;
3883 dxm = m_ch3d[i].x + xpos - 2 * m_coplax;
3884 const double rplan =
sqrt(dxm * dxm + dy * dy);
3885 if (
fabs(rplan) < Small)
continue;
3886 const double rplan3 =
pow(rplan, 3);
3887 exhelp += dxm / rplan3;
3888 eyhelp += dy / rplan3;
3889 ezhelp += dz / rplan3;
3890 vhelp -= 1. / rplan;
3894 dym = m_ch3d[i].y + ypos - 2. * m_coplay;
3895 const double rplan =
sqrt(dx * dx + dym * dym);
3896 if (
fabs(rplan) < Small)
continue;
3897 const double rplan3 =
pow(rplan, 3);
3898 exhelp += dx / rplan3;
3899 eyhelp += dym / rplan3;
3900 ezhelp += dz / rplan3;
3901 vhelp -= 1. / rplan;
3904 if (m_ynplax && m_ynplay) {
3905 const double rplan =
sqrt(dxm * dxm + dym * dym);
3906 if (
fabs(rplan) < Small)
continue;
3907 const double rplan3 =
pow(rplan, 3);
3908 exhelp -= dxm / rplan3;
3909 eyhelp -= dym / rplan3;
3910 ezhelp -= dz / rplan3;
3911 vhelp += 1. / rplan;
3914 ex -= m_ch3d[i].e * exhelp;
3915 ey -= m_ch3d[i].e * eyhelp;
3916 ez -= m_ch3d[i].e * ezhelp;
3917 volt += m_ch3d[i].e * vhelp;
3921void ComponentAnalyticField::Field3dB2X(
const double xpos,
const double ypos,
3922 const double zpos,
double& ex,
3923 double& ey,
double& ez,
double& volt) {
3934 const double rcut = 1.;
3936 double rr1, rr2, rm1, rm2, err, ezz;
3937 double exsum = 0., eysum = 0., ezsum = 0., vsum = 0.;
3938 double k0r, k1r, k0rm, k1rm;
3941 ex = ey = ez = volt = 0.;
3944 const unsigned int n3d = m_ch3d.size();
3945 for (
unsigned int i = 0; i < n3d; ++i) {
3947 if (xpos == m_ch3d[i].x && ypos == m_ch3d[i].y && zpos == m_ch3d[i].z)
3949 const double dx = xpos - m_ch3d[i].x;
3950 const double dy = ypos - m_ch3d[i].y;
3951 const double dz = zpos - m_ch3d[i].z;
3952 const double dxm = xpos + m_ch3d[i].x - 2 * m_coplax;
3954 if (dy * dy + dz * dz >
pow(rcut * 2 * m_sx, 2)) {
3956 exsum = eysum = ezsum = vsum = 0.;
3958 for (
unsigned int j = 1; j <= m_nTermBessel; ++j) {
3960 const double rr = Pi * j *
sqrt(dy * dy + dz * dz) / m_sx;
3961 const double zzp = Pi * j * dx / m_sx;
3962 const double zzn = Pi * j * dxm / m_sx;
3972 const double czzp =
cos(zzp);
3973 const double czzn =
cos(zzn);
3974 vsum += (1. / m_sx) * k0r * (czzp - czzn);
3975 err = (TwoPi * j / (m_sx * m_sx)) * k1r * (czzp - czzn);
3976 ezz = (TwoPi * j / (m_sx * m_sx)) * k0r * (
sin(zzp) -
sin(zzn));
3978 eysum += err * dy /
sqrt(dy * dy + dz * dz);
3979 ezsum += err * dz /
sqrt(dy * dy + dz * dz);
3985 for (
unsigned int j = 0; j <= m_nTermPoly; ++j) {
3987 rr1 =
sqrt(
pow(dx + j * 2 * m_sx, 2) + dy * dy + dz * dz);
3988 rr2 =
sqrt(
pow(dx - j * 2 * m_sx, 2) + dy * dy + dz * dz);
3989 rm1 =
sqrt(
pow(dxm - j * 2 * m_sx, 2) + dy * dy + dz * dz);
3990 rm2 =
sqrt(
pow(dxm + j * 2 * m_sx, 2) + dy * dy + dz * dz);
3991 const double rr13 =
pow(rr1, 3);
3992 const double rm13 =
pow(rm1, 3);
3995 vsum = 1. / rr1 - 1. / rm1;
3996 exsum = dx / rr13 - dxm / rm13;
3997 eysum = dy * (1. / rr13 - 1. / rm13);
3998 ezsum = dz * (1. / rr13 - 1. / rm13);
4001 const double rr23 =
pow(rr2, 3);
4002 const double rm23 =
pow(rm2, 3);
4004 vsum += 1. / rr1 + 1. / rr2 - 1. / rm1 - 1. / rm2;
4005 exsum += (dx + j * 2 * m_sx) / rr13 + (dx - j * 2 * m_sx) / rr23 -
4006 (dxm - j * 2 * m_sx) / rm13 - (dxm + j * 2 * m_sx) / rm23;
4007 eysum += dy * (1. / rr13 + 1. / rr23 - 1. / rm13 - 1. / rm23);
4008 ezsum += dz * (1. / rr13 + 1. / rr23 - 1. / rm13 - 1. / rm23);
4013 const double dym = ypos + m_ch3d[i].y - 2. * m_coplay;
4014 if (dym * dym + dz * dz >
pow(rcut * 2 * m_sx, 2)) {
4017 for (
unsigned int j = 1; j <= m_nTermBessel; ++j) {
4019 const double rrm = Pi * j *
sqrt(dym * dym + dz * dz) / m_sx;
4020 const double zzp = Pi * j * dx / m_sx;
4021 const double zzn = Pi * j * dxm / m_sx;
4031 const double czzp =
cos(zzp);
4032 const double czzn =
cos(zzn);
4033 vsum += (1. / m_sx) * k0rm * (czzp - czzn);
4034 err = (TwoPi / (m_sx * m_sx)) * k1rm * (czzp - czzn);
4035 ezz = (TwoPi / (m_sx * m_sx)) * k0rm * (
sin(zzp) -
sin(zzn));
4037 eysum += err * dym /
sqrt(dym * dym + dz * dz);
4038 ezsum += err * dz /
sqrt(dym * dym + dz * dz);
4043 for (
unsigned int j = 0; j <= m_nTermPoly; ++j) {
4045 rr1 =
sqrt(
pow(dx + j * 2 * m_sx, 2) + dym * dym + dz * dz);
4046 rr2 =
sqrt(
pow(dx - j * 2 * m_sx, 2) + dym * dym + dz * dz);
4047 rm1 =
sqrt(
pow(dxm - j * 2 * m_sx, 2) + dym * dym + dz * dz);
4048 rm2 =
sqrt(
pow(dxm + j * 2 * m_sx, 2) + dym * dym + dz * dz);
4049 const double rr13 =
pow(rr1, 3);
4050 const double rm13 =
pow(rm1, 3);
4053 vsum += -1. / rr1 + 1. / rm1;
4054 exsum += -dx / rr13 + dxm / rm13;
4055 eysum += -dym * (1. / rr13 - 1. / rm13);
4056 ezsum += -dz * (1. / rr13 - 1. / rm13);
4059 const double rr23 =
pow(rr2, 3);
4060 const double rm23 =
pow(rm2, 3);
4062 vsum += -1. / rr1 - 1. / rr2 + 1. / rm1 + 1. / rm2;
4063 exsum += -(dx + j * 2 * m_sx) / rr13 - (dx - j * 2 * m_sx) / rr23 +
4064 (dxm - j * 2 * m_sx) / rm13 + (dxm + j * 2 * m_sx) / rm23;
4065 eysum += -dym * (1. / rr13 + 1. / rr23 - 1. / rm13 - 1. / rm23);
4066 ezsum += -dz * (1. / rr13 + 1. / rr23 - 1. / rm13 - 1. / rm23);
4070 ex += m_ch3d[i].e * exsum;
4071 ey += m_ch3d[i].e * eysum;
4072 ez += m_ch3d[i].e * ezsum;
4073 volt += m_ch3d[i].e * vsum;
4077void ComponentAnalyticField::Field3dB2Y(
const double xpos,
const double ypos,
4078 const double zpos,
double& ex,
4079 double& ey,
double& ez,
double& volt) {
4090 const double rcut = 1.;
4092 double rr1, rr2, rm1, rm2, err, ezz;
4093 double exsum = 0., eysum = 0., ezsum = 0., vsum = 0.;
4094 double k0r, k1r, k0rm, k1rm;
4097 ex = ey = ez = volt = 0.;
4100 const unsigned int n3d = m_ch3d.size();
4101 for (
unsigned int i = 0; i < n3d; ++i) {
4103 if (xpos == m_ch3d[i].x && ypos == m_ch3d[i].y && zpos == m_ch3d[i].z)
4105 const double dx = xpos - m_ch3d[i].x;
4106 const double dy = ypos - m_ch3d[i].y;
4107 const double dz = zpos - m_ch3d[i].z;
4108 const double dym = ypos + m_ch3d[i].y - 2 * m_coplay;
4110 if (dx * dx + dz * dz >
pow(rcut * 2 * m_sy, 2)) {
4112 exsum = eysum = ezsum = vsum = 0.;
4114 for (
unsigned int j = 1; j <= m_nTermBessel; ++j) {
4116 const double rr = Pi * j *
sqrt(dx * dx + dz * dz) / m_sy;
4117 const double zzp = Pi * j * dy / m_sy;
4118 const double zzn = Pi * j * dym / m_sy;
4128 const double czzp =
cos(zzp);
4129 const double czzn =
cos(zzn);
4130 vsum += (1. / m_sy) * k0r * (czzp - czzn);
4131 err = (TwoPi * j / (m_sy * m_sy)) * k1r * (czzp - czzn);
4132 ezz = (TwoPi * j / (m_sy * m_sy)) * k0r * (
sin(zzp) -
sin(zzn));
4133 exsum += err * dx /
sqrt(dx * dx + dz * dz);
4134 ezsum += err * dz /
sqrt(dx * dx + dz * dz);
4141 for (
unsigned int j = 0; j <= m_nTermPoly; ++j) {
4143 rr1 =
sqrt(dx * dx + dz * dz +
pow(dy + j * 2 * m_sy, 2));
4144 rr2 =
sqrt(dx * dx + dz * dz +
pow(dy - j * 2 * m_sy, 2));
4145 rm1 =
sqrt(dx * dx + dz * dz +
pow(dym - j * 2 * m_sy, 2));
4146 rm2 =
sqrt(dx * dx + dz * dz +
pow(dym + j * 2 * m_sy, 2));
4147 const double rr13 =
pow(rr1, 3);
4148 const double rm13 =
pow(rm1, 3);
4151 vsum = 1. / rr1 - 1. / rm1;
4152 exsum = dx * (1. / rr13 - 1. / rm13);
4153 ezsum = dz * (1. / rr13 - 1. / rm13);
4154 eysum = dy / rr13 - dym / rm13;
4158 const double rr23 =
pow(rr2, 3);
4159 const double rm23 =
pow(rm2, 3);
4160 vsum += 1. / rr1 + 1. / rr2 - 1. / rm1 - 1. / rm2;
4161 exsum += dx * (1. / rr13 + 1. / rr23 - 1. / rm13 - 1. / rm23);
4162 ezsum += dz * (1. / rr13 + 1. / rr23 - 1. / rm13 - 1. / rm23);
4163 eysum += (dy + j * 2 * m_sy) / rr13 + (dy - j * 2 * m_sy) / rr23 -
4164 (dym - j * 2 * m_sy) / rm13 - (dym + j * 2 * m_sy) / rm23;
4169 const double dxm = xpos + m_ch3d[i].x - 2. * m_coplax;
4170 if (dxm * dxm + dz * dz >
pow(rcut * 2 * m_sy, 2)) {
4173 for (
unsigned int j = 1; j <= m_nTermBessel; ++j) {
4175 const double rrm = Pi * j *
sqrt(dxm * dxm + dz * dz) / m_sy;
4176 const double zzp = Pi * j * dy / m_sy;
4177 const double zzn = Pi * j * dym / m_sy;
4187 const double czzp =
cos(zzp);
4188 const double czzn =
cos(zzn);
4189 vsum += (1. / m_sy) * k0rm * (czzp - czzn);
4190 err = (TwoPi / (m_sy * m_sy)) * k1rm * (czzp - czzn);
4191 ezz = (TwoPi / (m_sy * m_sy)) * k0rm * (
sin(zzp) -
sin(zzn));
4192 exsum += err * dxm /
sqrt(dxm * dxm + dz * dz);
4193 ezsum += err * dz /
sqrt(dxm * dxm + dz * dz);
4199 for (
unsigned int j = 0; j <= m_nTermPoly; ++j) {
4201 rr1 =
sqrt(
pow(dy + j * 2 * m_sy, 2) + dxm * dxm + dz * dz);
4202 rr2 =
sqrt(
pow(dy - j * 2 * m_sy, 2) + dxm * dxm + dz * dz);
4203 rm1 =
sqrt(
pow(dym - j * 2 * m_sy, 2) + dxm * dxm + dz * dz);
4204 rm2 =
sqrt(
pow(dym + j * 2 * m_sy, 2) + dxm * dxm + dz * dz);
4205 const double rr13 =
pow(rr1, 3);
4206 const double rm13 =
pow(rm1, 3);
4209 vsum += -1. / rr1 + 1. / rm1;
4210 exsum += -dxm * (1. / rr13 - 1. / rm13);
4211 ezsum += -dz * (1. / rr13 - 1. / rm13);
4212 eysum += -dy / rr13 + dym / rm13;
4215 const double rr23 =
pow(rr2, 3);
4216 const double rm23 =
pow(rm2, 3);
4218 vsum += -1. / rr1 - 1. / rr2 + 1. / rm1 + 1. / rm2;
4219 exsum += -dxm * (1. / rr13 + 1. / rr23 - 1. / rm13 - 1. / rm23);
4220 ezsum += -dz * (1. / rr13 + 1. / rr23 - 1. / rm13 - 1. / rm23);
4221 eysum += -(dy + j * 2 * m_sy) / rr13 - (dy - j * 2 * m_sy) / rr23 +
4222 (dym - j * 2 * m_sy) / rm13 + (dym + j * 2 * m_sy) / rm23;
4226 ex += m_ch3d[i].e * exsum;
4227 ey += m_ch3d[i].e * eysum;
4228 ez += m_ch3d[i].e * ezsum;
4229 volt += m_ch3d[i].e * vsum;
4233void ComponentAnalyticField::Field3dD10(
const double xxpos,
const double yypos,
4234 const double zzpos,
double& eex,
4235 double& eey,
double& eez,
4247 const double rcut = 1.;
4249 double x3d, y3d, z3d;
4250 double exsum = 0., eysum = 0., ezsum = 0., vsum = 0.;
4251 double rr1, rr2, rm1, rm2, err, ezz;
4255 eex = eey = eez = volt = 0.;
4256 double ex = 0., ey = 0., ez = 0.;
4261 std::cerr <<
" Inappropriate potential function.\n";
4266 const double ssx = log(2. * m_cotube / m_w[0].d);
4267 const double cpl = log(0.5 * m_w[0].d);
4270 const double xpos = 0.5 * log(xxpos * xxpos + yypos * yypos);
4271 const double ypos = atan2(yypos, xxpos);
4272 const double zpos = zzpos;
4275 const unsigned int n3d = m_ch3d.size();
4276 for (
unsigned int i = 0; i < n3d; ++i) {
4277 for (
int ii = -1; ii <= 1; ++ii) {
4278 x3d = 0.5 * log(m_ch3d[i].x * m_ch3d[i].x + m_ch3d[i].y * m_ch3d[i].y);
4279 y3d = atan2(m_ch3d[i].y, m_ch3d[i].x + ii * TwoPi);
4281 const double dx = xpos - x3d;
4282 const double dy = ypos - y3d;
4283 const double dz = zpos - z3d;
4284 const double dxm = xpos + x3d - 2 * cpl;
4286 if (xpos == x3d && ypos == y3d && zpos == z3d)
continue;
4288 if (dy * dy + dz * dz >
pow(rcut * 2 * ssx, 2)) {
4290 exsum = eysum = ezsum = vsum = 0.;
4292 for (
unsigned int j = 1; j <= m_nTermBessel; ++j) {
4294 const double rr = Pi * j *
sqrt(dy * dy + dz * dz) / ssx;
4295 const double zzp = Pi * j * dx / ssx;
4296 const double zzn = Pi * j * dxm / ssx;
4306 const double czzp =
cos(zzp);
4307 const double czzn =
cos(zzn);
4308 vsum += (1. / ssx) * k0r * (czzp - czzn);
4309 err = (j * TwoPi / (ssx * ssx)) * k1r * (czzp - czzn);
4310 ezz = (j * TwoPi / (ssx * ssx)) * k0r * (
sin(zzp) -
sin(zzn));
4312 eysum += err * dy /
sqrt(dy * dy + dz * dz);
4313 ezsum += err * dz /
sqrt(dy * dy + dz * dz);
4319 for (
unsigned int j = 0; j < m_nTermPoly; ++j) {
4321 rr1 =
sqrt(
pow(dx + j * 2 * ssx, 2) + dy * dy + dz * dz);
4322 rr2 =
sqrt(
pow(dx - j * 2 * ssx, 2) + dy * dy + dz * dz);
4323 rm1 =
sqrt(
pow(dxm - j * 2 * ssx, 2) + dy * dy + dz * dz);
4324 rm2 =
sqrt(
pow(dxm + j * 2 * ssx, 2) + dy * dy + dz * dz);
4325 const double rr13 =
pow(rr1, 3);
4326 const double rm13 =
pow(rm1, 3);
4329 vsum = 1. / rr1 - 1. / rm1;
4330 exsum = dxm / rr13 - dxm / rm13;
4331 eysum = dy * (1. / rr13 - 1. / rm13);
4332 ezsum = dz * (1. / rr13 - 1. / rm13);
4335 const double rr23 =
pow(rr2, 3);
4336 const double rm23 =
pow(rm2, 3);
4338 vsum += 1. / rr1 + 1. / rr2 - 1. / rm1 - 1. / rm2;
4339 exsum += (dx + j * 2 * ssx) / rr13 + (dx - j * 2 * ssx) / rr23 -
4340 (dxm - j * 2 * ssx) / rm13 - (dxm + j * 2 * ssx) / rm23;
4341 eysum += dy * (1. / rr13 + 1. / rr23 - 1. / rm13 - 1. / rm23);
4342 ezsum += dz * (1. / rr13 + 1. / rr23 - 1. / rm13 - 1. / rm23);
4345 ex += m_ch3d[i].e * exsum;
4346 ey += m_ch3d[i].e * eysum;
4347 ez += m_ch3d[i].e * ezsum;
4353 eex =
exp(-xpos) * (ex *
cos(ypos) - ey *
sin(ypos));
4354 eey =
exp(-ypos) * (ex *
sin(ypos) + ey *
cos(ypos));
4358bool ComponentAnalyticField::PrepareSignals() {
4365 if (m_readout.empty()) {
4366 std::cerr <<
m_className <<
"::PrepareSignals:\n";
4367 std::cerr <<
" There are no readout groups defined.\n";
4368 std::cerr <<
" Calculation of weighting fields makes no sense.\n";
4374 std::cerr <<
m_className <<
"::PrepareSignals:\n";
4375 std::cerr <<
" Cell could not be set up.\n";
4376 std::cerr <<
" No calculation of weighting fields possible.\n";
4383 if (nFourier == 0) {
4384 m_scellTypeFourier = m_scellType;
4385 }
else if (m_scellType ==
"A " || m_scellType ==
"B1X" ||
4386 m_scellType ==
"B1Y" || m_scellType ==
"C1 ") {
4387 m_scellTypeFourier =
"A ";
4388 }
else if (m_scellType ==
"B2X" || m_scellType ==
"C2X") {
4389 m_scellTypeFourier =
"B2X";
4390 }
else if (m_scellType ==
"B2Y" || m_scellType ==
"C2Y") {
4391 m_scellTypeFourier =
"B2Y";
4392 }
else if (m_scellType ==
"C3 ") {
4393 m_scellTypeFourier =
"C3 ";
4394 }
else if (m_scellType ==
"D1 ") {
4395 m_scellTypeFourier =
"D1 ";
4396 }
else if (m_scellType ==
"D3 ") {
4397 m_scellTypeFourier =
"D3 ";
4400 std::cerr <<
m_className <<
"::PrepareSignals:\n";
4401 std::cerr <<
" No potentials available to handle cell type "
4402 << m_scellType <<
".\n";
4407 fperx = fpery =
false;
4408 if (nFourier == 0) {
4411 if (m_scellType ==
"B1X" || m_scellType ==
"C1 " || m_scellType ==
"C2Y") {
4414 if (m_scellType ==
"B1Y" || m_scellType ==
"C1 " || m_scellType ==
"C2X") {
4417 mfexp = int(0.1 + log(1. * nFourier) / log(2.));
4419 fperx = fpery =
false;
4423 mxmin = mymin = mxmax = mymax = 0;
4425 mxmin = std::min(0, 1 - nFourier / 2);
4426 mxmax = nFourier / 2;
4429 mymin = std::min(0, 1 - nFourier / 2);
4430 mymax = nFourier / 2;
4435 std::cout <<
m_className <<
"::PrepareSignals:\n";
4436 std::cout <<
" Cell type: " << m_scellType <<
"\n";
4437 std::cout <<
" Fourier cell type: " << m_scellTypeFourier <<
"\n";
4438 std::cout <<
" x convolutions: " << fperx <<
"\n";
4439 std::cout <<
" y convolutions: " << fpery <<
"\n";
4440 std::cout <<
" No of Fourier terms: " << nFourier <<
" (= 2**" << mfexp
4445 if (!SetupWireSignals()) {
4446 std::cerr <<
m_className <<
"::PrepareSignals:\n";
4447 std::cerr <<
" Preparing wire signal capacitance matrices failed.\n";
4451 if (!SetupPlaneSignals()) {
4452 std::cerr <<
m_className <<
"::PrepareSignals:\n";
4453 std::cerr <<
" Preparing plane charges failed.\n";
4463 const unsigned int nReadout = m_readout.size();
4464 for (
unsigned int i = 0; i < nReadout; ++i) {
4465 for (
unsigned int j = 0; j < m_nWires; ++j) {
4466 if (m_w[j].type == m_readout[i]) m_w[j].ind = i;
4468 for (
unsigned int j = 0; j < 5; ++j) {
4469 if (planes[j].type == m_readout[i]) planes[j].ind = i;
4470 const unsigned int nStrips1 = planes[j].strips1.size();
4471 for (
unsigned int k = 0; k < nStrips1; ++k) {
4472 if (planes[j].strips1[k].type == m_readout[i]) {
4473 planes[j].strips1[k].ind = i;
4476 const unsigned int nStrips2 = planes[j].strips2.size();
4477 for (
unsigned int k = 0; k < nStrips2; ++k) {
4478 if (planes[j].strips2[k].type == m_readout[i]) {
4479 planes[j].strips2[k].ind = i;
4482 const unsigned int nPixels = planes[j].pixels.size();
4483 for (
unsigned int k = 0; k < nPixels; ++k) {
4484 if (planes[j].pixels[k].type == m_readout[i]) {
4485 planes[j].pixels[k].ind = i;
4496bool ComponentAnalyticField::SetupWireSignals() {
4510 m_sigmat.resize(m_nWires);
4511 for (
unsigned int i = 0; i < m_nWires; ++i) {
4512 m_sigmat[i].clear();
4513 m_sigmat[i].resize(m_nWires);
4516 std::vector<std::vector<std::complex<double> > > fftmat;
4519 if (fperx || fpery) {
4520 fftmat.resize(nFourier);
4521 for (
int i = 0; i < nFourier; ++i) {
4522 fftmat[i].resize(m_nWires);
4526 if (fperx || fpery) {
4537 for (
int mx = mxmin; mx <= mxmax; ++mx) {
4538 for (
int my = mymin; my <= mymax; ++my) {
4540 if (m_scellTypeFourier ==
"A ") {
4542 }
else if (m_scellTypeFourier ==
"B2X") {
4544 }
else if (m_scellTypeFourier ==
"B2Y") {
4546 }
else if (m_scellTypeFourier ==
"C2X") {
4548 }
else if (m_scellTypeFourier ==
"C2Y") {
4550 }
else if (m_scellTypeFourier ==
"C3 ") {
4552 }
else if (m_scellTypeFourier ==
"D1 ") {
4554 }
else if (m_scellTypeFourier ==
"D3 ") {
4557 std::cerr <<
m_className <<
"::SetupWireSignals:\n";
4558 std::cerr <<
" Unknown signal cell type " << m_scellTypeFourier
4563 std::cout <<
m_className <<
"::SetupWireSignals:\n";
4564 std::cout <<
" Signal matrix MX = " << mx <<
", MY = " << my
4565 <<
" has been calculated.\n";
4567 if (fperx || fpery) {
4575 std::cout <<
m_className <<
"::SetupWireSignals:\n";
4576 std::cout <<
" Dump of signal matrix (" << mx <<
", " << my
4577 <<
") before inversion:\n";
4578 for (
unsigned int i = 0; i < m_nWires; i += 10) {
4579 for (
unsigned int j = 0; j < m_nWires; j += 10) {
4580 std::cout <<
" (Re-Block " << i / 10 <<
", " << j / 10 <<
")\n";
4581 for (
unsigned int ii = 0; ii < 10; ++ii) {
4582 if (i + ii >= m_nWires)
break;
4583 for (
unsigned int jj = 0; jj < 10; ++jj) {
4584 if (j + jj >= m_nWires)
break;
4585 std::cout << real(m_sigmat[i + ii][j + jj]) <<
" ";
4590 std::cout <<
" (Im-Block " << i / 10 <<
", " << j / 10 <<
")\n";
4591 for (
unsigned int ii = 0; ii < 10; ++ii) {
4592 if (i + ii >= m_nWires)
break;
4593 for (
unsigned int jj = 0; jj < 10; ++jj) {
4594 if (j + jj >= m_nWires)
break;
4595 std::cout << imag(m_sigmat[i + ii][j + jj]) <<
" ";
4602 std::cout <<
m_className <<
"::SetupWireSignals:\n";
4603 std::cout <<
" End of the uninverted capacitance matrix dump.\n";
4610 if ((fperx && !fpery) || (fpery && !fperx)) {
4611 for (
unsigned int i = 0; i < m_nWires; ++i) {
4612 for (
int m = -nFourier / 2; m < nFourier / 2; ++m) {
4615 for (
unsigned int j = 0; j < m_nWires; ++j) {
4616 fftmat[m + nFourier / 2][j] = m_sigmat[i][j];
4619 for (
unsigned int j = 0; j < m_nWires; ++j) {
4622 for (
int m = -nFourier / 2; m < nFourier / 2; ++m) {
4625 for (
unsigned int j = 0; j < m_nWires; ++j) {
4626 m_sigmat[i][j] = fftmat[m + nFourier / 2][j];
4634 if (fperx || fpery) {
4635 for (
unsigned int i = 0; i < m_nWires; ++i) {
4636 for (
int mx = mxmin; mx <= mxmax; ++mx) {
4637 for (
int my = mymin; my <= mymax; ++my) {
4640 for (
unsigned int j = 0; j < m_nWires; ++j) {
4641 fftmat[my + nFourier / 2 - 1][j] = m_sigmat[i][j];
4644 for (
unsigned int j = 0; j < m_nWires; ++j) {
4647 for (
int my = mymin; my <= mymax; ++my) {
4650 for (
unsigned int j = 0; j < m_nWires; ++j) {
4651 m_sigmat[i][j] = fftmat[my + nFourier / 2 - 1][j];
4657 for (
int my = mymin; my <= mymax; ++my) {
4658 for (
int mx = mxmin; mx <= mxmax; ++mx) {
4661 for (
unsigned int j = 0; j < m_nWires; ++j) {
4662 fftmat[mx + nFourier / 2 - 1][j] = m_sigmat[i][j];
4665 for (
unsigned int j = 0; j < m_nWires; ++j) {
4668 for (
int mx = mxmin; mx <= mxmax; ++mx) {
4671 for (
unsigned int j = 0; j < m_nWires; ++j) {
4672 m_sigmat[i][j] = fftmat[mx + nFourier / 2 - 1][j];
4682 for (
int mx = mxmin; mx <= mxmax; ++mx) {
4683 for (
int my = mymin; my <= mymax; ++my) {
4685 if (fperx || fpery) {
4690 if (m_nWires >= 1) {
4694 std::cerr <<
m_className <<
"::PrepareWireSignals:\n";
4695 std::cerr <<
" Inversion of signal matrix (" << mx <<
", " << my
4697 std::cerr <<
" No reliable results.\n";
4698 std::cerr <<
" Preparation of weighting fields is abandoned.\n";
4703 if (fperx || fpery) {
4712 if ((fperx && !fpery) || (fpery && !fperx)) {
4713 for (
unsigned int i = 0; i < m_nWires; ++i) {
4714 for (
int m = -nFourier / 2; m < nFourier / 2; ++m) {
4717 for (
unsigned int j = 0; j < m_nWires; ++j) {
4718 fftmat[m + nFourier / 2][j] = m_sigmat[i][j];
4721 for (
unsigned int j = 0; j < m_nWires; ++j) {
4724 for (
int m = -nFourier / 2; m < nFourier / 2; ++m) {
4727 for (
unsigned int j = 0; j < m_nWires; ++j) {
4728 m_sigmat[i][j] = fftmat[m + nFourier / 2][j] / double(nFourier);
4736 if (fperx && fpery) {
4737 for (
unsigned int i = 0; i < m_nWires; ++i) {
4738 for (
int mx = mxmin; mx <= mxmax; ++mx) {
4739 for (
int my = mymin; my <= mymax; ++my) {
4742 for (
unsigned int j = 0; j < m_nWires; ++j) {
4743 fftmat[my + nFourier / 2 - 1][j] = m_sigmat[i][j];
4746 for (
unsigned int j = 0; j < m_nWires; ++j) {
4749 for (
int my = mymin; my <= mymax; ++my) {
4752 for (
unsigned int j = 0; j < m_nWires; ++j) {
4754 fftmat[my + nFourier / 2 - 1][j] / double(nFourier);
4760 for (
int my = mymin; my <= mymax; ++my) {
4761 for (
int mx = mxmin; mx <= mxmax; ++mx) {
4764 for (
unsigned int j = 0; j < m_nWires; ++j) {
4765 fftmat[mx + nFourier / 2 - 1][j] = m_sigmat[i][j];
4768 for (
unsigned int j = 0; j < m_nWires; ++j) {
4771 for (
int mx = mxmin; mx <= mxmax; ++mx) {
4774 for (
unsigned int j = 0; j < m_nWires; ++j) {
4776 fftmat[mx + nFourier / 2 - 1][j] / double(nFourier);
4787 for (
int mx = mxmin; mx <= mxmax; ++mx) {
4788 for (
int my = mymin; my <= mymax; ++my) {
4789 std::cout <<
m_className <<
"::SetupWireSignals:\n";
4790 std::cout <<
" Dump of signal matrix (" << mx <<
", " << my
4791 <<
") after inversion:\n";
4792 for (
unsigned int i = 0; i < m_nWires; i += 10) {
4793 for (
unsigned int j = 0; j < m_nWires; j += 10) {
4794 std::cout <<
" (Re-Block " << i / 10 <<
", " << j / 10 <<
")\n";
4795 for (
unsigned int ii = 0; ii < 10; ++ii) {
4796 if (i + ii >= m_nWires)
break;
4797 for (
unsigned int jj = 0; jj < 10; ++jj) {
4798 if (j + jj >= m_nWires)
break;
4799 std::cout << real(m_sigmat[i + ii][j + jj]) <<
" ";
4804 std::cout <<
" (Im-Block " << i / 10 <<
", " << j / 10 <<
")\n";
4805 for (
int ii = 0; ii < 10; ++ii) {
4806 if (i + ii >= m_nWires)
break;
4807 for (
int jj = 0; jj < 10; ++jj) {
4808 if (j + jj >= m_nWires)
break;
4809 std::cout << imag(m_sigmat[i + ii][j + jj]) <<
" ";
4816 std::cout <<
m_className <<
"::SetupWireSignals:\n";
4817 std::cout <<
" End of the inverted capacitance matrix dump.\n";
4824bool ComponentAnalyticField::SetupPlaneSignals() {
4832 const int nPlanes = 5;
4833 m_qplane.assign(nPlanes, std::vector<double>(m_nWires, 0.));
4838 for (
int mx = mxmin; mx <= mxmax; ++mx) {
4839 for (
int my = mymin; my <= mymax; ++my) {
4848 m_qplane.assign(nPlanes, std::vector<double>(m_nWires, 0.));
4852 for (
unsigned int i = 0; i < m_nWires; ++i) {
4854 vw = -(m_coplan[1] - m_w[i].x) / (m_coplan[1] - m_coplan[0]);
4855 }
else if (m_perx) {
4856 vw = -(m_coplan[0] + m_sx - m_w[i].x) / m_sx;
4861 for (
unsigned int j = 0; j < m_nWires; ++j) {
4862 m_qplane[0][j] += real(m_sigmat[i][j]) * vw;
4869 for (
unsigned int i = 0; i < m_nWires; ++i) {
4871 vw = -(m_coplan[0] - m_w[i].x) / (m_coplan[0] - m_coplan[1]);
4872 }
else if (m_perx) {
4873 vw = -(m_w[i].x - m_coplan[1] + m_sx) / m_sx;
4878 for (
unsigned int j = 0; j < m_nWires; ++j) {
4879 m_qplane[1][j] += real(m_sigmat[i][j]) * vw;
4886 for (
unsigned int i = 0; i < m_nWires; ++i) {
4888 vw = -(m_coplan[3] - m_w[i].y) / (m_coplan[3] - m_coplan[2]);
4889 }
else if (m_pery) {
4890 vw = -(m_coplan[2] + m_sy - m_w[i].y) / m_sy;
4895 for (
unsigned int j = 0; j < m_nWires; ++j) {
4896 m_qplane[2][i] += real(m_sigmat[i][j]) * vw;
4903 for (
unsigned int i = 0; i < m_nWires; ++i) {
4905 vw = -(m_coplan[2] - m_w[i].y) / (m_coplan[2] - m_coplan[3]);
4906 }
else if (m_pery) {
4907 vw = -(m_w[i].y - m_coplan[3] + m_sy) / m_sy;
4912 for (
unsigned int j = 0; j < m_nWires; ++j) {
4913 m_qplane[3][i] += real(m_sigmat[i][j]) * vw;
4919 for (
unsigned int i = 0; i < m_nWires; ++i) {
4920 for (
unsigned int j = 0; j < m_nWires; ++j) {
4921 m_qplane[4][i] -= real(m_sigmat[i][j]);
4936 if (m_ynplan[0] && m_ynplan[1]) {
4937 planes[0].ewxcor = 1. / (m_coplan[1] - m_coplan[0]);
4938 planes[1].ewxcor = 1. / (m_coplan[0] - m_coplan[1]);
4939 }
else if (m_ynplan[0] && m_perx) {
4940 planes[0].ewxcor = 1. / m_sx;
4941 planes[1].ewxcor = 0.;
4942 }
else if (m_ynplan[1] && m_perx) {
4943 planes[0].ewxcor = 0.;
4944 planes[1].ewxcor = -1. / m_sx;
4946 planes[0].ewxcor = planes[1].ewxcor = 0.;
4948 planes[2].ewxcor = planes[3].ewxcor = planes[4].ewxcor = 0.;
4950 planes[0].ewycor = planes[1].ewycor = 0.;
4951 if (m_ynplan[2] && m_ynplan[3]) {
4952 planes[2].ewycor = 1. / (m_coplan[3] - m_coplan[2]);
4953 planes[3].ewycor = 1. / (m_coplan[2] - m_coplan[3]);
4954 }
else if (m_ynplan[2] && m_pery) {
4955 planes[2].ewycor = 1. / m_sy;
4956 planes[3].ewycor = 0.;
4957 }
else if (m_ynplan[3] && m_pery) {
4958 planes[2].ewycor = 0.;
4959 planes[3].ewycor = -1. / m_sy;
4961 planes[2].ewycor = planes[3].ewycor = 0.;
4964 planes[4].ewycor = 0.;
4968 std::cout <<
m_className <<
"::SetupPlaneSignals:\n";
4969 std::cout <<
" Charges for currents induced in the planes:\n";
4970 std::cout <<
" Wire x-Plane 1 x-Plane 2"
4971 <<
" y-Plane 1 y-Plane 2"
4973 for (
unsigned int i = 0; i < m_nWires; ++i) {
4974 std::cout <<
" " << i <<
" " << m_qplane[0][i] <<
" "
4975 << m_qplane[1][i] <<
" " << m_qplane[2][i] <<
" "
4976 << m_qplane[3][i] <<
" " << m_qplane[4][i] <<
"\n";
4978 std::cout <<
m_className <<
"::SetupPlaneSignals:\n";
4979 std::cout <<
" Bias fields:\n";
4980 std::cout <<
" Plane x-Bias [1/cm] y-Bias [1/cm]\n";
4981 for (
int i = 0; i < 4; ++i) {
4982 std::cout <<
" " << i <<
" " << planes[i].ewxcor <<
" "
4983 << planes[i].ewycor <<
"\n";
4990bool ComponentAnalyticField::IprA00(
const int mx,
const int my) {
4997 const double dx = mx * m_sx;
4998 const double dy = my * m_sy;
5001 for (
unsigned int i = 0; i < m_nWires; ++i) {
5003 if (dx != 0. || dy != 0.) {
5004 aa = dx * dx + dy * dy;
5006 aa = 0.25 * m_w[i].d * m_w[i].d;
5009 if (m_ynplax) aa /= 2. *
pow(m_w[i].x - m_coplax, 2) + dy * dy;
5010 if (m_ynplay) aa /= 2. *
pow(m_w[i].y - m_coplay, 2) + dx * dx;
5012 if (m_ynplax && m_ynplay)
5013 aa *= 4. * (
pow(m_w[i].x - m_coplax, 2) +
pow(m_w[i].y - m_coplay, 2));
5015 m_sigmat[i][i] = -0.5 * log(aa);
5016 for (
unsigned int j = i + 1; j < m_nWires; ++j) {
5017 aa =
pow(m_w[i].x + dx - m_w[j].x, 2) +
pow(m_w[i].y + dy - m_w[j].y, 2);
5020 aa /=
pow(2. * m_coplax - m_w[i].x - dx - m_w[j].x, 2) +
5021 pow(m_w[i].y + dy - m_w[j].y, 2);
5023 aa /=
pow(m_w[i].x + dx - m_w[j].x, 2) +
5024 pow(2. * m_coplay - m_w[i].y - dy - m_w[j].y, 2);
5026 if (m_ynplax && m_ynplay) {
5027 aa *=
pow(2. * m_coplax - m_w[i].x - dx - m_w[j].x, 2) +
5028 pow(2. * m_coplay - m_w[i].y - dy - m_w[j].y, 2);
5031 m_sigmat[i][j] = -0.5 * log(aa);
5032 m_sigmat[j][i] = m_sigmat[i][j];
5038bool ComponentAnalyticField::IprB2X(
const int my) {
5046 m_b2sin.resize(m_nWires);
5048 const double dy = my * m_sy;
5053 for (
unsigned int i = 0; i < m_nWires; ++i) {
5054 double xx = (Pi / m_sx) * (m_w[i].x - m_coplan[0]);
5056 aa =
pow(sinh(Pi * dy / m_sx) /
sin(xx), 2);
5058 aa =
pow((0.25 * m_w[i].d * Pi / m_sx) /
sin(xx), 2);
5062 const double yymirr = (Pi / m_sx) * (m_w[i].y - m_coplay);
5063 if (
fabs(yymirr) <= 20.) {
5064 const double sinhy = sinh(yymirr);
5065 const double sinxx =
sin(xx);
5066 aa *= (sinhy * sinhy + sinxx * sinxx) / (sinhy * sinhy);
5070 m_sigmat[i][i] = -0.5 * log(aa);
5072 for (
unsigned int j = i + 1; j < m_nWires; ++j) {
5073 const double yy = HalfPi * (m_w[i].y + dy - m_w[j].y) / m_sx;
5074 xx = HalfPi * (m_w[i].x - m_w[j].x) / m_sx;
5075 xxneg = HalfPi * (m_w[i].x + m_w[j].x - 2. * m_coplan[0]) / m_sx;
5076 if (
fabs(yy) <= 20.) {
5077 const double sinhy = sinh(yy);
5078 const double sinxx =
sin(xx);
5079 const double sinxxneg =
sin(xxneg);
5080 aa = (sinhy * sinhy + sinxx * sinxx) /
5081 (sinhy * sinhy + sinxxneg * sinxxneg);
5087 const double yymirr =
5088 HalfPi * (m_w[i].y + m_w[j].y - 2. * m_coplay) / m_sx;
5089 if (
fabs(yymirr) <= 20.) {
5090 const double sinhy = sinh(yymirr);
5091 const double sinxx =
sin(xx);
5092 const double sinxxneg =
sin(xxneg);
5093 aa *= (sinhy * sinhy + sinxxneg * sinxxneg) /
5094 (sinhy * sinhy + sinxx * sinxx);
5098 m_sigmat[i][j] = -0.5 * log(aa);
5099 m_sigmat[j][i] = m_sigmat[i][j];
5102 m_b2sin[i] =
sin(Pi * (m_coplan[0] - m_w[i].x) / m_sx);
5108bool ComponentAnalyticField::IprB2Y(
const int mx) {
5116 m_b2sin.resize(m_nWires);
5118 const double dx = mx * m_sx;
5120 double xx, yy, xxmirr, yyneg;
5123 for (
unsigned int i = 0; i < m_nWires; ++i) {
5124 yy = (Pi / m_sy) * (m_w[i].y - m_coplan[2]);
5126 aa =
pow(sinh(Pi * dx / m_sy) /
sin(yy), 2);
5128 aa =
pow((0.25 * m_w[i].d * Pi / m_sy) /
sin(yy), 2);
5132 xxmirr = (Pi / m_sy) * (m_w[i].x - m_coplax);
5133 if (
fabs(xxmirr) <= 20.) {
5134 aa *= (
pow(sinh(xxmirr), 2) +
pow(
sin(yy), 2)) /
pow(sinh(xxmirr), 2);
5138 m_sigmat[i][i] = -0.5 * log(aa);
5140 for (
unsigned int j = i + 1; j < m_nWires; ++j) {
5141 xx = HalfPi * (m_w[i].x + dx - m_w[j].x) / m_sy;
5142 yy = HalfPi * (m_w[i].y - m_w[j].y) / m_sy;
5143 yyneg = HalfPi * (m_w[i].y + m_w[j].y - 2. * m_coplan[2]) / m_sy;
5144 if (
fabs(xx) <= 20.) {
5145 aa = (
pow(sinh(xx), 2) +
pow(
sin(yy), 2)) /
5146 (
pow(sinh(xx), 2) +
pow(
sin(yyneg), 2));
5152 xxmirr = HalfPi * (m_w[i].x + m_w[j].x - 2. * m_coplax) / m_sy;
5153 if (
fabs(xxmirr) <= 20.) {
5154 aa *= (
pow(sinh(xxmirr), 2) +
pow(
sin(yyneg), 2)) /
5155 (
pow(sinh(xxmirr), 2) +
pow(
sin(yy), 2));
5159 m_sigmat[i][j] = -0.5 * log(aa);
5160 m_sigmat[j][i] = m_sigmat[i][j];
5163 m_b2sin[i] =
sin(Pi * (m_coplan[2] - m_w[i].y) / m_sy);
5168bool ComponentAnalyticField::IprC2X() {
5180 for (
unsigned int i = 0; i < m_nWires; ++i) {
5181 double cx = m_coplax - m_sx * int(round((m_coplax - m_w[i].x) / m_sx));
5182 for (
unsigned int j = 0; j < m_nWires; ++j) {
5185 temp = (m_w[i].x - cx) * (m_w[j].x - cx) * TwoPi / (m_sx * m_sy);
5189 Ph2Lim(0.5 * m_w[i].d) - Ph2(2. * (m_w[j].x - cx), 0.) - temp;
5192 Ph2(m_w[i].x - m_w[j].x, m_w[i].y - m_w[j].y) -
5193 Ph2(m_w[i].x + m_w[j].x - 2. * cx, m_w[i].y - m_w[j].y) - temp;
5200bool ComponentAnalyticField::IprC2Y() {
5212 for (
unsigned int i = 0; i < m_nWires; ++i) {
5213 const double cy = m_coplay - m_sy * int(round((m_coplay - m_w[i].y) / m_sy));
5214 for (
unsigned int j = 0; j < m_nWires; ++j) {
5217 temp = (m_w[i].y - cy) * (m_w[j].y - cy) * TwoPi / (m_sx * m_sy);
5221 Ph2Lim(0.5 * m_w[i].d) - Ph2(0., 2. * (m_w[j].y - cy)) - temp;
5224 Ph2(m_w[i].x - m_w[j].x, m_w[i].y - m_w[j].y) -
5225 Ph2(m_w[i].x - m_w[j].x, m_w[i].y + m_w[j].y - 2. * cy) - temp;
5232bool ComponentAnalyticField::IprC30() {
5243 for (
unsigned int i = 0; i < m_nWires; ++i) {
5244 const double cx = m_coplax - m_sx * int(round((m_coplax - m_w[i].x) / m_sx));
5245 const double cy = m_coplay - m_sy * int(round((m_coplay - m_w[i].y) / m_sy));
5246 for (
unsigned int j = 0; j < m_nWires; ++j) {
5248 m_sigmat[i][i] = Ph2Lim(0.5 * m_w[i].d) -
5249 Ph2(0., 2. * (m_w[i].y - cy)) -
5250 Ph2(2. * (m_w[i].x - cx), 0.) +
5251 Ph2(2. * (m_w[i].x - cx), 2. * (m_w[i].y - cy));
5254 Ph2(m_w[i].x - m_w[j].x, m_w[i].y - m_w[j].y) -
5255 Ph2(m_w[i].x - m_w[j].x, m_w[i].y + m_w[j].y - 2. * cy) -
5256 Ph2(m_w[i].x + m_w[j].x - 2. * cx, m_w[i].y - m_w[j].y) +
5257 Ph2(m_w[i].x + m_w[j].x - 2. * cx, m_w[i].y + m_w[j].y - 2. * cy);
5264bool ComponentAnalyticField::IprD10() {
5273 for (
unsigned int i = 0; i < m_nWires; ++i) {
5276 -log(0.5 * m_w[i].d /
5278 (m_w[i].x * m_w[i].x + m_w[i].y * m_w[i].y) / m_cotube));
5280 std::complex<double> zi(m_w[i].x, m_w[i].y);
5282 for (
unsigned int j = i + 1; j < m_nWires; ++j) {
5284 std::complex<double> zj(m_w[j].x, m_w[j].y);
5286 -log(abs((1. / m_cotube) * (zi - zj) / (1. - conj(zi) * zj / m_cotube2)));
5288 m_sigmat[j][i] = m_sigmat[i][j];
5294bool ComponentAnalyticField::IprD30() {
5302 wmap.resize(m_nWires);
5304 std::complex<double> wd;
5305 InitializeCoefficientTables();
5308 for (
int i = 0; i < int(m_nWires); ++i) {
5310 ConformalMap(std::complex<double>(m_w[i].x, m_w[i].y) / m_cotube, wmap[i],
5313 m_sigmat[i][i] = -log(abs((0.5 * m_w[i].d / m_cotube) * wd /
5314 (1. -
pow(abs(wmap[i]), 2))));
5316 for (
int j = 0; j < i - 1; ++j) {
5318 -log(abs((wmap[i] - wmap[j]) / (1. - conj(wmap[i]) * wmap[j])));
5320 m_sigmat[j][i] = m_sigmat[i][j];
5326bool ComponentAnalyticField::Wfield(
const double xpos,
const double ypos,
5327 const double zpos,
double& exsum,
5328 double& eysum,
double& ezsum,
double& vsum,
5329 const std::string& label,
5330 const bool opt)
const {
5338 exsum = eysum = ezsum = vsum = 0.;
5339 double ex = 0., ey = 0., ez = 0.;
5343 if (m_readout.empty())
return false;
5346 <<
" No weighting fields available.\n";
5350 if (label.empty())
return volt;
5351 std::vector<std::string>::const_iterator it =
5352 std::find(m_readout.begin(), m_readout.end(), label);
5353 if (it == m_readout.end())
return false;
5354 const int isw = it - m_readout.begin();
5357 for (
int mx = mxmin; mx <= mxmax; ++mx) {
5358 for (
int my = mymin; my <= mymax; ++my) {
5369 for (
int iw = m_nWires; iw--;) {
5371 if (m_w[iw].ind == isw) {
5373 if (m_scellTypeFourier ==
"A ") {
5374 WfieldWireA00(xpos, ypos, ex, ey, volt, mx, my, iw, opt);
5375 }
else if (m_scellTypeFourier ==
"B2X") {
5376 WfieldWireB2X(xpos, ypos, ex, ey, volt, my, iw, opt);
5377 }
else if (m_scellTypeFourier ==
"B2Y") {
5378 WfieldWireB2Y(xpos, ypos, ex, ey, volt, mx, iw, opt);
5379 }
else if (m_scellTypeFourier ==
"C2X") {
5380 WfieldWireC2X(xpos, ypos, ex, ey, volt, iw, opt);
5381 }
else if (m_scellTypeFourier ==
"C2Y") {
5382 WfieldWireC2Y(xpos, ypos, ex, ey, volt, iw, opt);
5383 }
else if (m_scellTypeFourier ==
"C3 ") {
5384 WfieldWireC30(xpos, ypos, ex, ey, volt, iw, opt);
5385 }
else if (m_scellTypeFourier ==
"D1 ") {
5386 WfieldWireD10(xpos, ypos, ex, ey, volt, iw, opt);
5387 }
else if (m_scellTypeFourier ==
"D3 ") {
5388 WfieldWireD30(xpos, ypos, ex, ey, volt, iw, opt);
5391 std::cerr <<
" Unknown signal field type " << m_scellTypeFourier
5392 <<
" received. Program error!\n";
5393 std::cerr <<
" Encountered for wire " << iw
5394 <<
", readout group = " << m_w[iw].ind <<
"\n";
5395 exsum = eysum = ezsum = vsum = 0.;
5401 if (opt) vsum += volt;
5414 for (
int ip = 0; ip < 5; ++ip) {
5416 if (planes[ip].ind == isw) {
5418 if (m_scellTypeFourier ==
"A ") {
5419 WfieldPlaneA00(xpos, ypos, ex, ey, volt, mx, my, ip, opt);
5420 }
else if (m_scellTypeFourier ==
"B2X") {
5421 WfieldPlaneB2X(xpos, ypos, ex, ey, volt, my, ip, opt);
5422 }
else if (m_scellTypeFourier ==
"B2Y") {
5423 WfieldPlaneB2Y(xpos, ypos, ex, ey, volt, mx, ip, opt);
5424 }
else if (m_scellTypeFourier ==
"C2X") {
5425 WfieldPlaneC2X(xpos, ypos, ex, ey, volt, ip, opt);
5426 }
else if (m_scellTypeFourier ==
"C2Y") {
5427 WfieldPlaneC2Y(xpos, ypos, ex, ey, volt, ip, opt);
5428 }
else if (m_scellTypeFourier ==
"D1 ") {
5429 WfieldPlaneD10(xpos, ypos, ex, ey, volt, ip, opt);
5430 }
else if (m_scellTypeFourier ==
"D3 ") {
5431 WfieldPlaneD30(xpos, ypos, ex, ey, volt, ip, opt);
5434 std::cerr <<
" Unkown field type " << m_scellTypeFourier
5435 <<
" received. Program error!\n";
5436 std::cerr <<
" Encountered for plane " << ip
5437 <<
", readout group = " << planes[ip].ind <<
"\n";
5438 exsum = eysum = ezsum = 0.;
5444 if (opt) vsum += volt;
5451 for (
int ip = 0; ip < 5; ++ip) {
5452 if (planes[ip].ind == isw) {
5453 exsum += planes[ip].ewxcor;
5454 eysum += planes[ip].ewycor;
5456 if (ip == 0 || ip == 1) {
5459 xx -= m_sx * int(round(xpos / m_sx));
5460 if (m_ynplan[0] && xx <= m_coplan[0]) xx += m_sx;
5461 if (m_ynplan[1] && xx >= m_coplan[1]) xx -= m_sx;
5463 vsum += 1. - planes[ip].ewxcor * (xx - m_coplan[ip]);
5464 }
else if (ip == 2 || ip == 3) {
5467 yy -= m_sy * int(round(ypos / m_sy));
5468 if (m_ynplan[2] && yy <= m_coplan[2]) yy += m_sy;
5469 if (m_ynplan[3] && yy >= m_coplan[3]) yy -= m_sy;
5471 vsum += 1. - planes[ip].ewycor * (yy - m_coplan[ip]);
5478 for (
unsigned int ip = 0; ip < 5; ++ip) {
5479 const unsigned int nStrips1 = planes[ip].strips1.size();
5480 for (
unsigned int istrip = 0; istrip < nStrips1; ++istrip) {
5481 if (planes[ip].strips1[istrip].ind == isw) {
5482 WfieldStripXy(xpos, ypos, zpos, ex, ey, ez, volt, ip, istrip, opt);
5486 if (opt) vsum += volt;
5489 const unsigned int nStrips2 = planes[ip].strips2.size();
5490 for (
unsigned int istrip = 0; istrip < nStrips2; ++istrip) {
5491 if (planes[ip].strips2[istrip].ind == isw) {
5492 WfieldStripZ(xpos, ypos, ex, ey, volt, ip, istrip, opt);
5495 if (opt) vsum += volt;
5498 const unsigned int nPixels = planes[ip].pixels.size();
5499 for (
unsigned int ipix = 0; ipix < nPixels; ++ipix) {
5500 if (planes[ip].pixels[ipix].ind != isw)
continue;
5501 WfieldPixel(xpos, ypos, zpos, ex, ey, ez, volt, ip, ipix, opt);
5505 if (opt) vsum += volt;
5511void ComponentAnalyticField::WfieldWireA00(
const double xpos,
const double ypos,
5512 double& ex,
double& ey,
double& volt,
5513 const int mx,
const int my,
5514 const int isw,
const bool opt)
const {
5528 ex = ey = volt = 0.;
5530 double xxmirr = 0., yymirr = 0.;
5532 for (
int i = m_nWires; i--;) {
5534 const double xx = xpos - m_w[i].x - mx * m_sx;
5535 const double yy = ypos - m_w[i].y - my * m_sy;
5537 double r2 = xx * xx + yy * yy;
5538 if (r2 <= 0.)
continue;
5539 double exhelp = xx / r2;
5540 double eyhelp = yy / r2;
5543 xxmirr = xpos + m_w[i].x - 2. * m_coplax;
5544 const double r2plan = xxmirr * xxmirr + yy * yy;
5545 if (r2plan <= 0.)
continue;
5546 exhelp -= xxmirr / r2plan;
5547 eyhelp -= yy / r2plan;
5552 yymirr = ypos + m_w[i].y - 2. * m_coplay;
5553 const double r2plan = xx * xx + yymirr * yymirr;
5554 if (r2plan <= 0.)
continue;
5555 exhelp -= xx / r2plan;
5556 eyhelp -= yymirr / r2plan;
5560 if (m_ynplax && m_ynplay) {
5561 const double r2plan = xxmirr * xxmirr + yymirr * yymirr;
5562 if (r2plan <= 0.)
continue;
5563 exhelp += xxmirr / r2plan;
5564 eyhelp += yymirr / r2plan;
5568 if (opt) volt -= 0.5 * real(m_sigmat[isw][i]) * log(r2);
5569 ex += real(m_sigmat[isw][i]) * exhelp;
5570 ey += real(m_sigmat[isw][i]) * eyhelp;
5574void ComponentAnalyticField::WfieldWireB2X(
const double xpos,
const double ypos,
5575 double& ex,
double& ey,
double& volt,
5576 const int my,
const int isw,
5577 const bool opt)
const {
5589 ex = ey = volt = 0.;
5591 const double tx = HalfPi / m_sx;
5593 for (
unsigned int i = 0; i < m_nWires; ++i) {
5594 const double xx = tx * (xpos - m_w[i].x);
5595 const double yy = tx * (ypos - m_w[i].y - my * m_sy);
5596 const double xxneg = tx * (xpos + m_w[i].x - 2 * m_coplan[0]);
5597 const std::complex<double> zz(xx, yy);
5598 const std::complex<double> zzneg(xxneg, yy);
5600 std::complex<double> ecompl(0., 0.);
5602 if (
fabs(yy) <= 20.) {
5603 ecompl = -m_b2sin[i] / (
sin(zz) *
sin(zzneg));
5605 const double sinhy = sinh(yy);
5606 const double sinxx =
sin(xx);
5607 const double sinxxneg =
sin(xxneg);
5608 r2 = (sinhy * sinhy + sinxx * sinxx) /
5609 (sinhy * sinhy + sinxxneg * sinxxneg);
5614 const double yymirr = tx * (ypos + m_w[i].y - 2. * m_coplay);
5615 const std::complex<double> zzmirr(xx, yymirr);
5616 const std::complex<double> zznmirr(xxneg, yymirr);
5617 if (
fabs(yymirr) <= 20.) {
5618 ecompl += m_b2sin[i] / (
sin(zzmirr) *
sin(zznmirr));
5620 const double sinhy = sinh(yymirr);
5621 const double sinxx =
sin(xx);
5622 const double sinxxneg =
sin(xxneg);
5623 const double r2plan = (sinhy * sinhy + sinxx * sinxx) /
5624 (sinhy * sinhy + sinxxneg * sinxxneg);
5630 ex += real(m_sigmat[isw][i]) * real(ecompl);
5631 ey -= real(m_sigmat[isw][i]) * imag(ecompl);
5632 if (opt) volt -= 0.5 * real(m_sigmat[isw][i]) * log(r2);
5638void ComponentAnalyticField::WfieldWireB2Y(
const double xpos,
const double ypos,
5639 double& ex,
double& ey,
double& volt,
5640 const int mx,
const int isw,
5641 const bool opt)
const {
5652 const std::complex<double> icons = std::complex<double>(0., 1.);
5655 ex = ey = volt = 0.;
5657 const double ty = HalfPi / m_sy;
5659 for (
unsigned int i = 0; i < m_nWires; ++i) {
5660 const double xx = ty * (xpos - m_w[i].x - mx * m_sx);
5661 const double yy = ty * (ypos - m_w[i].y);
5662 const double yyneg = ty * (ypos + m_w[i].y - 2. * m_coplan[2]);
5663 const std::complex<double> zz(xx, yy);
5664 const std::complex<double> zzneg(xx, yyneg);
5666 std::complex<double>ecompl(0., 0.);
5668 if (
fabs(xx) <= 20.) {
5669 ecompl = icons * m_b2sin[i] / (
sin(icons * zz) *
sin(icons * zzneg));
5671 const double sinhx = sinh(xx);
5672 const double sinyy =
sin(yy);
5673 const double sinyyneg =
sin(yyneg);
5674 r2 = (sinhx * sinhx + sinyy * sinyy) /
5675 (sinhx * sinhx + sinyyneg * sinyyneg);
5680 const double xxmirr = ty * (xpos + m_w[i].x - 2 * m_coplax);
5681 const std::complex<double> zzmirr(xxmirr, yy);
5682 const std::complex<double> zznmirr(xxmirr, yyneg);
5683 if (
fabs(xxmirr) <= 20.) {
5685 icons * m_b2sin[i] / (
sin(icons * zzmirr) *
sin(icons * zznmirr));
5687 const double sinhx = sinh(xxmirr);
5688 const double sinyy =
sin(yy);
5689 const double sinyyneg =
sin(yyneg);
5690 const double r2plan = (sinhx * sinhx + sinyy * sinyy) /
5691 (sinhx * sinhx + sinyyneg * sinyyneg);
5697 ex += real(m_sigmat[isw][i]) * real(ecompl);
5698 ey -= real(m_sigmat[isw][i]) * imag(ecompl);
5699 if (opt) volt -= 0.5 * real(m_sigmat[isw][i]) * log(r2);
5705void ComponentAnalyticField::WfieldWireC2X(
const double xpos,
const double ypos,
5706 double& ex,
double& ey,
double& volt,
5707 const int isw,
const bool opt)
const {
5716 const std::complex<double> icons = std::complex<double>(0., 1.);
5717 std::complex<double> zsin, zcof, zu, zunew, zterm1, zterm2, zeta;
5720 std::complex<double> wsum1 = 0.;
5721 std::complex<double> wsum2 = 0.;
5726 for (
unsigned int i = 0; i < m_nWires; ++i) {
5728 zeta = m_zmult * std::complex<double>(xpos - m_w[i].x, ypos - m_w[i].y);
5729 if (imag(zeta) > 15.) {
5730 wsum1 -= real(m_sigmat[isw][i]) * icons;
5731 if (opt) volt -= real(m_sigmat[isw][i]) * (
fabs(imag(zeta)) - CLog2);
5732 }
else if (imag(zeta) < -15.) {
5733 wsum1 += real(m_sigmat[isw][i]) * icons;
5734 if (opt) volt -= real(m_sigmat[isw][i]) * (
fabs(imag(zeta)) - CLog2);
5737 zcof = 4. * zsin * zsin - 2.;
5738 zu = -m_p1 - zcof * m_p2;
5739 zunew = 1. - zcof * zu - m_p2;
5740 zterm1 = (zunew + zu) * zsin;
5741 zu = -3. * m_p1 - zcof * 5. * m_p2;
5742 zunew = 1. - zcof * zu - 5. * m_p2;
5743 zterm2 = (zunew - zu) *
cos(zeta);
5744 wsum1 += real(m_sigmat[isw][i]) * (zterm2 / zterm1);
5745 if (opt) volt -= real(m_sigmat[isw][i]) * log(abs(zterm1));
5748 double cx = m_coplax - m_sx * int(round((m_coplax - m_w[i].x) / m_sx));
5750 s += real(m_sigmat[isw][i]) * (m_w[i].x - cx);
5753 std::complex<double>(2. * cx - xpos - m_w[i].x, ypos - m_w[i].y);
5754 if (imag(zeta) > +15.) {
5755 wsum2 -= real(m_sigmat[isw][i]) * icons;
5756 if (opt) volt += real(m_sigmat[isw][i]) * (
fabs(imag(zeta)) - CLog2);
5757 }
else if (imag(zeta) < -15.) {
5758 wsum2 += real(m_sigmat[isw][i]) * icons;
5759 if (opt) volt += real(m_sigmat[isw][i]) * (
fabs(imag(zeta)) - CLog2);
5762 zcof = 4. * zsin * zsin - 2.;
5763 zu = -m_p1 - zcof * m_p2;
5764 zunew = 1. - zcof * zu - m_p2;
5765 zterm1 = (zunew + zu) * zsin;
5766 zu = -3. * m_p1 - zcof * 5. * m_p2;
5767 zunew = 1. - zcof * zu - 5. * m_p2;
5768 zterm2 = (zunew - zu) *
cos(zeta);
5769 wsum2 += real(m_sigmat[isw][i]) * (zterm2 / zterm1);
5770 if (opt) volt += real(m_sigmat[isw][i]) * log(abs(zterm1));
5773 if (opt && m_mode == 0) {
5774 volt -= TwoPi * real(m_sigmat[isw][i]) * (xpos - cx) * (m_w[i].x - cx) /
5779 ex = real(m_zmult * (wsum1 + wsum2));
5780 ey = -imag(m_zmult * (wsum1 - wsum2));
5782 if (m_mode == 0) ex += s * TwoPi / (m_sx * m_sy);
5785void ComponentAnalyticField::WfieldWireC2Y(
const double xpos,
const double ypos,
5786 double& ex,
double& ey,
double& volt,
5787 const int isw,
const bool opt)
const {
5796 const std::complex<double> icons = std::complex<double>(0., 1.);
5797 std::complex<double> zsin, zcof, zu, zunew, zterm1, zterm2, zeta;
5800 std::complex<double> wsum1 = 0.;
5801 std::complex<double> wsum2 = 0.;
5806 for (
unsigned int i = 0; i < m_nWires; ++i) {
5808 zeta = m_zmult * std::complex<double>(xpos - m_w[i].x, ypos - m_w[i].y);
5809 if (imag(zeta) > +15.) {
5810 wsum1 -= real(m_sigmat[isw][i]) * icons;
5811 if (opt) volt -= real(m_sigmat[isw][i]) * (
fabs(imag(zeta)) - CLog2);
5812 }
else if (imag(zeta) < -15.) {
5813 wsum1 += real(m_sigmat[isw][i]) * icons;
5814 if (opt) volt -= real(m_sigmat[isw][i]) * (
fabs(imag(zeta)) - CLog2);
5817 zcof = 4. * zsin * zsin - 2.;
5818 zu = -m_p1 - zcof * m_p2;
5819 zunew = 1. - zcof * zu - m_p2;
5820 zterm1 = (zunew + zu) * zsin;
5821 zu = -3. * m_p1 - zcof * 5. * m_p2;
5822 zunew = 1. - zcof * zu - 5. * m_p2;
5823 zterm2 = (zunew - zu) *
cos(zeta);
5824 wsum1 += real(m_sigmat[isw][i]) * (zterm2 / zterm1);
5825 if (opt) volt -= real(m_sigmat[isw][i]) * log(abs(zterm1));
5828 double cy = m_coplay - m_sy * int(round((m_coplay - m_w[i].y) / m_sy));
5830 s += real(m_sigmat[isw][i]) * (m_w[i].y - cy);
5833 std::complex<double>(xpos - m_w[i].x, 2. * cy - ypos - m_w[i].y);
5834 if (imag(zeta) > +15.) {
5835 wsum2 -= real(m_sigmat[isw][i]) * icons;
5836 if (opt) volt += real(m_sigmat[isw][i]) * (
fabs(imag(zeta)) - CLog2);
5837 }
else if (imag(zeta) < -15.) {
5838 wsum2 += real(m_sigmat[isw][i]) * icons;
5839 if (opt) volt += real(m_sigmat[isw][i]) * (
fabs(imag(zeta)) - CLog2);
5842 zcof = 4. * zsin * zsin - 2.;
5843 zu = -m_p1 - zcof * m_p2;
5844 zunew = 1. - zcof * zu - m_p2;
5845 zterm1 = (zunew + zu) * zsin;
5846 zu = -3. * m_p1 - zcof * 5. * m_p2;
5847 zunew = 1. - zcof * zu - 5. * m_p2;
5848 zterm2 = (zunew - zu) *
cos(zeta);
5849 wsum2 += real(m_sigmat[isw][i]) * (zterm2 / zterm1);
5850 if (opt) volt += real(m_sigmat[isw][i]) * log(abs(zterm1));
5853 if (opt && m_mode == 1) {
5854 volt -= TwoPi * real(m_sigmat[isw][i]) * (ypos - cy) * (m_w[i].y - cy) /
5859 ex = real(m_zmult * (wsum1 - wsum2));
5860 ey = -imag(m_zmult * (wsum1 + wsum2));
5862 if (m_mode == 1) ey += s * TwoPi / (m_sx * m_sy);
5865void ComponentAnalyticField::WfieldWireC30(
const double xpos,
const double ypos,
5866 double& ex,
double& ey,
double& volt,
5867 const int isw,
const bool opt)
const {
5876 const std::complex<double> icons = std::complex<double>(0., 1.);
5877 std::complex<double> zsin, zcof, zu, zunew, zterm1, zterm2, zeta;
5880 std::complex<double> wsum1 = 0.;
5881 std::complex<double> wsum2 = 0.;
5882 std::complex<double> wsum3 = 0.;
5883 std::complex<double> wsum4 = 0.;
5887 for (
unsigned int i = 0; i < m_nWires; ++i) {
5889 zeta = m_zmult * std::complex<double>(xpos - m_w[i].x, ypos - m_w[i].y);
5890 if (imag(zeta) > +15.) {
5891 wsum1 -= real(m_sigmat[isw][i]) * icons;
5892 if (opt) volt -= real(m_sigmat[isw][i]) * (
fabs(imag(zeta)) - CLog2);
5893 }
else if (imag(zeta) < -15.) {
5894 wsum1 += real(m_sigmat[isw][i]) * icons;
5895 if (opt) volt -= real(m_sigmat[isw][i]) * (
fabs(imag(zeta)) - CLog2);
5898 zcof = 4. * zsin * zsin - 2.;
5899 zu = -m_p1 - zcof * m_p2;
5900 zunew = 1. - zcof * zu - m_p2;
5901 zterm1 = (zunew + zu) * zsin;
5902 zu = -3. * m_p1 - zcof * 5. * m_p2;
5903 zunew = 1. - zcof * zu - 5. * m_p2;
5904 zterm2 = (zunew - zu) *
cos(zeta);
5905 wsum1 += real(m_sigmat[isw][i]) * (zterm2 / zterm1);
5906 if (opt) volt -= real(m_sigmat[isw][i]) * log(abs(zterm1));
5909 double cx = m_coplax - m_sx * int(round((m_coplax - m_w[i].x) / m_sx));
5912 std::complex<double>(2. * cx - xpos - m_w[i].x, ypos - m_w[i].y);
5913 if (imag(zeta) > +15.) {
5914 wsum2 -= real(m_sigmat[isw][i]) * icons;
5915 if (opt) volt += real(m_sigmat[isw][i]) * (
fabs(imag(zeta)) - CLog2);
5916 }
else if (imag(zeta) < -15.) {
5917 wsum2 += real(m_sigmat[isw][i]) * icons;
5918 if (opt) volt += real(m_sigmat[isw][i]) * (
fabs(imag(zeta)) - CLog2);
5921 zcof = 4. * zsin * zsin - 2.;
5922 zu = -m_p1 - zcof * m_p2;
5923 zunew = 1. - zcof * zu - m_p2;
5924 zterm1 = (zunew + zu) * zsin;
5925 zu = -3. * m_p1 - zcof * 5. * m_p2;
5926 zunew = 1. - zcof * zu - 5. * m_p2;
5927 zterm2 = (zunew - zu) *
cos(zeta);
5928 wsum2 += real(m_sigmat[isw][i]) * (zterm2 / zterm1);
5929 if (opt) volt += real(m_sigmat[isw][i]) * log(abs(zterm1));
5932 double cy = m_coplay - m_sy * int(round((m_coplay - m_w[i].y) / m_sy));
5935 std::complex<double>(xpos - m_w[i].x, 2. * cy - ypos - m_w[i].y);
5936 if (imag(zeta) > +15.) {
5937 wsum3 -= real(m_sigmat[isw][i]) * icons;
5938 if (opt) volt += real(m_sigmat[isw][i]) * (
fabs(imag(zeta)) - CLog2);
5939 }
else if (imag(zeta) < -15.) {
5940 wsum3 += real(m_sigmat[isw][i]) * icons;
5941 if (opt) volt += real(m_sigmat[isw][i]) * (
fabs(imag(zeta)) - CLog2);
5944 zcof = 4. * zsin * zsin - 2.;
5945 zu = -m_p1 - zcof * m_p2;
5946 zunew = 1. - zcof * zu - m_p2;
5947 zterm1 = (zunew + zu) * zsin;
5948 zu = -3. * m_p1 - zcof * 5. * m_p2;
5949 zunew = 1. - zcof * zu - 5. * m_p2;
5950 zterm2 = (zunew - zu) *
cos(zeta);
5951 wsum3 += real(m_sigmat[isw][i]) * (zterm2 / zterm1);
5952 if (opt) volt += real(m_sigmat[isw][i]) * log(abs(zterm1));
5955 zeta = m_zmult * std::complex<double>(2. * cx - xpos - m_w[i].x,
5956 2. * cy - ypos - m_w[i].y);
5957 if (imag(zeta) > +15.) {
5958 wsum4 -= real(m_sigmat[isw][i]) * icons;
5959 if (opt) volt -= real(m_sigmat[isw][i]) * (
fabs(imag(zeta)) - CLog2);
5960 }
else if (imag(zeta) < -15.) {
5961 wsum4 += real(m_sigmat[isw][i]) * icons;
5962 if (opt) volt -= real(m_sigmat[isw][i]) * (
fabs(imag(zeta)) - CLog2);
5965 zcof = 4. * zsin * zsin - 2.;
5966 zu = -m_p1 - zcof * m_p2;
5967 zunew = 1. - zcof * zu - m_p2;
5968 zterm1 = (zunew + zu) * zsin;
5969 zu = -3. * m_p1 - zcof * 5. * m_p2;
5970 zunew = 1. - zcof * zu - 5. * m_p2;
5971 zterm2 = (zunew - zu) *
cos(zeta);
5972 wsum4 += real(m_sigmat[isw][i]) * (zterm2 / zterm1);
5973 if (opt) volt -= real(m_sigmat[isw][i]) * log(abs(zterm1));
5977 ex = real(m_zmult * (wsum1 + wsum2 - wsum3 - wsum4));
5978 ey = -imag(m_zmult * (wsum1 - wsum2 + wsum3 - wsum4));
5981void ComponentAnalyticField::WfieldWireD10(
const double xpos,
const double ypos,
5982 double& ex,
double& ey,
double& volt,
5983 const int isw,
const bool opt)
const {
5996 ex = ey = volt = 0.;
5999 std::complex<double> zpos = std::complex<double>(xpos, ypos);
6000 std::complex<double> zi;
6001 std::complex<double> wi;
6003 for (
int i = m_nWires; i--;) {
6005 zi = std::complex<double>(m_w[i].x, m_w[i].y);
6008 volt -= real(m_sigmat[isw][i]) *
6009 log(abs(m_cotube * (zpos - zi) / (m_cotube2 - zpos * conj(zi))));
6012 wi = 1. / conj(zpos - zi) + zi / (m_cotube2 - conj(zpos) * zi);
6013 ex += real(m_sigmat[isw][i]) * real(wi);
6014 ey += real(m_sigmat[isw][i]) * imag(wi);
6018void ComponentAnalyticField::WfieldWireD30(
const double xpos,
const double ypos,
6019 double& ex,
double& ey,
double& volt,
6020 const int isw,
const bool opt)
const {
6032 ex = ey = volt = 0.;
6034 std::complex<double> whelp;
6037 std::complex<double> wpos, wdpos;
6038 ConformalMap(std::complex<double>(xpos, ypos) / m_cotube, wpos, wdpos);
6040 for (
int i = m_nWires; i--;) {
6043 volt -= real(m_sigmat[isw][i]) *
6044 log(abs((wpos - wmap[i]) / (1. - wpos * conj(wmap[i]))));
6047 whelp = wdpos * (1. -
pow(abs(wmap[i]), 2)) /
6048 ((wpos - wmap[i]) * (1. - conj(wmap[i]) * wpos));
6049 ex += real(m_sigmat[isw][i]) * real(whelp);
6050 ey -= real(m_sigmat[isw][i]) * imag(whelp);
6056void ComponentAnalyticField::WfieldPlaneA00(
const double xpos,
6057 const double ypos,
double& ex,
6058 double& ey,
double& volt,
6059 const int mx,
const int my,
6061 const bool opt)
const {
6073 ex = ey = volt = 0.;
6075 double xxmirr = 0., yymirr = 0.;
6077 for (
int i = m_nWires; i--;) {
6079 const double xx = xpos - m_w[i].x - mx * m_sx;
6080 const double yy = ypos - m_w[i].y - my * m_sy;
6082 double r2 = xx * xx + yy * yy;
6083 if (r2 <= 0.)
continue;
6084 double exhelp = xx / r2;
6085 double eyhelp = yy / r2;
6088 xxmirr = xpos + m_w[i].x - 2 * m_coplax;
6089 const double r2plan = xxmirr * xxmirr + yy * yy;
6090 if (r2plan <= 0.)
continue;
6091 exhelp -= xxmirr / r2plan;
6092 eyhelp -= yy / r2plan;
6097 yymirr = ypos + m_w[i].y - 2 * m_coplay;
6098 const double r2plan = xx * xx + yymirr * yymirr;
6099 if (r2plan <= 0.)
continue;
6100 exhelp -= xx / r2plan;
6101 eyhelp -= yymirr / r2plan;
6105 if (m_ynplax && m_ynplay) {
6106 const double r2plan = xxmirr * xxmirr + yymirr * yymirr;
6107 if (r2plan <= 0.)
continue;
6108 exhelp += xxmirr / r2plan;
6109 eyhelp += yymirr / r2plan;
6113 if (opt) volt -= 0.5 * m_qplane[iplane][i] * log(r2);
6114 ex += m_qplane[iplane][i] * exhelp;
6115 ey += m_qplane[iplane][i] * eyhelp;
6119void ComponentAnalyticField::WfieldPlaneB2X(
const double xpos,
6120 const double ypos,
double& ex,
6121 double& ey,
double& volt,
6122 const int my,
const int iplane,
6123 const bool opt)
const {
6135 ex = ey = volt = 0.;
6137 const double tx = HalfPi / m_sx;
6139 for (
unsigned int i = 0; i < m_nWires; ++i) {
6140 const double xx = tx * (xpos - m_w[i].x);
6141 const double yy = tx * (ypos - m_w[i].y - my * m_sy);
6142 const double xxneg = tx * (xpos + m_w[i].x - 2 * m_coplan[0]);
6143 const std::complex<double> zz(xx, yy);
6144 const std::complex<double> zzneg(xxneg, yy);
6146 std::complex<double> ecompl(0., 0.);
6148 if (
fabs(yy) <= 20.) {
6149 ecompl = -m_b2sin[i] / (
sin(zz) *
sin(zzneg));
6151 const double sinhy = sinh(yy);
6152 const double sinxx =
sin(xx);
6153 const double sinxxneg =
sin(xxneg);
6154 r2 = (sinhy * sinhy + sinxx * sinxx) /
6155 (sinhy * sinhy + sinxxneg * sinxxneg);
6160 const double yymirr = tx * (ypos + m_w[i].y - 2 * m_coplay);
6161 const std::complex<double> zzmirr(yy, yymirr);
6162 const std::complex<double> zznmirr(xxneg, yymirr);
6163 if (
fabs(yymirr) <= 20.) {
6164 ecompl += m_b2sin[i] / (
sin(zzmirr) *
sin(zznmirr));
6166 const double sinhy = sinh(yymirr);
6167 const double sinxx =
sin(xx);
6168 const double sinxxneg =
sin(xxneg);
6169 const double r2plan = (sinhy * sinhy + sinxx * sinxx) /
6170 (sinhy * sinhy + sinxxneg * sinxxneg);
6176 ex += m_qplane[iplane][i] * real(ecompl);
6177 ey -= m_qplane[iplane][i] * imag(ecompl);
6178 if (opt) volt -= 0.5 * m_qplane[iplane][i] * log(r2);
6184void ComponentAnalyticField::WfieldPlaneB2Y(
const double xpos,
6185 const double ypos,
double& ex,
6186 double& ey,
double& volt,
6187 const int mx,
const int iplane,
6188 const bool opt)
const {
6199 const std::complex<double> icons = std::complex<double>(0., 1.);
6202 ex = ey = volt = 0.;
6204 const double ty = HalfPi / m_sy;
6206 for (
unsigned int i = 0; i < m_nWires; ++i) {
6207 const double xx = ty * (xpos - m_w[i].x - mx * m_sx);
6208 const double yy = ty * (ypos - m_w[i].y);
6209 const double yyneg = ty * (ypos + m_w[i].y - 2 * m_coplan[2]);
6210 const std::complex<double> zz(xx, yy);
6211 const std::complex<double> zzneg(xx, yyneg);
6213 std::complex<double> ecompl(0., 0.);
6215 if (
fabs(xx) <= 20.) {
6216 ecompl = icons * m_b2sin[i] / (
sin(icons * zz) *
sin(icons * zzneg));
6218 const double sinhx = sinh(xx);
6219 const double sinyy =
sin(yy);
6220 const double sinyyneg =
sin(yyneg);
6221 r2 = (sinhx * sinhx + sinyy * sinyy) /
6222 (sinhx * sinhx + sinyyneg * sinyyneg);
6227 const double xxmirr = ty * (xpos + m_w[i].x - 2 * m_coplax);
6228 const std::complex<double> zzmirr(xxmirr, yy);
6229 const std::complex<double> zznmirr(xxmirr, yyneg);
6230 if (
fabs(xxmirr) <= 20.) {
6231 ecompl -= m_b2sin[i] / (
sin(icons * zzmirr) *
sin(icons * zznmirr));
6233 const double sinhx = sinh(xxmirr);
6234 const double sinyy =
sin(yy);
6235 const double sinyyneg =
sin(yyneg);
6236 const double r2plan = (sinhx * sinhx + sinyy * sinyy) /
6237 (sinhx * sinhx + sinyyneg * sinyyneg);
6243 ex += m_qplane[iplane][i] * real(ecompl);
6244 ey -= m_qplane[iplane][i] * imag(ecompl);
6245 if (opt) volt -= 0.5 * m_qplane[iplane][i] * log(r2);
6251void ComponentAnalyticField::WfieldPlaneC2X(
const double xpos,
6252 const double ypos,
double& ex,
6253 double& ey,
double& volt,
6255 const bool opt)
const {
6264 const std::complex<double> icons = std::complex<double>(0., 1.);
6266 std::complex<double> zsin, zcof, zu, zunew, zterm1, zterm2, zeta;
6268 std::complex<double> wsum1 = 0.;
6269 std::complex<double> wsum2 = 0.;
6274 for (
unsigned int i = 0; i < m_nWires; ++i) {
6276 zeta = m_zmult * std::complex<double>(xpos - m_w[i].x, ypos - m_w[i].y);
6277 if (imag(zeta) > +15.) {
6278 wsum1 -= m_qplane[iplane][i] * icons;
6279 if (opt) volt -= m_qplane[iplane][i] * (
fabs(imag(zeta)) - CLog2);
6280 }
else if (imag(zeta) < -15.) {
6281 wsum1 += m_qplane[iplane][i] * icons;
6282 if (opt) volt -= m_qplane[iplane][i] * (
fabs(imag(zeta)) - CLog2);
6285 zcof = 4. * zsin * zsin - 2.;
6286 zu = -m_p1 - zcof * m_p2;
6287 zunew = 1. - zcof * zu - m_p2;
6288 zterm1 = (zunew + zu) * zsin;
6289 zu = -3. * m_p1 - zcof * 5. * m_p2;
6290 zunew = 1. - zcof * zu - 5. * m_p2;
6291 zterm2 = (zunew - zu) *
cos(zeta);
6292 wsum1 += m_qplane[iplane][i] * (zterm2 / zterm1);
6293 if (opt) volt -= m_qplane[iplane][i] * log(abs(zterm1));
6296 double cx = m_coplax - m_sx * int(round((m_coplax - m_w[i].x) / m_sx));
6298 s += m_qplane[iplane][i] * (m_w[i].x - cx);
6301 std::complex<double>(2. * cx - xpos - m_w[i].x, ypos - m_w[i].y);
6302 if (imag(zeta) > 15.) {
6303 wsum2 -= m_qplane[iplane][i] * icons;
6304 if (opt) volt += m_qplane[iplane][i] * (
fabs(imag(zeta)) - CLog2);
6305 }
else if (imag(zeta) < -15.) {
6306 wsum2 += m_qplane[iplane][i] * icons;
6307 if (opt) volt += m_qplane[iplane][i] * (
fabs(imag(zeta)) - CLog2);
6310 zcof = 4. * zsin * zsin - 2.;
6311 zu = -m_p1 - zcof * m_p2;
6312 zunew = 1. - zcof * zu - m_p2;
6313 zterm1 = (zunew + zu) * zsin;
6314 zu = -3. * m_p1 - zcof * 5. * m_p2;
6315 zunew = 1. - zcof * zu - 5. * m_p2;
6316 zterm2 = (zunew - zu) *
cos(zeta);
6317 wsum2 += m_qplane[iplane][i] * (zterm2 / zterm1);
6318 if (opt) volt += m_qplane[iplane][i] * log(abs(zterm1));
6320 if (opt && m_mode == 0) {
6321 volt -= TwoPi * m_qplane[iplane][i] * (xpos - cx) * (m_w[i].x - cx) /
6326 ex = real(m_zmult * (wsum1 + wsum2));
6327 ey = -imag(m_zmult * (wsum1 - wsum2));
6329 if (m_mode == 0) ex += s * TwoPi / (m_sx * m_sy);
6332void ComponentAnalyticField::WfieldPlaneC2Y(
const double xpos,
6333 const double ypos,
double& ex,
6334 double& ey,
double& volt,
6336 const bool opt)
const {
6345 const std::complex<double> icons = std::complex<double>(0., 1.);
6347 std::complex<double> zsin, zcof, zu, zunew, zterm1, zterm2, zeta;
6349 std::complex<double> wsum1 = 0.;
6350 std::complex<double> wsum2 = 0.;
6355 for (
unsigned int i = 0; i < m_nWires; ++i) {
6357 zeta = m_zmult * std::complex<double>(xpos - m_w[i].x, ypos - m_w[i].y);
6358 if (imag(zeta) > +15.) {
6359 wsum1 -= m_qplane[iplane][i] * icons;
6360 if (opt) volt -= m_qplane[iplane][i] * (
fabs(imag(zeta)) - CLog2);
6361 }
else if (imag(zeta) < -15.) {
6362 wsum1 += m_qplane[iplane][i] * icons;
6363 if (opt) volt -= m_qplane[iplane][i] * (
fabs(imag(zeta)) - CLog2);
6366 zcof = 4. * zsin * zsin - 2.;
6367 zu = -m_p1 - zcof * m_p2;
6368 zunew = 1. - zcof * zu - m_p2;
6369 zterm1 = (zunew + zu) * zsin;
6370 zu = -3. * m_p1 - zcof * 5. * m_p2;
6371 zunew = 1. - zcof * zu - 5. * m_p2;
6372 zterm2 = (zunew - zu) *
cos(zeta);
6373 wsum1 += m_qplane[iplane][i] * (zterm2 / zterm1);
6374 if (opt) volt -= m_qplane[iplane][i] * log(abs(zterm1));
6377 double cy = m_coplay - m_sy * int(round((m_coplay - m_w[i].y) / m_sy));
6379 s += m_qplane[iplane][i] * (m_w[i].y - cy);
6382 std::complex<double>(xpos - m_w[i].x, 2. * cy - ypos - m_w[i].y);
6383 if (imag(zeta) > 15.) {
6384 wsum2 -= m_qplane[iplane][i] * icons;
6385 if (opt) volt += m_qplane[iplane][i] * (
fabs(imag(zeta)) - CLog2);
6386 }
else if (imag(zeta) < -15.) {
6387 wsum2 += m_qplane[iplane][i] * icons;
6388 if (opt) volt += m_qplane[iplane][i] * (
fabs(imag(zeta)) - CLog2);
6391 zcof = 4. * zsin * zsin - 2.;
6392 zu = -m_p1 - zcof * m_p2;
6393 zunew = 1. - zcof * zu - m_p2;
6394 zterm1 = (zunew + zu) * zsin;
6395 zu = -3. * m_p1 - zcof * 5. * m_p2;
6396 zunew = 1. - zcof * zu - 5. * m_p2;
6397 zterm2 = (zunew - zu) *
cos(zeta);
6398 wsum2 += m_qplane[iplane][i] * (zterm2 / zterm1);
6399 if (opt) volt += m_qplane[iplane][i] * log(abs(zterm1));
6402 if (opt && m_mode == 1) {
6403 volt -= TwoPi * m_qplane[iplane][i] * (ypos - cy) * (m_w[i].y - cy) /
6408 ex = real(m_zmult * (wsum1 - wsum2));
6409 ey = -imag(m_zmult * (wsum1 + wsum2));
6411 if (m_mode == 1) ey += s * TwoPi / (m_sx * m_sy);
6414void ComponentAnalyticField::WfieldPlaneC30(
const double xpos,
6415 const double ypos,
double& ex,
6416 double& ey,
double& volt,
6418 const bool opt)
const {
6427 const std::complex<double> icons = std::complex<double>(0., 1.);
6429 std::complex<double> zsin, zcof, zu, zunew, zterm1, zterm2, zeta;
6431 std::complex<double> wsum1 = 0.;
6432 std::complex<double> wsum2 = 0.;
6433 std::complex<double> wsum3 = 0.;
6434 std::complex<double> wsum4 = 0.;
6438 for (
unsigned int i = 0; i < m_nWires; ++i) {
6440 zeta = m_zmult * std::complex<double>(xpos - m_w[i].x, ypos - m_w[i].y);
6441 if (imag(zeta) > +15.) {
6442 wsum1 -= m_qplane[iplane][i] * icons;
6443 if (opt) volt -= m_qplane[iplane][i] * (
fabs(imag(zeta)) - CLog2);
6444 }
else if (imag(zeta) < -15.) {
6445 wsum1 += m_qplane[iplane][i] * icons;
6446 if (opt) volt -= m_qplane[iplane][i] * (
fabs(imag(zeta)) - CLog2);
6449 zcof = 4. * zsin * zsin - 2.;
6450 zu = -m_p1 - zcof * m_p2;
6451 zunew = 1. - zcof * zu - m_p2;
6452 zterm1 = (zunew + zu) * zsin;
6453 zu = -3. * m_p1 - zcof * 5. * m_p2;
6454 zunew = 1. - zcof * zu - 5. * m_p2;
6455 zterm2 = (zunew - zu) *
cos(zeta);
6456 wsum1 += m_qplane[iplane][i] * zterm2 / zterm1;
6457 if (opt) volt -= m_qplane[iplane][i] * log(abs(zterm1));
6460 double cx = m_coplax - m_sx * int(round((m_coplax - m_w[i].x) / m_sx));
6463 std::complex<double>(2. * cx - xpos - m_w[i].x, ypos - m_w[i].y);
6464 if (imag(zeta) > 15.) {
6465 wsum2 -= m_qplane[iplane][i] * icons;
6466 if (opt) volt += m_qplane[iplane][i] * (
fabs(imag(zeta)) - CLog2);
6467 }
else if (imag(zeta) < -15.) {
6468 wsum2 += m_qplane[iplane][i] * icons;
6469 if (opt) volt += m_qplane[iplane][i] * (
fabs(imag(zeta)) - CLog2);
6472 zcof = 4. * zsin * zsin - 2.;
6473 zu = -m_p1 - zcof * m_p2;
6474 zunew = 1. - zcof * zu - m_p2;
6475 zterm1 = (zunew + zu) * zsin;
6476 zu = -3. * m_p1 - zcof * 5. * m_p2;
6477 zunew = 1. - zcof * zu - 5. * m_p2;
6478 zterm2 = (zunew - zu) *
cos(zeta);
6479 wsum2 += m_qplane[iplane][i] * zterm2 / zterm1;
6480 if (opt) volt += m_qplane[iplane][i] * log(abs(zterm1));
6483 double cy = m_coplay - m_sy * int(round((m_coplay - m_w[i].y) / m_sy));
6486 std::complex<double>(xpos - m_w[i].x, 2. * cy - ypos - m_w[i].y);
6487 if (imag(zeta) > 15.) {
6488 wsum3 -= m_qplane[iplane][i] * icons;
6489 if (opt) volt += m_qplane[iplane][i] * (
fabs(imag(zeta)) - CLog2);
6490 }
else if (imag(zeta) < -15.) {
6491 wsum3 += m_qplane[iplane][i] * icons;
6492 if (opt) volt += m_qplane[iplane][i] * (
fabs(imag(zeta)) - CLog2);
6495 zcof = 4. * zsin * zsin - 2.;
6496 zu = -m_p1 - zcof * m_p2;
6497 zunew = 1. - zcof * zu - m_p2;
6498 zterm1 = (zunew + zu) * zsin;
6499 zu = -3. * m_p1 - zcof * 5. * m_p2;
6500 zunew = 1. - zcof * zu - 5. * m_p2;
6501 zterm2 = (zunew - zu) *
cos(zeta);
6502 wsum3 += m_qplane[iplane][i] * zterm2 / zterm1;
6503 if (opt) volt += m_qplane[iplane][i] * log(abs(zterm1));
6506 zeta = m_zmult * std::complex<double>(2. * cx - xpos - m_w[i].x,
6507 2. * cy - ypos - m_w[i].y);
6508 if (imag(zeta) > 15.) {
6509 wsum4 -= m_qplane[iplane][i] * icons;
6510 if (opt) volt -= m_qplane[iplane][i] * (
fabs(imag(zeta)) - CLog2);
6511 }
else if (imag(zeta) < -15.) {
6512 wsum4 += m_qplane[iplane][i] * icons;
6513 if (opt) volt -= m_qplane[iplane][i] * (
fabs(imag(zeta)) - CLog2);
6516 zcof = 4. * zsin * zsin - 2.;
6517 zu = -m_p1 - zcof * m_p2;
6518 zunew = 1. - zcof * zu - m_p2;
6519 zterm1 = (zunew + zu) * zsin;
6520 zu = -3. * m_p1 - zcof * 5. * m_p2;
6521 zunew = 1. - zcof * zu - 5. * m_p2;
6522 zterm2 = (zunew - zu) *
cos(zeta);
6523 wsum4 += m_qplane[iplane][i] * zterm2 / zterm1;
6524 if (opt) volt -= m_qplane[iplane][i] * log(abs(zterm1));
6527 ex = real(m_zmult * (wsum1 + wsum2 - wsum3 - wsum4));
6528 ey = -imag(m_zmult * (wsum1 - wsum2 + wsum3 - wsum4));
6531void ComponentAnalyticField::WfieldPlaneD10(
const double xpos,
6532 const double ypos,
double& ex,
6533 double& ey,
double& volt,
6535 const bool opt)
const {
6547 ex = ey = volt = 0.;
6550 std::complex<double> zpos = std::complex<double>(xpos, ypos);
6551 std::complex<double> zi;
6552 std::complex<double> wi;
6554 for (
int i = m_nWires; i--;) {
6556 zi = std::complex<double>(m_w[i].x, m_w[i].y);
6559 volt -= m_qplane[iplane][i] *
6560 log(abs(m_cotube * (zpos - zi) / (m_cotube2 - zpos * conj(zi))));
6563 wi = 1. / conj(zpos - zi) + zi / (m_cotube2 - conj(zpos) * zi);
6564 ex += m_qplane[iplane][i] * real(wi);
6565 ey += m_qplane[iplane][i] * imag(wi);
6569void ComponentAnalyticField::WfieldPlaneD30(
const double xpos,
6570 const double ypos,
double& ex,
6571 double& ey,
double& volt,
6573 const bool opt)
const {
6585 ex = ey = volt = 0.;
6587 std::complex<double> whelp;
6590 std::complex<double> wpos, wdpos;
6591 ConformalMap(std::complex<double>(xpos, ypos) / m_cotube, wpos, wdpos);
6593 for (
unsigned int i = 0; i < m_nWires; ++i) {
6596 volt -= m_qplane[iplane][i] *
6597 log(abs((wpos - wmap[i]) / (1. - wpos * conj(wmap[i]))));
6600 whelp = wdpos * (1. -
pow(abs(wmap[i]), 2)) /
6601 ((wpos - wmap[i]) * (1. - conj(wmap[i]) * wpos));
6602 ex += m_qplane[iplane][i] * real(whelp);
6603 ey -= m_qplane[iplane][i] * imag(whelp);
6609void ComponentAnalyticField::WfieldStripZ(
const double xpos,
const double ypos,
6610 double& ex,
double& ey,
double& volt,
6611 const int ip,
const int is,
6612 const bool opt)
const {
6620 ex = ey = volt = 0.;
6622 strip theStrip = planes[ip].strips2[is];
6624 double xw = 0., yw = 0.;
6627 xw = -ypos + 0.5 * (theStrip.smin + theStrip.smax);
6628 yw = xpos - m_coplan[ip];
6631 xw = ypos - 0.5 * (theStrip.smin + theStrip.smax);
6632 yw = m_coplan[ip] - xpos;
6635 xw = xpos - 0.5 * (theStrip.smin + theStrip.smax);
6636 yw = ypos - m_coplan[ip];
6639 xw = -xpos + 0.5 * (theStrip.smin + theStrip.smax);
6640 yw = m_coplan[ip] - ypos;
6646 const double w = 0.5 *
fabs(theStrip.smax - theStrip.smin);
6647 const double g = theStrip.gap;
6650 if (yw <= 0. || yw > g)
return;
6653 const double s =
sin(Pi * yw / g);
6654 const double c =
cos(Pi * yw / g);
6655 const double e1 =
exp(Pi * (w - xw) / g);
6656 const double e2 =
exp(-Pi * (w + xw) / g);
6657 const double ce12 =
pow(c - e1, 2);
6658 const double ce22 =
pow(c - e2, 2);
6660 if (c == e1 || c == e2)
return;
6663 volt = atan((c - e2) / s) - atan((c - e1) / s);
6667 const double ewx = (s / g) * (e1 / (ce12 + s * s) - e2 / (ce22 + s * s));
6668 const double ewy = ((c / (c - e2) + s * s / ce22) / (1. + s * s / ce22) -
6669 (c / (c - e1) + s * s / ce12) / (1. + s * s / ce12)) /
6693void ComponentAnalyticField::WfieldStripXy(
const double xpos,
const double ypos,
6694 const double zpos,
double& ex,
6695 double& ey,
double& ez,
double& volt,
6696 const int ip,
const int is,
6697 const bool opt)
const {
6705 ex = ey = ez = volt = 0.;
6707 strip theStrip = planes[ip].strips1[is];
6709 double xw = 0., yw = 0.;
6712 xw = -zpos + 0.5 * (theStrip.smin + theStrip.smax);
6713 yw = xpos - m_coplan[ip];
6716 xw = zpos - 0.5 * (theStrip.smin + theStrip.smax);
6717 yw = m_coplan[ip] - xpos;
6720 xw = zpos - 0.5 * (theStrip.smin + theStrip.smax);
6721 yw = ypos - m_coplan[ip];
6724 xw = -zpos + 0.5 * (theStrip.smin + theStrip.smax);
6725 yw = m_coplan[ip] - ypos;
6732 const double w = 0.5 *
fabs(theStrip.smax - theStrip.smin);
6733 const double g = theStrip.gap;
6736 if (yw <= 0. || yw > g)
return;
6739 const double s =
sin(Pi * yw / g);
6740 const double c =
cos(Pi * yw / g);
6741 const double e1 =
exp(Pi * (w - xw) / g);
6742 const double e2 =
exp(-Pi * (w + xw) / g);
6743 const double ce12 =
pow(c - e1, 2);
6744 const double ce22 =
pow(c - e2, 2);
6746 if (c == e1 || c == e2)
return;
6749 volt = atan((c - e2) / s) - atan((c - e1) / s);
6753 const double ewx = (s / g) * (e1 / (ce12 + s * s) - e2 / (ce22 + s * s));
6754 const double ewy = ((c / (c - e2) + s * s / ce22) / (1. + s * s / ce22) -
6755 (c / (c - e1) + s * s / ce12) / (1. + s * s / ce12)) / g;
6782void ComponentAnalyticField::WfieldPixel(
const double xpos,
const double ypos,
6783 const double zpos,
double& ex,
6784 double& ey,
double& ez,
double& volt,
6785 const int ip,
const int is,
6786 const bool opt)
const {
6797 ex = ey = ez = volt = 0.;
6799 const pixel& thePixel = planes[ip].pixels[is];
6800 const double d = thePixel.gap;
6802 double x = 0., y = 0., z = 0.;
6805 const double ps = 0.5 * (thePixel.smin + thePixel.smax);
6806 const double pz = 0.5 * (thePixel.zmin + thePixel.zmax);
6807 const double wx = thePixel.smax - thePixel.smin;
6808 const double wy = thePixel.zmax - thePixel.zmin;
6813 z = xpos - m_coplan[ip];
6818 z = -xpos + m_coplan[ip];
6823 z = ypos - m_coplan[ip];
6828 z = -ypos + m_coplan[ip];
6833 if (z < 0.) std::cerr <<
" z = " << z << std::endl;
6840 const double x1 = x - 0.5 *
wx;
6841 const double x2 = x + 0.5 *
wx;
6842 const double y1 = y - 0.5 *
wy;
6843 const double y2 = y + 0.5 *
wy;
6844 const double x1s = x1 * x1;
6845 const double x2s = x2 * x2;
6846 const double y1s = y1 * y1;
6847 const double y2s = y2 * y2;
6850 const double maxError = 1.e-5;
6851 const double d3 = d * d * d;
6852 const unsigned int nz =
6853 std::ceil(
sqrt(wx * wy / (8 * Pi * d3 * maxError)));
6854 const unsigned int nx =
6855 std::ceil(
sqrt(wy * z / (4 * Pi * d3 * maxError)));
6856 const unsigned int ny =
6857 std::ceil(
sqrt(wx * z / (4 * Pi * d3 * maxError)));
6858 const unsigned int nn = std::max(ny, std::max(nx, nz));
6859 for (
unsigned int i = 1; i <= nn; ++i) {
6860 const double u1 = 2 * i * d - z;
6861 const double u2 = 2 * i * d + z;
6862 const double u1s = u1 * u1;
6863 const double u2s = u2 * u2;
6864 const double u1x1y1 =
sqrt(x1s + y1s + u1s);
6865 const double u1x1y2 =
sqrt(x1s + y2s + u1s);
6866 const double u1x2y1 =
sqrt(x2s + y1s + u1s);
6867 const double u1x2y2 =
sqrt(x2s + y2s + u1s);
6868 const double u2x1y1 =
sqrt(x1s + y1s + u2s);
6869 const double u2x1y2 =
sqrt(x1s + y2s + u2s);
6870 const double u2x2y1 =
sqrt(x2s + y1s + u2s);
6871 const double u2x2y2 =
sqrt(x2s + y2s + u2s);
6876 u1 * y1 / ((u1s + x2s) * u1x2y1) - u1 * y1 / ((u1s + x1s) * u1x1y1) +
6877 u1 * y2 / ((u1s + x1s) * u1x1y2) - u1 * y2 / ((u1s + x2s) * u1x2y2);
6881 u2 * y1 / ((u2s + x2s) * u2x2y1) - u2 * y1 / ((u2s + x1s) * u2x1y1) +
6882 u2 * y2 / ((u2s + x1s) * u2x1y2) - u2 * y2 / ((u2s + x2s) * u2x2y2);
6887 u1 * x1 / ((u1s + y2s) * u1x1y2) - u1 * x1 / ((u1s + y1s) * u1x1y1) +
6888 u1 * x2 / ((u1s + y1s) * u1x2y1) - u1 * x2 / ((u1s + y2s) * u1x2y2);
6892 u2 * x1 / ((u2s + y2s) * u2x1y2) - u2 * x1 / ((u2s + y1s) * u2x1y1) +
6893 u2 * x2 / ((u2s + y1s) * u2x2y1) - u2 * x2 / ((u2s + y2s) * u2x2y2);
6897 ez += x1 * y1 * (x1s + y1s + 2 * u1s) /
6898 ((x1s + u1s) * (y1s + u1s) * u1x1y1) +
6899 x2 * y2 * (x2s + y2s + 2 * u1s) /
6900 ((x2s + u1s) * (y2s + u1s) * u1x2y2) -
6901 x1 * y2 * (x1s + y2s + 2 * u1s) /
6902 ((x1s + u1s) * (y2s + u1s) * u1x1y2) -
6903 x2 * y1 * (x2s + y1s + 2 * u1s) /
6904 ((x2s + u1s) * (y1s + u1s) * u1x2y1);
6907 ez += x1 * y1 * (x1s + y1s + 2 * u2s) /
6908 ((x1s + u2s) * (y1s + u2s) * u2x1y1) +
6909 x2 * y2 * (x2s + y2s + 2 * u2s) /
6910 ((x2s + u2s) * (y2s + u2s) * u2x2y2) -
6911 x1 * y2 * (x1s + y2s + 2 * u2s) /
6912 ((x1s + u2s) * (y2s + u2s) * u2x1y2) -
6913 x2 * y1 * (x2s + y1s + 2 * u2s) /
6914 ((x2s + u2s) * (y1s + u2s) * u2x2y1);
6917 volt -= atan(x1 * y1 / (u1 * u1x1y1)) + atan(x2 * y2 / (u1 * u1x2y2)) -
6918 atan(x1 * y2 / (u1 * u1x1y2)) - atan(x2 * y1 / (u1 * u1x2y1));
6919 volt += atan(x1 * y1 / (u2 * u2x1y1)) + atan(x2 * y2 / (u2 * u2x2y2)) -
6920 atan(x1 * y2 / (u2 * u2x1y2)) - atan(x2 * y1 / (u2 * u2x2y1));
6923 const double zs = z * z;
6924 const double x1y1 =
sqrt(x1s + y1s + zs);
6925 const double x1y2 =
sqrt(x1s + y2s + zs);
6926 const double x2y1 =
sqrt(x2s + y1s + zs);
6927 const double x2y2 =
sqrt(x2s + y2s + zs);
6929 ex += z * y1 / ((zs + x2s) * x2y1) - z * y1 / ((zs + x1s) * x1y1) +
6930 z * y2 / ((zs + x1s) * x1y2) - z * y2 / ((zs + x2s) * x2y2);
6933 ey += z * x1 / ((zs + y2s) * x1y2) - z * x1 / ((zs + y1s) * x1y1) +
6934 z * x2 / ((zs + y1s) * x2y1) - z * x2 / ((zs + y2s) * x2y2);
6937 ez += x1 * y1 * (x1s + y1s + 2 * zs) / ((x1s + zs) * (y1s + zs) * x1y1) +
6938 x2 * y2 * (x2s + y2s + 2 * zs) / ((x2s + zs) * (y2s + zs) * x2y2) -
6939 x1 * y2 * (x1s + y2s + 2 * zs) / ((x1s + zs) * (y2s + zs) * x1y2) -
6940 x2 * y1 * (x2s + y1s + 2 * zs) / ((x2s + zs) * (y1s + zs) * x2y1);
6947 volt += atan(x1 * y1 / (z * x1y1)) + atan(x2 * y2 / (z * x2y2)) -
6948 atan(x1 * y2 / (z * x1y2)) - atan(x2 * y1 / (z * x2y1));
6953 const double fx = ex;
6954 const double fy = ey;
6955 const double fz = ez;
bool GetVoltageRange(double &pmin, double &pmax)
Calculate the voltage range [V].
void SetPeriodicityX(const double s)
Set the periodic length [cm] in the x-direction.
bool GetPeriodicityY(double &s)
void AddTube(const double radius, const double voltage, const int nEdges, const std::string &label)
Add a tube.
void AddPlaneX(const double x, const double voltage, const std::string &label)
Add a plane at constant x.
bool IsInTrapRadius(const double q0, const double x0, const double y0, const double z0, double &xw, double &yx, double &rw)
void AddPixelOnPlaneY(const double y, const double xmin, const double xmax, const double zmin, const double zmax, const std::string &label, const double gap=-1.)
void AddCharge(const double x, const double y, const double z, const double q)
Add a point charge.
bool GetWire(const unsigned int i, double &x, double &y, double &diameter, double &voltage, std::string &label, double &length, double &charge, int &ntrap) const
ComponentAnalyticField()
Constructor.
void AddStripOnPlaneX(const char direction, const double x, const double smin, const double smax, const std::string &label, const double gap=-1.)
unsigned int GetNumberOfPlanesY() const
void AddWire(const double x, const double y, const double diameter, const double voltage, const std::string &label, const double length=100., const double tension=50., const double rho=19.3, const int ntrap=5)
Add a wire at (x, y) .
unsigned int GetNumberOfPlanesX() const
bool GetPeriodicityX(double &s)
bool GetTube(double &r, double &voltage, int &nEdges, std::string &label) const
bool GetPlaneX(const unsigned int i, double &x, double &voltage, std::string &label) const
void AddPlaneY(const double y, const double voltage, const std::string &label)
Add a plane at constant y.
void AddPixelOnPlaneX(const double x, const double ymin, const double ymax, const double zmin, const double zmax, const std::string &label, const double gap=-1.)
bool GetBoundingBox(double &x0, double &y0, double &z0, double &x1, double &y1, double &z1)
Get the bounding box coordinates.
void AddStripOnPlaneY(const char direction, const double y, const double smin, const double smax, const std::string &label, const double gap=-1.)
void AddReadout(const std::string &label)
Setup the weighting field for a given group of wires or planes.
void SetPeriodicityY(const double s)
Set the periodic length [cm] in the y-direction.
bool IsWireCrossed(const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, double &xc, double &yc, double &zc)
void PrintCharges() const
bool GetPlaneY(const unsigned int i, double &y, double &voltage, std::string &label) const
bool m_zMirrorPeriodic
Mirror periodicity in z.
bool m_yAxiallyPeriodic
Axial periodicity in y.
GeometryBase * m_geometry
Pointer to the geometry.
bool m_zRotationSymmetry
Rotation symmetry around z-axis.
bool m_yRotationSymmetry
Rotation symmetry around y-axis.
bool m_zAxiallyPeriodic
Axial periodicity in z.
bool m_xRotationSymmetry
Rotation symmetry around x-axis.
bool m_yPeriodic
Simple periodicity in y.
bool m_yMirrorPeriodic
Mirror periodicity in y.
bool m_xPeriodic
Simple periodicity in x.
bool m_zPeriodic
Simple periodicity in z.
std::string m_className
Class name.
bool m_xAxiallyPeriodic
Axial periodicity in x.
bool m_debug
Switch on/off debugging messages.
bool m_xMirrorPeriodic
Mirror periodicity in x.
virtual bool GetBoundingBox(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)=0
double BesselK0L(const double xx)
double BesselK1S(const double xx)
void Cinv(const int n, std::vector< std::vector< std::complex< double > > > &a, int &ifail)
double BesselK0S(const double xx)
void Deqinv(const int n, std::vector< std::vector< double > > &a, int &ifail, std::vector< double > &b)
double BesselK1L(const double xx)
DoubleAc cos(const DoubleAc &f)
DoubleAc pow(const DoubleAc &f, double p)
DoubleAc exp(const DoubleAc &f)
DoubleAc fabs(const DoubleAc &f)
DoubleAc sin(const DoubleAc &f)
DoubleAc sqrt(const DoubleAc &f)