Garfield++ 4.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Macros | Functions
Isles.c File Reference
#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "Isles.h"

Go to the source code of this file.

Macros

#define DEFINE_ISLESGLOBAL
 
#define SHIFT   2.0
 
#define ARMAX   10000.0
 
#define ARMIN   0.0001
 
#define XNSegApprox   100
 
#define ZNSegApprox   100
 

Functions

int ExactRecSurf (double X, double Y, double Z, double xlo, double zlo, double xhi, double zhi, double *Potential, Vector3D *Flux)
 
int ApproxRecSurf (double X, double Y, double Z, double xlo, double zlo, double xhi, double zhi, int xseg, int zseg, double *Potential, Vector3D *Flux)
 
int ExactTriSurf (double zMax, double X, double Y, double Z, double *Potential, Vector3D *Flux)
 
int ApproxTriSurf (double zMax, double X, double Y, double Z, int nbxseg, int nbzseg, double *Potential, Vector3D *Flux)
 
double ExactCentroidalP_W (double rW, double lW)
 
double ExactAxialP_W (double rW, double lW, double Z)
 
double ExactAxialFZ_W (double rW, double lW, double Z)
 
double ApproxP_W (double rW, double lW, double X, double Y, double Z, int zseg)
 
double ApproxFX_W (double rW, double lW, double X, double Y, double Z, int zseg)
 
double ApproxFY_W (double rW, double lW, double X, double Y, double Z, int zseg)
 
double ApproxFZ_W (double rW, double lW, double X, double Y, double Z, int zseg)
 
int ApproxWire (double rW, double lW, double X, double Y, double Z, int zseg, double *potential, Vector3D *Flux)
 
double ImprovedP_W (double rW, double lW, double X, double Y, double Z)
 
double ImprovedFX_W (double rW, double lW, double X, double Y, double Z)
 
double ImprovedFY_W (double rW, double lW, double X, double Y, double Z)
 
double ImprovedFZ_W (double rW, double lW, double X, double Y, double Z)
 
int ImprovedWire (double rW, double lW, double X, double Y, double Z, double *potential, Vector3D *Flux)
 
double ExactThinP_W (double rW, double lW, double X, double Y, double Z)
 
double ExactThinFX_W (double rW, double lW, double X, double Y, double Z)
 
double ExactThinFY_W (double rW, double lW, double X, double Y, double Z)
 
double ExactThinFZ_W (double rW, double lW, double X, double Y, double Z)
 
int ExactThinWire (double rW, double lW, double X, double Y, double Z, double *potential, Vector3D *Flux)
 
double PointKnChPF (Point3D SourcePt, Point3D FieldPt, Vector3D *globalF)
 
double LineKnChPF (Point3D LineStart, Point3D LineStop, double radius, Point3D FieldPt, Vector3D *globalF)
 
double AreaKnChPF (int NbVertices, Point3D *Vertex, Point3D FieldPt, Vector3D *globalF)
 
double VolumeKnChPF (int, Point3D *, Point3D, Vector3D *globalF)
 
int Sign (double x)
 

Macro Definition Documentation

◆ ARMAX

#define ARMAX   10000.0

Definition at line 15 of file Isles.c.

◆ ARMIN

#define ARMIN   0.0001

Definition at line 16 of file Isles.c.

◆ DEFINE_ISLESGLOBAL

#define DEFINE_ISLESGLOBAL

Definition at line 5 of file Isles.c.

◆ SHIFT

#define SHIFT   2.0

Definition at line 14 of file Isles.c.

◆ XNSegApprox

#define XNSegApprox   100

Definition at line 20 of file Isles.c.

◆ ZNSegApprox

#define ZNSegApprox   100

Definition at line 21 of file Isles.c.

Function Documentation

◆ ApproxFX_W()

double ApproxFX_W ( double  rW,
double  lW,
double  X,
double  Y,
double  Z,
int  zseg 
)

Definition at line 1833 of file Isles.c.

1834 {
1835 if (DebugISLES) printf("In ApproxFX_W ...\n");
1836 ++ApproxCntr;
1837
1838 double dz = lW / zseg;
1839 double area = 2.0 * ST_PI * rW * dz;
1840
1841 double Fx = 0.0;
1842 double z0 = -0.5 * lW + 0.5 * dz;
1843 for (int k = 1; k <= zseg; ++k) {
1844 double zk = z0 + (k - 1) * dz;
1845 double dist = sqrt(X * X + Y * Y + (Z - zk) * (Z - zk));
1846 double dist3 = pow(dist, 3.0);
1847 if (fabs(dist) >= MINDIST) {
1848 Fx += (area * X / dist3);
1849 }
1850 } // zseg
1851
1852 return (Fx);
1853} // ApproxFX_W ends
ApproxCntr
ISLESGLOBAL int ApproxCntr
Definition: Isles.h:34
DebugISLES
ISLESGLOBAL int DebugISLES
Definition: Isles.h:35
ST_PI
#define ST_PI
Definition: Isles.h:24
MINDIST
#define MINDIST
Definition: Isles.h:17
Heed::pow
DoubleAc pow(const DoubleAc &f, double p)
Definition: DoubleAc.cpp:337
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314

◆ ApproxFY_W()

double ApproxFY_W ( double  rW,
double  lW,
double  X,
double  Y,
double  Z,
int  zseg 
)

Definition at line 1855 of file Isles.c.

1856 {
1857 if (DebugISLES) printf("In ApproxFY_W ...\n");
1858 ++ApproxCntr;
1859
1860 double dz = lW / zseg;
1861 double area = 2.0 * ST_PI * rW * dz;
1862
1863 double Fy = 0.0;
1864 double z0 = -0.5 * lW + 0.5 * dz;
1865 for (int k = 1; k <= zseg; ++k) {
1866 double zk = z0 + (k - 1) * dz;
1867 double dist = sqrt(X * X + Y * Y + (Z - zk) * (Z - zk));
1868 double dist3 = pow(dist, 3.0);
1869 if (fabs(dist) >= MINDIST) {
1870 Fy += (area * X / dist3);
1871 }
1872 } // zseg
1873
1874 return (Fy);
1875} // ApproxFY_W ends

◆ ApproxFZ_W()

double ApproxFZ_W ( double  rW,
double  lW,
double  X,
double  Y,
double  Z,
int  zseg 
)

Definition at line 1877 of file Isles.c.

1878 {
1879 if (DebugISLES) printf("In ApproxFZ_W ...\n");
1880 ++ApproxCntr;
1881
1882 double dz = lW / zseg;
1883 double area = 2.0 * ST_PI * rW * dz;
1884
1885 double Fz = 0.0;
1886 double z0 = -0.5 * lW + 0.5 * dz;
1887 for (int k = 1; k <= zseg; ++k) {
1888 double zk = z0 + (k - 1) * dz;
1889 double dist = sqrt(X * X + Y * Y + (Z - zk) * (Z - zk));
1890 double dist3 = pow(dist, 3.0);
1891 if (fabs(dist) >= MINDIST) {
1892 Fz += (area * X / dist3);
1893 }
1894 } // zseg
1895
1896 return (Fz);
1897} // ApproxFZ_W ends

◆ ApproxP_W()

double ApproxP_W ( double  rW,
double  lW,
double  X,
double  Y,
double  Z,
int  zseg 
)

Definition at line 1812 of file Isles.c.

1812 {
1813 // Approximate potential due to a wire segment
1814 if (DebugISLES) printf("In ApproxP_W ...\n");
1815 ++ApproxCntr;
1816
1817 double dz = lW / zseg;
1818 double area = 2.0 * ST_PI * rW * dz;
1819
1820 double Pot = 0.0;
1821 double z0 = -0.5 * lW + 0.5 * dz;
1822 for (int k = 1; k <= zseg; ++k) {
1823 double zk = z0 + (k - 1) * dz;
1824 double dist = sqrt(X * X + Y * Y + (Z - zk) * (Z - zk));
1825 if (fabs(dist) >= MINDIST) {
1826 Pot += area / dist;
1827 }
1828 } // zseg
1829
1830 return (Pot);
1831} // ApproxP_W ends

◆ ApproxRecSurf()

int ApproxRecSurf ( double  X,
double  Y,
double  Z,
double  xlo,
double  zlo,
double  xhi,
double  zhi,
int  xseg,
int  zseg,
double *  Potential,
Vector3D *  Flux 
)

Definition at line 706 of file Isles.c.

708 {
709
710 if (DebugISLES) {
711 printf("In ApproxRecSurf ...\n");
712 }
713
714 ++ApproxCntr;
715
716 double dx = (xhi - xlo) / xseg;
717 double dz = (zhi - zlo) / zseg;
718 double xel = (xhi - xlo) / xseg;
719 double zel = (zhi - zlo) / zseg;
720 double diag = sqrt(dx * dx + dz * dz);
721 double area = xel * zel;
722
723 double Pot = 0., XFlux = 0., YFlux = 0., ZFlux = 0.;
724
725 if (area > MINDIST2) { // else not necessary
726 for (int i = 1; i <= xseg; ++i) {
727 double xi = xlo + (dx / 2.0) + (i - 1) * dx;
728 for (int k = 1; k <= zseg; ++k) {
729 double zk = zlo + (dz / 2.0) + (k - 1) * dz;
730
731 double dist = sqrt((X - xi) * (X - xi) + Y * Y + (Z - zk) * (Z - zk));
732 if (DebugISLES) printf("dist: %lg\n", dist);
733 if (dist >= diag) {
734 Pot += area / dist;
735 } else if (dist <= MINDIST) {
736 // Self influence
737 Pot += 2.0 * (xel * log((zel + sqrt(xel * xel + zel * zel)) / xel) +
738 zel * log((xel + sqrt(xel * xel + zel * zel)) / zel));
739 } else {
740 // in the intermediate region where diag > dist > MINDIST
741 Pot += area / diag; // replace by expression of self-influence
742 if (DebugISLES) printf("Special Pot: %lg\n", area / diag);
743 }
744
745 if (dist >= diag) {
746 double f = area / (dist * dist * dist);
747 XFlux += f * (X - xi);
748 YFlux += f * Y;
749 ZFlux += f * (Z - zk);
750 } else {
751 double f = area / (diag * diag * diag);
752 XFlux += f * (X - xi);
753 YFlux += f * Y;
754 ZFlux += f * (Z - zk);
755 if (DebugISLES) {
756 printf("Special XFlux: %lg, YFlux: %lg, ZFlux: %lg\n",
757 f * (X - xi), f * Y, f * (Z - zk));
758 }
759 } // else dist >= diag
760 } // zseg
761 } // xseg
762 } // if area > MINDIST2
763
764 *Potential = Pot;
765 Flux->X = XFlux;
766 Flux->Y = YFlux;
767 Flux->Z = ZFlux;
768
769 return 0;
770} // ApproxRecSurf ends
MINDIST2
#define MINDIST2
Definition: Isles.h:18
Vector3D::Z
double Z
Definition: Vector.h:31
Vector3D::Y
double Y
Definition: Vector.h:30
Vector3D::X
double X
Definition: Vector.h:29

Referenced by ExactRecSurf().

◆ ApproxTriSurf()

int ApproxTriSurf ( double  zMax,
double  X,
double  Y,
double  Z,
int  nbxseg,
int  nbzseg,
double *  Potential,
Vector3D *  Flux 
)

Definition at line 1660 of file Isles.c.

1661 {
1662
1663 if (DebugISLES) printf("In ApproxTriSurf ...\n");
1664 ++ApproxCntr;
1665
1666 if (DebugISLES) {
1667 printf("zMax: %lg, X: %lg, Y: %lg, Z: %lg\n", zMax, X, Y, Z);
1668 printf("nbxseg: %d, nbzseg: %d\n", nbxseg, nbzseg);
1669 }
1670
1671 double dx = (1.0 - 0.0) / nbxseg;
1672 double dz = (zMax - 0.0) / nbzseg;
1673 double diag = sqrt(dx * dx + dz * dz);
1674 if (DebugISLES) printf("dx: %lg, dz: %lg, diag: %lg\n", dx, dz, diag);
1675 if ((dx < MINDIST) || (dz < MINDIST)) {
1676 printf("sub-element size too small in ApproxTriSurf.\n");
1677 return -1;
1678 }
1679
1680 double grad = (zMax - 0.0) / (1.0 - 0.0);
1681 if (DebugISLES) printf("grad: %lg\n", grad);
1682
1683 double Pot = 0., XFlux = 0., YFlux = 0., ZFlux = 0.;
1684 double Area = 0.0; // Total area of the element - useful for cross-check
1685 for (int i = 1; i <= nbxseg; ++i) {
1686 double xbgn = (i - 1) * dx;
1687 double zlimit_xbgn = zMax - grad * xbgn;
1688 double xend = i * dx;
1689 double zlimit_xend = zMax - grad * xend;
1690 if (DebugISLES)
1691 printf("i: %d, xbgn: %lg, zlimit_xbgn: %lg, xend: %lg, zlimit_xend:%lg\n",
1692 i, xbgn, zlimit_xbgn, xend, zlimit_xend);
1693
1694 for (int k = 1; k <= nbzseg; ++k) {
1695 double zbgn = (k - 1) * dz;
1696 double zend = k * dz;
1697 if (DebugISLES) printf("k: %d, zbgn: %lg, zend: %lg\n", k, zbgn, zend);
1698 int type_subele = 0;
1699 double area = 0.;
1700 double xi = 0., zk = 0.;
1701 if (zbgn >= zlimit_xbgn) {
1702 // completely outside the element - no effect of sub-element
1703 type_subele = 0;
1704 area = 0.0;
1705 } else if (zend <= zlimit_xend) {
1706 // completely within the element - rectangular sub-element
1707 type_subele = 1;
1708 xi = xbgn + 0.5 * dx;
1709 zk = zbgn + 0.5 * dz;
1710 area = dx * dz;
1711 } else if ((zbgn <= zlimit_xend) && (zend >= zlimit_xend)) {
1712 // partially inside the triangle
1713 type_subele = 2;
1714 // refer to the trapezoid figure
1715 double a = (zlimit_xend - zbgn);
1716 double b = (zlimit_xbgn - zbgn);
1717 double h = dx;
1718
1719 if (fabs(a) <= MINDIST) {
1720 // centroid of a triangle
1721 type_subele = 3;
1722 xi = xbgn + (h / 3.0);
1723 zk = zbgn + (b / 3.0);
1724 area = 0.5 * b * h;
1725 } else {
1726 // centroid of the trapezoid
1727 xi = xbgn + (h * (2. * a + b)) / (3. * (a + b));
1728 zk = zbgn + (a * a + a * b + b * b) / (3. * (a + b));
1729 area = 0.5 * h * (a + b);
1730 }
1731 } else {
1732 // takes care of round-off issues
1733 type_subele = 4;
1734 area = 0.0;
1735 }
1736 if (DebugISLES) printf("type_subele: %d, area: %lg\n", type_subele, area);
1737
1738 Area += area;
1739 if (area <= MINDIST2) continue;
1740 double dist = sqrt((X - xi) * (X - xi) + Y * Y + (Z - zk) * (Z - zk));
1741 if (DebugISLES) printf("dist: %lg\n", dist);
1742 if (dist >= diag) {
1743 Pot += area / dist;
1744 double f = area / (dist * dist * dist);
1745 XFlux += f * (X - xi);
1746 YFlux += f * Y;
1747 ZFlux += f * (Z - zk);
1748 } else {
1749 Pot += area / diag; // replace by expression of self-influence
1750 if (DebugISLES) printf("Special Pot: %lg\n", area / diag);
1751 double f = area / (diag * diag * diag);
1752 XFlux += f * (X - xi);
1753 YFlux += f * Y;
1754 ZFlux += f * (Z - zk);
1755 if (DebugISLES) {
1756 printf("Special XFlux: %lg, YFlux: %lg, ZFlux: %lg\n",
1757 f * (X - xi), f * Y, f * (Z - zk));
1758 }
1759 }
1760 } // nbzseg
1761 } // nbxseg
1762
1763 *Potential = Pot;
1764 Flux->X = XFlux;
1765 Flux->Y = YFlux;
1766 Flux->Z = ZFlux;
1767
1768 return 0;
1769} // ApproxTriSurf ends
Garfield::Polygon::Area
double Area(const std::vector< double > &xp, const std::vector< double > &yp)
Determine the (signed) area of a polygon.
Definition: Polygon.cc:274

Referenced by ExactTriSurf().

◆ ApproxWire()

int ApproxWire ( double  rW,
double  lW,
double  X,
double  Y,
double  Z,
int  zseg,
double *  potential,
Vector3D *  Flux 
)

Definition at line 1899 of file Isles.c.

1900 {
1901 if (DebugISLES) printf("In ApproxWire ...\n");
1902
1903 ++ApproxCntr;
1904
1905 double dz = lW / zseg;
1906 double area = 2.0 * ST_PI * rW * dz;
1907
1908 double Pot = 0., Fx = 0., Fy = 0., Fz = 0.;
1909 double z0 = -0.5 * lW + 0.5 * dz;
1910 for (int k = 1; k <= zseg; ++k) {
1911 double zk = z0 + (k - 1) * dz;
1912 double dist = sqrt(X * X + Y * Y + (Z - zk) * (Z - zk));
1913 double dist3 = pow(dist, 3.0);
1914 if (fabs(dist) >= MINDIST) {
1915 Pot += area / dist;
1916 Fx += (area * X / dist3);
1917 Fy += (area * Y / dist3);
1918 Fz += (area * Z / dist3);
1919 }
1920 } // zseg
1921
1922 *potential = Pot;
1923 Flux->X = Fx;
1924 Flux->Y = Fy;
1925 Flux->Z = Fz;
1926
1927 return 0;
1928}

◆ AreaKnChPF()

double AreaKnChPF ( int  NbVertices,
Point3D *  Vertex,
Point3D  FieldPt,
Vector3D *  globalF 
)

Definition at line 2350 of file Isles.c.

2351 {
2352 if (NbVertices != 4) {
2353 printf(
2354 "Only rectangular areas with known charges allowed at present ...\n");
2355 exit(-1);
2356 }
2357
2358 double xvert[4], yvert[4], zvert[4];
2359 xvert[0] = Vertex[0].X;
2360 xvert[1] = Vertex[1].X;
2361 xvert[2] = Vertex[2].X;
2362 xvert[3] = Vertex[3].X;
2363 yvert[0] = Vertex[0].Y;
2364 yvert[1] = Vertex[1].Y;
2365 yvert[2] = Vertex[2].Y;
2366 yvert[3] = Vertex[3].Y;
2367 zvert[0] = Vertex[0].Z;
2368 zvert[1] = Vertex[1].Z;
2369 zvert[2] = Vertex[2].Z;
2370 zvert[3] = Vertex[3].Z;
2371
2372 double xfld = FieldPt.X;
2373 double yfld = FieldPt.Y;
2374 double zfld = FieldPt.Z;
2375
2376 double xorigin = 0.25 * (xvert[3] + xvert[2] + xvert[1] + xvert[0]);
2377 double yorigin = 0.25 * (yvert[3] + yvert[2] + yvert[1] + yvert[0]);
2378 double zorigin = 0.25 * (zvert[3] + zvert[2] + zvert[1] + zvert[0]);
2379
2380 // lengths of the sides
2381 double SurfLX = sqrt((xvert[1] - xvert[0]) * (xvert[1] - xvert[0]) +
2382 (yvert[1] - yvert[0]) * (yvert[1] - yvert[0]) +
2383 (zvert[1] - zvert[0]) * (zvert[1] - zvert[0]));
2384 double SurfLZ = sqrt((xvert[2] - xvert[1]) * (xvert[2] - xvert[1]) +
2385 (yvert[2] - yvert[1]) * (yvert[2] - yvert[1]) +
2386 (zvert[2] - zvert[1]) * (zvert[2] - zvert[1]));
2387
2388 // Direction cosines - CHECK!
2389 //
2390 // 3 2
2391 // ________
2392 // | |
2393 // | |
2394 // | |
2395 // -------- -> +ve x-axis
2396 // 0 1
2397 // |
2398 // |
2399 // V +ve z-axis
2400 // The resulting +ve y-axis is out of the paper towards the reader.
2401 DirnCosn3D DirCos;
2402 DirCos.XUnit.X = (xvert[1] - xvert[0]) / SurfLX;
2403 DirCos.XUnit.Y = (yvert[1] - yvert[0]) / SurfLX;
2404 DirCos.XUnit.Z = (zvert[1] - zvert[0]) / SurfLX;
2405 DirCos.ZUnit.X = (xvert[0] - xvert[3]) / SurfLZ;
2406 DirCos.ZUnit.Y = (yvert[0] - yvert[3]) / SurfLZ;
2407 DirCos.ZUnit.Z = (zvert[0] - zvert[3]) / SurfLZ;
2408 DirCos.YUnit = Vector3DCrossProduct(&DirCos.ZUnit, &DirCos.XUnit);
2409
2410 Point3D localPt;
2411 // Through InitialVector[], field point gets translated to ECS origin.
2412 // Axes direction are, however, still global which when rotated to ECS
2413 // system, yields FinalVector[].
2414 { // Rotate point3D from global to local system to get localPt.
2415 double InitialVector[4];
2416 double TransformationMatrix[4][4] = {{0.0, 0.0, 0.0, 0.0},
2417 {0.0, 0.0, 0.0, 0.0},
2418 {0.0, 0.0, 0.0, 0.0},
2419 {0.0, 0.0, 0.0, 1.0}};
2420 double FinalVector[4];
2421
2422 InitialVector[0] = xfld - xorigin;
2423 InitialVector[1] = yfld - yorigin;
2424 InitialVector[2] = zfld - zorigin;
2425 InitialVector[3] = 1.0;
2426
2427 TransformationMatrix[0][0] = DirCos.XUnit.X;
2428 TransformationMatrix[0][1] = DirCos.XUnit.Y;
2429 TransformationMatrix[0][2] = DirCos.XUnit.Z;
2430 TransformationMatrix[1][0] = DirCos.YUnit.X;
2431 TransformationMatrix[1][1] = DirCos.YUnit.Y;
2432 TransformationMatrix[1][2] = DirCos.YUnit.Z;
2433 TransformationMatrix[2][0] = DirCos.ZUnit.X;
2434 TransformationMatrix[2][1] = DirCos.ZUnit.Y;
2435 TransformationMatrix[2][2] = DirCos.ZUnit.Z;
2436
2437 for (int i = 0; i < 4; ++i) {
2438 FinalVector[i] = 0.0;
2439 for (int j = 0; j < 4; ++j) {
2440 FinalVector[i] += TransformationMatrix[i][j] * InitialVector[j];
2441 }
2442 }
2443
2444 localPt.X = FinalVector[0];
2445 localPt.Y = FinalVector[1];
2446 localPt.Z = FinalVector[2];
2447 } // Point3D rotated
2448
2449 double Pot;
2450 Vector3D localF;
2451 double xpt = localPt.X;
2452 double ypt = localPt.Y;
2453 double zpt = localPt.Z;
2454
2455 double a = SurfLX;
2456 double b = SurfLZ;
2457 double diag = sqrt(a * a + b * b); // diagonal of the element
2458
2459 // distance of field point from element centroid
2460 double dist = sqrt(xpt * xpt + ypt * ypt + zpt * zpt);
2461 double dist3 = dist * dist * dist;
2462
2463 if (dist >= FarField * diag) // all are distances and, hence, +ve
2464 {
2465 double dA = a * b;
2466 Pot = dA / dist;
2467 localF.X = xpt * dA / dist3;
2468 localF.Y = ypt * dA / dist3;
2469 localF.Z = zpt * dA / dist3;
2470 } else {
2471 // normalize distances by `a' while sending - likely to improve accuracy
2472 int fstatus =
2473 ExactRecSurf(xpt / a, ypt / a, zpt / a, -1.0 / 2.0, -(b / a) / 2.0,
2474 1.0 / 2.0, (b / a) / 2.0, &Pot, &localF);
2475 if (fstatus) // non-zero
2476 {
2477 printf("problem in computing Potential of rectangular element ... \n");
2478 printf("a: %lg, b: %lg, X: %lg, Y: %lg, Z: %lg\n", a, b, xpt, ypt, zpt);
2479 // printf("returning ...\n");
2480 // return -1; void function at present
2481 }
2482 Pot *= a; // rescale Potential - cannot be done outside because of if(dist)
2483 }
2484
2485 (*globalF) = RotateVector3D(&localF, &DirCos, local2global);
2486 return (Pot);
2487} // AreaKnChPF ends
ExactRecSurf
int ExactRecSurf(double X, double Y, double Z, double xlo, double zlo, double xhi, double zhi, double *Potential, Vector3D *Flux)
Definition: Isles.c:40
FarField
#define FarField
Definition: Isles.h:21
RotateVector3D
Vector3D RotateVector3D(Vector3D *A, DirnCosn3D *DC, int Sense)
Definition: Vector.c:397
Vector3DCrossProduct
Vector3D Vector3DCrossProduct(Vector3D *A, Vector3D *B)
Definition: Vector.c:246
local2global
#define local2global
Definition: Vector.h:14
NbVertices
neBEMGLOBAL int * NbVertices
Definition: neBEM.h:70
DirnCosn3D::ZUnit
Vector3D ZUnit
Definition: Vector.h:38
DirnCosn3D::YUnit
Vector3D YUnit
Definition: Vector.h:37
DirnCosn3D::XUnit
Vector3D XUnit
Definition: Vector.h:36
Point3D
Definition: Vector.h:21
Point3D::X
double X
Definition: Vector.h:22
Point3D::Z
double Z
Definition: Vector.h:24
Point3D::Y
double Y
Definition: Vector.h:23

Referenced by KnChPFAtPoint().

◆ ExactAxialFZ_W()

double ExactAxialFZ_W ( double  rW,
double  lW,
double  Z 
)

Definition at line 1801 of file Isles.c.

1801 {
1802 if (DebugISLES) printf("In ExactAxialFZ_W ...\n");
1803 double h = 0.5 * lW;
1804 double Fz = 2.0 * ST_PI *
1805 (sqrt(h * h + 2.0 * Z * h + Z * Z + rW * rW) -
1806 sqrt(h * h - 2.0 * Z * h + Z * Z + rW * rW)) /
1807 sqrt(h * h - 2.0 * Z * h + Z * Z + rW * rW) /
1808 sqrt(h * h + 2.0 * Z * h + Z * Z + rW * rW);
1809 return (rW * Fz);
1810} // ExactAxialFZ ends

◆ ExactAxialP_W()

double ExactAxialP_W ( double  rW,
double  lW,
double  Z 
)

Definition at line 1789 of file Isles.c.

1789 {
1790 if (DebugISLES) printf("In ExactAxialP_W ...\n");
1791 // Expressions from PF00Z.m
1792 const double h = 0.5 * lW;
1793 const double r2 = rW * rW;
1794 const double a = h + Z;
1795 const double b = h - Z;
1796 double Pot = log((a + sqrt(a * a + r2)) * (b + sqrt(b * b + r2)) / r2);
1797 return 2. * ST_PI * rW * Pot;
1798} // ExactAxialP ends

Referenced by LineKnChPF(), WirePF(), WirePot(), and WirePrimPF().

◆ ExactCentroidalP_W()

double ExactCentroidalP_W ( double  rW,
double  lW 
)

Definition at line 1779 of file Isles.c.

1779 {
1780 // Self-Influence: wire5.m
1781 if (DebugISLES) printf("In ExactCentroidalP_W ...\n");
1782 const double h = 0.5 * lW;
1783 double dtmp1 = hypot(rW, h);
1784 dtmp1 = log((dtmp1 + h) / (dtmp1 - h));
1785 return 2.0 * ST_PI * dtmp1 * rW;
1786} // ExactCentroidalP ends

Referenced by LineKnChPF(), WirePF(), WirePot(), and WirePrimPF().

◆ ExactRecSurf()

int ExactRecSurf ( double  X,
double  Y,
double  Z,
double  xlo,
double  zlo,
double  xhi,
double  zhi,
double *  Potential,
Vector3D *  Flux 
)

Definition at line 40 of file Isles.c.

41 {
42
43 if (DebugISLES) printf("In ExactRecSurf ...\n");
44
45 ++IslesCntr;
46 ++ExactCntr;
47
48 ApproxFlag = 0;
49
50 if ((fabs(xhi - xlo) < 3.0 * MINDIST) || (fabs(zhi - zlo) < 3.0 * MINDIST)) {
51 fprintf(stdout, "Element size too small! ... returning ...\n");
52 return -1;
53 }
54
55 double ar = fabs((zhi - zlo) / (xhi - xlo));
56 if (ar > ARMAX || ar < ARMIN) {
57 fprintf(stdout, "Element too thin! ... returning ...\n");
58 return -2;
59 }
60
61 if (fabs(X) < MINDIST) X = 0.0;
62 if (fabs(Y) < MINDIST) Y = 0.0;
63 if (fabs(Z) < MINDIST) Z = 0.0;
64
65 double dxlo = X - xlo; // zero at the X=xlo edge
66 if (fabs(dxlo) < MINDIST) dxlo = 0.0;
67 double dxhi = X - xhi; // zero at the X=xhi edge
68 if (fabs(dxhi) < MINDIST) dxhi = 0.0;
69 double dzlo = Z - zlo; // zero at the Z=zlo edge
70 if (fabs(dzlo) < MINDIST) dzlo = 0.0;
71 double dzhi = Z - zhi; // zero at the Z=zhi edge
72 if (fabs(dzhi) < MINDIST) dzhi = 0.0;
73
74 // These four parameters can never be zero except at the four corners where
75 // one of them can become zero. For example, at X=xlo, Y=0, Z=zlo, D11
76 // is zero but the others are nonzero.
77 double D11 = sqrt(dxlo * dxlo + Y * Y + dzlo * dzlo);
78 if (fabs(D11) < MINDIST) D11 = 0.0;
79 double D21 = sqrt(dxhi * dxhi + Y * Y + dzlo * dzlo);
80 if (fabs(D21) < MINDIST) D21 = 0.0;
81 double D12 = sqrt(dxlo * dxlo + Y * Y + dzhi * dzhi);
82 if (fabs(D12) < MINDIST) D12 = 0.0;
83 double D22 = sqrt(dxhi * dxhi + Y * Y + dzhi * dzhi);
84 if (fabs(D22) < MINDIST) D22 = 0.0;
85
86 // Parameters related to the Y terms
87 int S1 = Sign(dzlo);
88 int S2 = Sign(dzhi);
89 double modY = fabs(Y);
90 int SY = Sign(Y);
91 double I1 = dxlo * modY;
92 double I2 = dxhi * modY;
93 double R1 = Y * Y + dzlo * dzlo;
94 double R2 = Y * Y + dzhi * dzhi;
95 if (fabs(I1) < MINDIST2) I1 = 0.0;
96 if (fabs(I2) < MINDIST2) I2 = 0.0;
97 if (fabs(R1) < MINDIST2) R1 = 0.0;
98 if (fabs(R2) < MINDIST2) R2 = 0.0;
99
100
101 if (DebugISLES) {
102 fprintf(stdout, "X: %.16lg, Y: %.16lg, Z: %.16lg\n", X, Y, Z);
103 fprintf(stdout, "xlo: %.16lg, zlo: %.16lg, xhi: %.16lg, zhi: %.16lg\n", xlo,
104 zlo, xhi, zhi);
105 fprintf(stdout, "dxlo: %.16lg, dzlo: %.16lg, dxhi: %.16lg, dzhi: %.16lg\n",
106 dxlo, dzlo, dxhi, dzhi);
107 fprintf(stdout, "D11: %.16lg, D12: %.16lg, D21: %.16lg, D22: %.16lg\n", D11,
108 D12, D21, D22);
109 fprintf(stdout, "S1: %d, S2: %d, modY: %.16lg\n", S1, S2, modY);
110 fprintf(stdout, "I1: %.16lg, I2: %.16lg, R1: %.16lg, R2: %.16lg\n", I1, I2,
111 R1, R2);
112 fprintf(stdout, "MINDIST: %.16lg, MINDIST2: %.16lg, SHIFT: %.16lg\n",
113 MINDIST, MINDIST2, SHIFT);
114 fflush(stdout);
115 }
116
117 // Check for possible numerical difficuties and take care.
118 // Presently the idea is to shift the field point slightly to a 'safe'
119 // position. Note that a shift in Y does not work because the singularities
120 // are associated with D11-s and dxlo-s and their combinations. A Y shift can
121 // sometimes alleviate the problem, but it cannot gurantee a permanent1s
122 // solution.
123 // A better approach is to evaluate analytic expressions truly valid for
124 // these critical regions, especially to take care of the >= and <=
125 // comparisons. For the == case, both of the previous comparisons are true and
126 // it is hard to justify one choice over the other. On the other hand, there
127 // are closed form analytic solutions for the == cases, and the problem does
128 // not stem from round-off errors as in the > or <, but not quite == cases.
129 // Note that shifts that are exactly equal to MINDIST, can give rise to
130 // un-ending recursions, leading to Segmentation fault.
131 // Corners
132 // Four corners
133 // Averages over four point evaluations may be carried out (skipped at
134 // present) This, however, may be difficult since we'll have to make sure that
135 // the shift towards the element does not bring the point too close to the
136 // same, or another difficult-to-evaluate situation. One of the ways to ensure
137 // this is to make SHIFT large enough, but that is unreasnoable and will
138 // introduce large amount of error.
139 if ((fabs(D11) <= MINDIST)) {
140 // close to xlo, 0, zlo
141 if (DebugISLES) printf("fabs(D11) <= MINDIST ... ");
142 double X1 = X;
143 double Z1 = Z;
144 if ((X >= xlo) && (Z >= zlo)) {
145 // point on the element
146 if (DebugISLES) printf("Case 1\n");
147 X1 = X + SHIFT * MINDIST;
148 Z1 = Z + SHIFT * MINDIST;
149 } else if ((X <= xlo) && (Z >= zlo)) {
150 // field point outside the element
151 if (DebugISLES) printf("Case 2 ...\n");
152 X1 = X - SHIFT * MINDIST;
153 Z1 = Z + SHIFT * MINDIST;
154 } else if ((X >= xlo) && (Z <= zlo)) {
155 // field point outside the element
156 if (DebugISLES) printf("Case 3 ...\n");
157 X1 = X + SHIFT * MINDIST;
158 Z1 = Z - SHIFT * MINDIST;
159 } else if ((X <= xlo) && (Z <= zlo)) {
160 // field point outside the element
161 if (DebugISLES) printf("Case 4 ...\n");
162 X1 = X - SHIFT * MINDIST;
163 Z1 = Z - SHIFT * MINDIST;
164 }
165 double Pot1;
166 Vector3D Flux1;
167 ExactRecSurf(X1, Y, Z1, xlo, zlo, xhi, zhi, &Pot1, &Flux1);
168 *Potential = Pot1;
169 Flux->X = Flux1.X;
170 Flux->Y = Flux1.Y;
171 Flux->Z = Flux1.Z;
172 return 0;
173 }
174 if ((fabs(D21) <= MINDIST)) {
175 // close to xhi, 0, zlo
176 if (DebugISLES) printf("fabs(D21) <= MINDIST ... ");
177 double X1 = X;
178 double Z1 = Z;
179 if ((X >= xhi) && (Z >= zlo)) {
180 // point outside the element
181 if (DebugISLES) printf("Case 1 ...\n");
182 X1 = X + SHIFT * MINDIST;
183 Z1 = Z + SHIFT * MINDIST;
184 } else if ((X <= xhi) && (Z >= zlo)) {
185 // point on the element
186 if (DebugISLES) printf("Case 2 ...\n");
187 X1 = X - SHIFT * MINDIST;
188 Z1 = Z + SHIFT * MINDIST;
189 } else if ((X >= xhi) && (Z <= zlo)) {
190 // field point outside the element
191 if (DebugISLES) printf("Case 3 ...\n");
192 X1 = X + SHIFT * MINDIST;
193 Z1 = Z - SHIFT * MINDIST;
194 } else if ((X <= xhi) && (Z <= zlo)) {
195 // field point outside the element
196 if (DebugISLES) printf("Case 4 ...\n");
197 X1 = X - SHIFT * MINDIST;
198 Z1 = Z - SHIFT * MINDIST;
199 }
200 double Pot1;
201 Vector3D Flux1;
202 ExactRecSurf(X1, Y, Z1, xlo, zlo, xhi, zhi, &Pot1, &Flux1);
203 *Potential = Pot1;
204 Flux->X = Flux1.X;
205 Flux->Y = Flux1.Y;
206 Flux->Z = Flux1.Z;
207 return 0;
208 }
209 if ((fabs(D12) <= MINDIST)) {
210 // close to xlo, 0, zhi
211 if (DebugISLES) printf("fabs(D12) <= MINDIST ... ");
212 double X1 = X;
213 double Z1 = Z;
214 if ((X >= xlo) && (Z >= zhi)) {
215 // point outside the element
216 if (DebugISLES) printf("Case 1 ...\n");
217 X1 = X + SHIFT * MINDIST;
218 Z1 = Z + SHIFT * MINDIST;
219 } else if ((X <= xlo) && (Z >= zhi)) {
220 // field point outside the element
221 if (DebugISLES) printf("Case 2 ...\n");
222 X1 = X - SHIFT * MINDIST;
223 Z1 = Z + SHIFT * MINDIST;
224 } else if ((X >= xlo) && (Z <= zhi)) {
225 // field point on the element
226 if (DebugISLES) printf("Case 3 ...\n");
227 X1 = X + SHIFT * MINDIST;
228 Z1 = Z - SHIFT * MINDIST;
229 } else if ((X <= xlo) && (Z <= zhi)) {
230 // field point outside the element
231 if (DebugISLES) printf("Case 4 ...\n");
232 X1 = X - SHIFT * MINDIST;
233 Z1 = Z - SHIFT * MINDIST;
234 }
235 double Pot1;
236 Vector3D Flux1;
237 ExactRecSurf(X1, Y, Z1, xlo, zlo, xhi, zhi, &Pot1, &Flux1);
238 *Potential = Pot1;
239 Flux->X = Flux1.X;
240 Flux->Y = Flux1.Y;
241 Flux->Z = Flux1.Z;
242 return 0;
243 }
244 if ((fabs(D22) <= MINDIST)) {
245 // close to xhi, 0, zhi
246 if (DebugISLES) printf("fabs(D22) <= MINDIST ... ");
247 double X1 = X;
248 double Z1 = Z;
249 if ((X >= xhi) && (Z >= zhi)) {
250 // point outside the element
251 if (DebugISLES) printf("Case 1 ...\n");
252 X1 = X + SHIFT * MINDIST;
253 Z1 = Z + SHIFT * MINDIST;
254 } else if ((X <= xhi) && (Z >= zhi)) {
255 // field point outside the element
256 if (DebugISLES) printf("Case 2 ...\n");
257 X1 = X - SHIFT * MINDIST;
258 Z1 = Z + SHIFT * MINDIST;
259 } else if ((X >= xhi) && (Z <= zhi)) {
260 // field point outside the element
261 if (DebugISLES) printf("Case 3 ...\n");
262 X1 = X + SHIFT * MINDIST;
263 Z1 = Z - SHIFT * MINDIST;
264 } else if ((X <= xhi) && (Z <= zhi)) {
265 // field point on the element
266 if (DebugISLES) printf("Case 4 ...\n");
267 X1 = X - SHIFT * MINDIST;
268 Z1 = Z - SHIFT * MINDIST;
269 }
270 double Pot1;
271 Vector3D Flux1;
272 ExactRecSurf(X1, Y, Z1, xlo, zlo, xhi, zhi, &Pot1, &Flux1);
273 *Potential = Pot1;
274 Flux->X = Flux1.X;
275 Flux->Y = Flux1.Y;
276 Flux->Z = Flux1.Z;
277 return 0;
278 }
279 // Four edges - average over two points on two sides of the edge may be ok.
280 // Here also we'll have to make sure that the shift
281 // towards the element does not bring the point too close to the same, or
282 // another difficult-to-evaluate situation. One of the ways to ensure this
283 // is to make SHIFT large enough, but that is unreasonable and will introduce
284 // large amount of error.
285 if (fabs(dxlo) < MINDIST) {
286 // edge at x=xlo || to Z - axis
287 if (DebugISLES) printf("fabs(dxlo) < MINDIST ... ");
288 double X1 = X;
289 double X2 = X;
290 if (X >= xlo) {
291 // field point on +ve side of YZ plane
292 if (DebugISLES) printf("Case 1 ...\n");
293 X1 += SHIFT * MINDIST;
294 X2 -= SHIFT * MINDIST;
295 } else {
296 // field point on -ve side of YZ plane
297 if (DebugISLES) printf("Case 2 ...\n");
298 X1 -= SHIFT * MINDIST;
299 X2 += SHIFT * MINDIST;
300 }
301 double Pot1, Pot2;
302 Vector3D Flux1, Flux2;
303 ExactRecSurf(X1, Y, Z, xlo, zlo, xhi, zhi, &Pot1, &Flux1);
304 ExactRecSurf(X2, Y, Z, xlo, zlo, xhi, zhi, &Pot2, &Flux2);
305 *Potential = 0.5 * (Pot1 + Pot2);
306 Flux->X = 0.5 * (Flux1.X + Flux2.X);
307 Flux->Y = 0.5 * (Flux1.Y + Flux2.Y);
308 Flux->Z = 0.5 * (Flux1.Z + Flux2.Z);
309 return 0;
310 }
311 if (fabs(dzlo) < MINDIST) {
312 // edge at z=zlo, || to X axis
313 if (DebugISLES) printf("fabs(dzlo) < MINDIST ... ");
314 double Z1 = Z;
315 double Z2 = Z;
316 if (Z >= zlo) {
317 // field point on +ve side of XY plane
318 if (DebugISLES) printf("Case 1 ...\n");
319 Z1 += SHIFT * MINDIST;
320 Z2 -= SHIFT * MINDIST;
321 } else {
322 // field point on -ve side of XY plane
323 if (DebugISLES) printf("Case 2 ...\n");
324 Z1 -= SHIFT * MINDIST;
325 Z2 += SHIFT * MINDIST;
326 }
327 double Pot1, Pot2;
328 Vector3D Flux1, Flux2;
329 ExactRecSurf(X, Y, Z1, xlo, zlo, xhi, zhi, &Pot1, &Flux1);
330 ExactRecSurf(X, Y, Z2, xlo, zlo, xhi, zhi, &Pot2, &Flux2);
331 *Potential = 0.5 * (Pot1 + Pot2);
332 Flux->X = 0.5 * (Flux1.X + Flux2.X);
333 Flux->Y = 0.5 * (Flux1.Y + Flux2.Y);
334 Flux->Z = 0.5 * (Flux1.Z + Flux2.Z);
335 return 0;
336 }
337 if (fabs(dxhi) < MINDIST) {
338 // edge at x=xhi, || to Z axis
339 if (DebugISLES) printf("fabs(dxhi) < MINDIST ... ");
340 double X1 = X;
341 double X2 = X;
342 if (X >= xhi) {
343 // field point on +ve side of YZ plane
344 if (DebugISLES) printf("Case 1 ...\n");
345 X1 += SHIFT * MINDIST;
346 X2 -= SHIFT * MINDIST;
347 } else {
348 // field point on -ve side of YZ plane
349 if (DebugISLES) printf("Case 2 ...\n");
350 X1 -= SHIFT * MINDIST;
351 X2 += SHIFT * MINDIST;
352 }
353 double Pot1, Pot2;
354 Vector3D Flux1, Flux2;
355 ExactRecSurf(X1, Y, Z, xlo, zlo, xhi, zhi, &Pot1, &Flux1);
356 ExactRecSurf(X2, Y, Z, xlo, zlo, xhi, zhi, &Pot2, &Flux2);
357 *Potential = 0.5 * (Pot1 + Pot2);
358 Flux->X = 0.5 * (Flux1.X + Flux2.X);
359 Flux->Y = 0.5 * (Flux1.Y + Flux2.Y);
360 Flux->Z = 0.5 * (Flux1.Z + Flux2.Z);
361 return 0;
362 }
363 if (fabs(dzhi) < MINDIST) {
364 // edge at z=zhi || to X axis
365 if (DebugISLES) printf("fabs(dzhi) < MINDIST ... ");
366 double Z1 = Z;
367 double Z2 = Z;
368 if (Z >= zhi) {
369 // field point on +ve side of XY plane
370 if (DebugISLES) printf("Case 1 ...\n");
371 Z1 += SHIFT * MINDIST;
372 Z2 -= SHIFT * MINDIST;
373 } else {
374 // field point on -ve side of XY plane
375 if (DebugISLES) printf("Case 2 ...\n");
376 Z1 -= SHIFT * MINDIST;
377 Z2 += SHIFT * MINDIST;
378 }
379 double Pot1, Pot2;
380 Vector3D Flux1, Flux2;
381 ExactRecSurf(X, Y, Z1, xlo, zlo, xhi, zhi, &Pot1, &Flux1);
382 ExactRecSurf(X, Y, Z2, xlo, zlo, xhi, zhi, &Pot2, &Flux2);
383 *Potential = 0.5 * (Pot1 + Pot2);
384 Flux->X = 0.5 * (Flux1.X + Flux2.X);
385 Flux->Y = 0.5 * (Flux1.Y + Flux2.Y);
386 Flux->Z = 0.5 * (Flux1.Z + Flux2.Z);
387 return 0;
388 }
389
390 // Logarithmic weak singularities are possible.
391 // Checks to be performed for 0 or -ve denominators and also
392 // 0 and +ve numerators.
393 // Interestingly, 0/0 does not cause a problem.
394 double DZTerm1 = log((D11 - dzlo) / (D12 - dzhi));
395 double DZTerm2 = log((D21 - dzlo) / (D22 - dzhi));
396 double DXTerm1 = log((D11 - dxlo) / (D21 - dxhi));
397 double DXTerm2 = log((D12 - dxlo) / (D22 - dxhi));
398
399 if (DebugISLES) {
400 fprintf(
401 stdout,
402 "DZTerm1: %.16lg, DZTerm2: %.16lg, DXTerm1: %.16lg, DXTerm2: %.16lg\n",
403 DZTerm1, DZTerm2, DXTerm1, DXTerm2);
404 }
405 // Four conditions based on the logarithmic terms
406 if (isnan(DZTerm1) || isinf(DZTerm1)) {
407 ++FailureCntr;
408 --ExactCntr;
409 ApproxFlag = 1;
410 fprintf(fIsles, "DZTerm1 problem ... approximating: %d.\n", ApproxFlag);
411 if (DebugISLES)
412 fprintf(stdout, "DZTerm1 problem ... approximating: %d.\n", ApproxFlag);
413 return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox, ZNSegApprox,
414 Potential, Flux));
415 }
416 if (isnan(DZTerm2) || isinf(DZTerm2)) {
417 ++FailureCntr;
418 --ExactCntr;
419 ApproxFlag = 2;
420 fprintf(fIsles, "DZTerm2 problem ... approximating: %d.\n", ApproxFlag);
421 if (DebugISLES)
422 fprintf(stdout, "DZTerm2 problem ... approximating: %d.\n", ApproxFlag);
423 return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox, ZNSegApprox,
424 Potential, Flux));
425 }
426 if (isnan(DXTerm1) || isinf(DXTerm1)) {
427 ++FailureCntr;
428 --ExactCntr;
429 ApproxFlag = 3;
430 fprintf(fIsles, "DXTerm1 problem ... approximating: %d.\n", ApproxFlag);
431 if (DebugISLES)
432 fprintf(stdout, "DXTerm1 problem ... approximating: %d.\n", ApproxFlag);
433 return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox, ZNSegApprox,
434 Potential, Flux));
435 }
436 if (isnan(DXTerm2) || isinf(DXTerm2)) {
437 ++FailureCntr;
438 --ExactCntr;
439 ApproxFlag = 4;
440 fprintf(fIsles, "DXTerm2 problem ... approximating: %d.\n", ApproxFlag);
441 if (DebugISLES)
442 fprintf(stdout, "DXTerm2 problem ... approximating: %d.\n", ApproxFlag);
443 return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox, ZNSegApprox,
444 Potential, Flux));
445 }
446 // Four criticalities based on the tanhyperbolic terms
447 // Since gsl_complex_arctanh_real can have any real number as its argument,
448 // all the criticalities are related to gsl_complex_arctanh where the
449 // imaginary component is present. So, fabs(I1) > mindist and
450 // fabs(I2) > mindist are only being tested.
451 if (S1 != 0) {
452 if (fabs(I1) > MINDIST2) {
453 if (D11 * fabs(dzlo) < MINDIST2) {
454 ++FailureCntr;
455 --ExactCntr;
456 ApproxFlag = 5;
457 fprintf(fIsles, "S1-I1 problem ... approximating: %d.\n", ApproxFlag);
458 if (DebugISLES)
459 fprintf(stdout, "S1-I1 problem ... approximating: %d.\n", ApproxFlag);
460 return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox,
461 ZNSegApprox, Potential, Flux));
462 }
463 }
464 if (fabs(I2) > MINDIST2) {
465 if (D21 * fabs(dzlo) < MINDIST2) {
466 ++FailureCntr;
467 --ExactCntr;
468 ApproxFlag = 6;
469 fprintf(fIsles, "S1-I2 problem ... approximating: %d.\n", ApproxFlag);
470 if (DebugISLES)
471 fprintf(stdout, "S1-I2 problem ... approximating: %d.\n", ApproxFlag);
472 return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox,
473 ZNSegApprox, Potential, Flux));
474 }
475 }
476 }
477 if (S2 != 0) {
478 if (fabs(I1) > MINDIST2) {
479 if (D12 * fabs(dzhi) < MINDIST2) {
480 ++FailureCntr;
481 --ExactCntr;
482 ApproxFlag = 7;
483 fprintf(fIsles, "S2-I1 problem ... approximating: %d.\n", ApproxFlag);
484 if (DebugISLES)
485 fprintf(stdout, "S2-I1 problem ... approximating: %d.\n", ApproxFlag);
486 return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox,
487 ZNSegApprox, Potential, Flux));
488 }
489 }
490 if (fabs(I2) > MINDIST2) {
491 if (D22 * fabs(dzhi) < MINDIST2) {
492 ++FailureCntr;
493 --ExactCntr;
494 ApproxFlag = 8;
495 fprintf(fIsles, "S2-I2 problem ... approximating: %d.\n", ApproxFlag);
496 if (DebugISLES)
497 fprintf(stdout, "S2-I2 problem ... approximating: %d.\n", ApproxFlag);
498 return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox,
499 ZNSegApprox, Potential, Flux));
500 }
501 }
502 }
503
504 double sumTanTerms = 0.;
505 // The possibility of singularities for dzhi or dzlo (division by zero)
506 // is overridden by the fact that S1 or S2 becomes zero in such cases
507 // and the singularity is avoided.
508 if (S1 != 0) {
509 if (fabs(I1) > MINDIST2) {
510 double a = D11 * D11 * dzlo * dzlo;
511 double b = R1 * R1 + I1 * I1;
512 double tmp = -S1 * atan(2 * I1 * D11 * fabs(dzlo) / (a - b));
513 if (b > a) {
514 if ((X > xlo && Z > zlo) || (X < xlo && Z < zlo)) {
515 tmp -= ST_PI;
516 } else if ((X < xlo && Z > zlo) || (X > xlo && Z < zlo)) {
517 tmp += ST_PI;
518 }
519 }
520 sumTanTerms += tmp;
521 }
522
523 if (fabs(I2) > MINDIST2) {
524 double a = D21 * D21 * dzlo * dzlo;
525 double b = R1 * R1 + I2 * I2;
526 double tmp = -S1 * atan(2 * I2 * D21 * fabs(dzlo) / (a - b));
527 if (b > a) {
528 if ((X > xhi && Z > zlo) || (X < xhi && Z < zlo)) {
529 tmp -= ST_PI;
530 } else if ((X < xhi && Z > zlo) || (X > xhi && Z < zlo)) {
531 tmp += ST_PI;
532 }
533 }
534 sumTanTerms -= tmp;
535 }
536 }
537
538 if (S2 != 0) {
539 if (fabs(I1) > MINDIST2) {
540 double a = D12 * D12 * dzhi * dzhi;
541 double b = R2 * R2 + I1 * I1;
542 double tmp = -S2 * atan(2 * I1 * D12 * fabs(dzhi) / (a - b));
543 if (b > a) {
544 if ((X > xlo && Z > zhi) || (X < xlo && Z < zhi)) {
545 tmp -= ST_PI;
546 } else if ((X < xlo && Z > zhi) || (X > xlo && Z < zhi)) {
547 tmp += ST_PI;
548 }
549 }
550 sumTanTerms -= tmp;
551 }
552
553 if (fabs(I2) > MINDIST2) {
554 double a = D22 * D22 * dzhi * dzhi;
555 double b = R2 * R2 + I2 * I2;
556 double tmp = -S2 * atan(2 * I2 * D22 * fabs(dzhi) / (a - b));
557 if (b > a) {
558 if ((X > xhi && Z > zhi) || (X < xhi && Z < zhi)) {
559 tmp -= ST_PI;
560 } else if ((X < xhi && Z > zhi) || (X > xhi && Z < zhi)) {
561 tmp += ST_PI;
562 }
563 }
564 sumTanTerms += tmp;
565 }
566 }
567
568 sumTanTerms *= -0.5;
569 double Pot = -dxlo * DZTerm1 + dxhi * DZTerm2 + modY * sumTanTerms -
570 dzlo * DXTerm1 + dzhi * DXTerm2;
571 double Fx = DZTerm1 - DZTerm2;
572 double Fy = -SY * sumTanTerms;
573 double Fz = DXTerm1 - DXTerm2;
574 if (DebugISLES) {
575 printf("XTerms: %.16lg, YTerms: %.16lg, ZTerms: %.16lg\n",
576 -dxlo * DZTerm1 + dxhi * DZTerm2, modY * sumTanTerms,
577 -dzlo * DXTerm1 + dzhi * DXTerm2);
578 printf("Pot: %lf, Fx: %lf, Fy: %lf, Fz: %lf\n", Pot, Fx, Fy, Fz);
579 fflush(stdout);
580 }
581
582 // constants of integration
583 // The only logic for the Fy constant seems to be the fact that the
584 // potential has a negative of this constant
585 if (((X > (xlo + MINDIST)) && (X < (xhi - MINDIST))) &&
586 ((Z > (zlo + MINDIST)) && (Z < (zhi - MINDIST)))) {
587 Pot -= 2.0 * modY * ST_PI;
588 if (SY != 0)
589 Fy += 2.0 * (double)SY * ST_PI;
590 else
591 Fy = 2.0 * ST_PI;
592 }
593 if (DebugISLES) {
594 printf("Constants of integration added for potential and Fy.\n");
595 printf("Pot: %lf, Fx: %lf, Fy: %lf, Fz: %lf\n", Pot, Fx, Fy, Fz);
596 fflush(stdout);
597 }
598
599 // Error situations handled before returning the values
600 if ((Pot < 0.0) || (isnan(Pot) || isinf(Pot))) {
601 fprintf(fIsles, "\n--- Approximation in ExactRecSurf ---\n");
602 fprintf(fIsles, "Negative, nan or inf potential ... approximating!\n");
603 if (DebugISLES) {
604 fprintf(stdout, "\n--- Approximation in ExactRecSurf ---\n");
605 fprintf(stdout, "Negative, nan or inf potential ... approximating!\n");
606 }
607 fprintf(fIsles, "X: %.16lg, Y: %.16lg, Z: %.16lg\n", X, Y, Z);
608 fprintf(fIsles, "xlo: %.16lg, zlo: %.16lg, xhi: %.16lg, zhi: %.16lg\n", xlo,
609 zlo, xhi, zhi);
610 fprintf(fIsles, "dxlo: %.16lg, dzlo: %.16lg, dxhi: %.16lg, dzhi: %.16lg\n",
611 dxlo, dzlo, dxhi, dzhi);
612 fprintf(fIsles, "D11: %.16lg, D12: %.16lg, D21: %.16lg, D22: %.16lg\n", D11,
613 D12, D21, D22);
614 fprintf(fIsles, "modY: %.16lg\n", modY);
615 fprintf(fIsles, "I1: %.16lg, I2: %.16lg, R1: %.16lg, R2: %.16lg\n", I1, I2,
616 R1, R2);
617 fprintf(fIsles, "S1: %d, S2: %d\n", S1, S2);
618 fprintf(fIsles, "Computed Pot: %.16lg\n", Pot);
619 ApproxFlag = 13;
620 ++FailureCntr;
621 --ExactCntr;
622 fprintf(fIsles, "Approximating: %d.\n", ApproxFlag);
623 return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox, ZNSegApprox,
624 Potential, Flux));
625 }
626 if (isnan(Fx) || isinf(Fx)) {
627 fprintf(fIsles, "\n--- Approximation in ExactRecSurf ---\n");
628 fprintf(fIsles, "Nan or inf Fx ... approximating!\n");
629 if (DebugISLES) {
630 fprintf(stdout, "\n--- Approximation in ExactRecSurf ---\n");
631 fprintf(stdout, "Nan or inf Fx ... approximating!\n");
632 }
633 fprintf(fIsles, "X: %.16lg, Y: %.16lg, Z: %.16lg\n", X, Y, Z);
634 fprintf(fIsles, "xlo: %.16lg, zlo: %.16lg, xhi: %.16lg, zhi: %.16lg\n", xlo,
635 zlo, xhi, zhi);
636 fprintf(fIsles, "dxlo: %.16lg, dzlo: %.16lg, dxhi: %.16lg, dzhi: %.16lg\n",
637 dxlo, dzlo, dxhi, dzhi);
638 fprintf(fIsles, "D11: %.16lg, D12: %.16lg, D21: %.16lg, D22: %.16lg\n", D11,
639 D12, D21, D22);
640 fprintf(fIsles, "Computed Fx: %.16lg\n", Fx);
641 ApproxFlag = 14;
642 ++FailureCntr;
643 --ExactCntr;
644 fprintf(fIsles, "Approximating: %d.\n", ApproxFlag);
645 return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox, ZNSegApprox,
646 Potential, Flux));
647 }
648 if (isnan(Fy) || isinf(Fy)) {
649 fprintf(fIsles, "\n--- Approximation in ExactRecSurf ---\n");
650 fprintf(fIsles, "Nan or inf Fy ... approximating!\n");
651 if (DebugISLES) {
652 fprintf(stdout, "\n--- Approximation in ExactRecSurf ---\n");
653 fprintf(stdout, "Nan or inf Fy ... approximating!\n");
654 }
655 fprintf(fIsles, "X: %.16lg, Y: %.16lg, Z: %.16lg\n", X, Y, Z);
656 fprintf(fIsles, "xlo: %.16lg, zlo: %.16lg, xhi: %.16lg, zhi: %.16lg\n", xlo,
657 zlo, xhi, zhi);
658 fprintf(fIsles, "dxlo: %.16lg, dzlo: %.16lg, dxhi: %.16lg, dzhi: %.16lg\n",
659 dxlo, dzlo, dxhi, dzhi);
660 fprintf(fIsles, "D11: %.16lg, D12: %.16lg, D21: %.16lg, D22: %.16lg\n", D11,
661 D12, D21, D22);
662 fprintf(fIsles, "S1: %d, S2: %d, SY: %d, modY: %.16lg\n", S1, S2, SY, modY);
663 fprintf(fIsles, "I1: %.16lg, I2: %.16lg, R1: %.16lg, R2: %.16lg\n", I1, I2,
664 R1, R2);
665 fprintf(fIsles, "Computed Fy: %.16lg\n", Fy);
666 ApproxFlag = 15;
667 ++FailureCntr;
668 --ExactCntr;
669 fprintf(fIsles, "Approximating: %d.\n", ApproxFlag);
670 return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox, ZNSegApprox,
671 Potential, Flux));
672 }
673 if (isnan(Fz) || isinf(Fz)) {
674 fprintf(fIsles, "\n--- Approximation in ExactRecSurf ---\n");
675 fprintf(fIsles, "Nan or inf Fz ... approximating!\n");
676 if (DebugISLES) {
677 fprintf(stdout, "\n--- Approximation in ExactRecSurf ---\n");
678 fprintf(stdout, "Nan or inf Fz ... approximating!\n");
679 }
680 fprintf(fIsles, "X: %.16lg, Y: %.16lg, Z: %.16lg\n", X, Y, Z);
681 fprintf(fIsles, "xlo: %.16lg, zlo: %.16lg, xhi: %.16lg, zhi: %.16lg\n", xlo,
682 zlo, xhi, zhi);
683 fprintf(fIsles, "dxlo: %.16lg, dzlo: %.16lg, dxhi: %.16lg, dzhi: %.16lg\n",
684 dxlo, dzlo, dxhi, dzhi);
685 fprintf(fIsles, "D11: %.16lg, D12: %.16lg, D21: %.16lg, D22: %.16lg\n", D11,
686 D12, D21, D22);
687 fprintf(fIsles, "Computed Fz: %.16lg\n", Fz);
688 ApproxFlag = 16;
689 ++FailureCntr;
690 --ExactCntr;
691 fprintf(fIsles, "Approximating: %d.\n", ApproxFlag);
692 return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox, ZNSegApprox,
693 Potential, Flux));
694 }
695
696 *Potential = Pot;
697 Flux->X = Fx;
698 Flux->Y = Fy;
699 Flux->Z = Fz;
700
701 if (DebugISLES) printf("Going out of ExactRecSurf ...\n");
702
703 return 0;
704} // ExactRecSurf ends
ARMIN
#define ARMIN
Definition: Isles.c:16
ZNSegApprox
#define ZNSegApprox
Definition: Isles.c:21
ARMAX
#define ARMAX
Definition: Isles.c:15
Sign
int Sign(double x)
Definition: Isles.c:2506
SHIFT
#define SHIFT
Definition: Isles.c:14
ApproxRecSurf
int ApproxRecSurf(double X, double Y, double Z, double xlo, double zlo, double xhi, double zhi, int xseg, int zseg, double *Potential, Vector3D *Flux)
Definition: Isles.c:706
XNSegApprox
#define XNSegApprox
Definition: Isles.c:20
IslesCntr
ISLESGLOBAL int IslesCntr
Definition: Isles.h:34
fIsles
ISLESGLOBAL FILE * fIsles
Definition: Isles.h:36
FailureCntr
ISLESGLOBAL int FailureCntr
Definition: Isles.h:34
ApproxFlag
ISLESGLOBAL int ApproxFlag
Definition: Isles.h:35
ExactCntr
ISLESGLOBAL int ExactCntr
Definition: Isles.h:34

Referenced by AreaKnChPF(), ExactRecSurf(), RecFlux(), RecPF(), RecPot(), and RecPrimPF().

◆ ExactThinFX_W()

double ExactThinFX_W ( double  rW,
double  lW,
double  X,
double  Y,
double  Z 
)

Definition at line 2103 of file Isles.c.

2103 {
2104 if (DebugISLES) {
2105 printf("In ExactThinFX_W ...\n");
2106 printf("rW: %lg, lW: %lg, X: %lg, Y: %lg, Z: %lg\n", rW, lW, X, Y, Z);
2107 }
2108 // Expressions from PFXYZ_thin.m
2109 double h = 0.5 * lW;
2110 double Fx = 2.0 * X *
2111 (sqrt(X * X + Y * Y + Z * Z + h * h + 2.0 * Z * h) * h -
2112 sqrt(X * X + Y * Y + Z * Z + h * h + 2.0 * Z * h) * Z +
2113 sqrt(X * X + Y * Y + Z * Z - 2.0 * Z * h + h * h) * h +
2114 sqrt(X * X + Y * Y + Z * Z - 2.0 * Z * h + h * h) * Z) /
2115 (X * X + Y * Y) / sqrt(X * X + Y * Y + Z * Z - 2.0 * Z * h + h * h) /
2116 sqrt(X * X + Y * Y + Z * Z + h * h + 2.0 * Z * h) * ST_PI;
2117 return (rW * Fx);
2118} // ExactThinFX_W ends

Referenced by LineKnChPF(), and WireFlux().

◆ ExactThinFY_W()

double ExactThinFY_W ( double  rW,
double  lW,
double  X,
double  Y,
double  Z 
)

Definition at line 2121 of file Isles.c.

2121 {
2122 if (DebugISLES) {
2123 printf("In ExactThinFY_W ...\n");
2124 printf("rW: %lg, lW: %lg, X: %lg, Y: %lg, Z: %lg\n", rW, lW, X, Y, Z);
2125 }
2126 // Expressions from PFXYZ_thin.m
2127 double h = 0.5 * lW;
2128 double Fy = 2.0 * Y *
2129 (sqrt(X * X + Y * Y + Z * Z + h * h + 2.0 * Z * h) * h -
2130 sqrt(X * X + Y * Y + Z * Z + h * h + 2.0 * Z * h) * Z +
2131 sqrt(X * X + Y * Y + Z * Z - 2.0 * Z * h + h * h) * h +
2132 sqrt(X * X + Y * Y + Z * Z - 2.0 * Z * h + h * h) * Z) /
2133 (X * X + Y * Y) / sqrt(X * X + Y * Y + Z * Z - 2.0 * Z * h + h * h) /
2134 sqrt(X * X + Y * Y + Z * Z + h * h + 2.0 * Z * h) * ST_PI;
2135 return (rW * Fy);
2136} // ExactThinFY_W ends

Referenced by LineKnChPF(), and WireFlux().

◆ ExactThinFZ_W()

double ExactThinFZ_W ( double  rW,
double  lW,
double  X,
double  Y,
double  Z 
)

Definition at line 2139 of file Isles.c.

2139 {
2140 if (DebugISLES) {
2141 printf("In ExactThinFZ_W ...\n");
2142 printf("rW: %lg, lW: %lg, X: %lg, Y: %lg, Z: %lg\n", rW, lW, X, Y, Z);
2143 }
2144
2145 // Expressions from PFXYZ_thin.m
2146 double h = 0.5 * lW;
2147 double Fz = 2.0 *
2148 (sqrt(X * X + Y * Y + Z * Z + h * h + 2.0 * Z * h) -
2149 sqrt(X * X + Y * Y + Z * Z - 2.0 * Z * h + h * h)) /
2150 sqrt(X * X + Y * Y + Z * Z - 2.0 * Z * h + h * h) /
2151 sqrt(X * X + Y * Y + Z * Z + h * h + 2.0 * Z * h) * ST_PI;
2152 return (rW * Fz);
2153} // ExactThinFZ_W ends

Referenced by LineKnChPF(), WireFlux(), WirePF(), and WirePrimPF().

◆ ExactThinP_W()

double ExactThinP_W ( double  rW,
double  lW,
double  X,
double  Y,
double  Z 
)

Definition at line 2088 of file Isles.c.

2088 {
2089 if (DebugISLES) {
2090 printf("In ExactThinP_W ...\n");
2091 printf("rW: %lg, lW: %lg, X: %lg, Y: %lg, Z: %lg\n", rW, lW, X, Y, Z);
2092 }
2093 // Expressions from PFXYZ_thin.m
2094 const double h = 0.5 * lW;
2095 const double r2 = X * X + Y * Y;
2096 const double a = h + Z;
2097 const double b = h - Z;
2098 double Pot = log((a + sqrt(r2 + a * a)) * (b + sqrt(r2 + b * b)) / r2);
2099 return 2. * ST_PI * rW * Pot;
2100} // ExactThinP_W ends

Referenced by LineKnChPF(), and WirePot().

◆ ExactThinWire()

int ExactThinWire ( double  rW,
double  lW,
double  X,
double  Y,
double  Z,
double *  potential,
Vector3D *  Flux 
)

Definition at line 2155 of file Isles.c.

2156 {
2157 if (DebugISLES) {
2158 printf("In ExactThinWire_W ...\n");
2159 printf("rW: %lg, lW: %lg, X: %lg, Y: %lg, Z: %lg\n", rW, lW, X, Y, Z);
2160 }
2161 const double h = 0.5 * lW;
2162 const double r2 = X * X + Y * Y;
2163 const double a = h + Z;
2164 const double b = h - Z;
2165 const double c = sqrt(r2 + a * a);
2166 const double d = sqrt(r2 + b * b);
2167 const double s = 2. * ST_PI * rW;
2168 *potential = s * log((a + c) * (b + d) / r2);
2169 double f = s * (c * b + d * a) / (r2 * c * d);
2170 Flux->X = f * X;
2171 Flux->Y = f * Y;
2172 Flux->Z = s * (c - d) / (c * d);
2173 return 0;
2174}

Referenced by WirePF(), and WirePrimPF().

◆ ExactTriSurf()

int ExactTriSurf ( double  zMax,
double  X,
double  Y,
double  Z,
double *  Potential,
Vector3D *  Flux 
)

Definition at line 784 of file Isles.c.

785 {
786
787 if (DebugISLES) printf("In ExactTriSurf ...\n");
788
789 ++IslesCntr;
790 ++ExactCntr;
791 // Reset the flag to indicate approximate evaluation of potential.
792 ApproxFlag = 0;
793
794 // We do not need to check for X, since element extent is always 0 - 1.
795 if (zMax < 3.0 * SHIFT * MINDIST) {
796 // should allow enough space for Z corrections
797 // One SHIFT should not lead to another criticality
798 fprintf(stdout, "Element size too small! ... returning ...\n");
799 return -1;
800 }
801
802 if (zMax > ARMAX || zMax < ARMIN) {
803 fprintf(stdout, "Element too thin! ... returning ...\n");
804 return -1;
805 }
806
807 // These three parameters can never be zero except at the three corners where
808 // one of them can become zero. For example, at X=0, Y=0, Z=0, D11
809 // is zero but the others are nonzero.
810 if (fabs(X) < MINDIST) X = 0.0;
811 if (fabs(Y) < MINDIST) Y = 0.0;
812 if (fabs(Z) < MINDIST) Z = 0.0;
813 double modY = fabs(Y);
814 if (modY < MINDIST) modY = 0.0;
815 int S1 = Sign(Z);
816 int SY = Sign(Y);
817
818 // distances from corners (0,0,0), (1,0,0) and (0,0,zMax)
819 double D11 = sqrt(X * X + Y * Y + Z * Z);
820 if (D11 < MINDIST) D11 = 0.0;
821 double D21 = sqrt((X - 1.0) * (X - 1.0) + Y * Y + Z * Z);
822 if (D21 < MINDIST) D21 = 0.0;
823 double D12 = sqrt(X * X + Y * Y + (Z - zMax) * (Z - zMax));
824 if (D12 < MINDIST) D12 = 0.0;
825
826 double G = zMax * (X - 1.0) + Z;
827 if (fabs(G) < MINDIST) G = 0.0;
828 double E1 = (X - zMax * (Z - zMax));
829 if (fabs(E1) < MINDIST) E1 = 0.0;
830 double E2 = (X - 1.0 - zMax * Z);
831 if (fabs(E2) < MINDIST) E2 = 0.0;
832 double H1 = Y * Y + (Z - zMax) * G;
833 if (fabs(H1) < MINDIST2) H1 = 0.0;
834 double H2 = Y * Y + Z * G;
835 if (fabs(H2) < MINDIST2) H2 = 0.0;
836 double R1 = Y * Y + Z * Z;
837 if (fabs(R1) < MINDIST2) R1 = 0.0;
838 double I1 = modY * X;
839 if (fabs(I1) < MINDIST2) I1 = 0.0;
840 double I2 = modY * (X - 1.0);
841 if (fabs(I2) < MINDIST2) I2 = 0.0;
842 double Hypot = sqrt(1.0 + zMax * zMax);
843 if (Hypot < MINDIST) { // superfluous
844 fprintf(stdout, "Hypotenuse too small! ... returning ...\n");
845 return -1;
846 }
847 double Zhyp = zMax * (1.0 - X); // Z on hypotenuse (extended) for given X
848
849 if (DebugISLES) {
850 printf("\n\nzMax: %.16lg, X: %.16lg, Y: %.16lg, Z: %.16lg\n", zMax, X, Y,
851 Z);
852 printf("D11: %.16lg, D21: %.16lg, D12: %.16lg, Hypot: %.16lg\n", D11, D21,
853 D12, Hypot);
854 printf("modY: %.16lg, G: %.16lg\n", modY, G);
855 printf("E1: %.16lg, E2: %.16lg, H1: %.16lg, H2: %.16lg\n", E1, E2, H1, H2);
856 printf("S1: %d, SY: %d, R1: %.16lg, I1: %.16lg, I2: %.16lg\n", S1, SY, R1,
857 I1, I2);
858 printf("H1+G*D12: %.16lg, E1-zMax*D12: %.16lg\n", H1 + G * D12,
859 E1 - zMax * D12);
860 printf("H2+G*D21: %.16lg, E2-zMax*D21: %.16lg\n", H2 + G * D21,
861 E2 - zMax * D21);
862 printf("D11*fabs(Z): %.16lg, D21*fabs(Z): %.16lg\n", D11 * fabs(Z),
863 D21 * fabs(Z));
864 printf("Hypot*D12 - E1: %.16lg\n", Hypot * D12 - E1);
865 printf("Hypot*D21 - E2: %.16lg\n\n", Hypot * D21 - E2);
866 fprintf(stdout, "MINDIST: %.16lg, MINDIST2: %.16lg, SHIFT: %.16lg\n",
867 MINDIST, MINDIST2, SHIFT);
868 fflush(stdout);
869 }
870
871 // Check for possible numerical difficulties and take care.
872 // Presently the idea is to shift the field point slightly to a 'safe'
873 // position. Note that a shift in Y does not work because the singularities
874 // are associated with D11-s and dxlo-s and their combinations. A Y shift can
875 // sometimes alleviate the problem, but it cannot gurantee a permanent1s
876 // solution.
877 // A better approach is to evaluate analytic expressions truly valid for
878 // these critical regions, especially to take care of the >= and <=
879 // comparisons. For the == case, both of the previous comparisons are true and
880 // it is hard to justify one choice over the other. On the other hand, there
881 // are closed form analytic solutions for the == cases, and the problem does
882 // not stem from round-off errors as in the > or <, but not quite == cases.
883 // Possible singularity at D21 corner where X=1, Z=0
884 // Possible singularity at D12 corner where X=0, Z=zMax
885 // Check for possible numerical difficuties and take care
886 // Note that modY = 0 cannot be a condition of criticality since considering
887 // that would mean omitting even the barycenter properties.
888 // Also note that shifts that are exactly equal to MINDIST, can give rise to
889 // un-ending recursions, leading to Segmentation fault.
890 // Special points and combinations
891
892 // Three corners
893 if ((fabs(D11) <= MINDIST)) {
894 if (DebugISLES) printf("D11 <= MINDIST\n");
895 double X1 = X;
896 double Z1 = Z;
897 if ((X >= 0.0) && (Z >= 0.0)) {
898 // point on the element
899 if (DebugISLES) printf("Case1=>\n");
900 X1 = X + SHIFT * MINDIST;
901 Z1 = Z + SHIFT * MINDIST;
902 } else if ((X <= 0.0) && (Z >= 0.0)) {
903 // field point outside the element
904 if (DebugISLES) printf("Case2=>\n");
905 X1 = X - SHIFT * MINDIST;
906 Z1 = Z + SHIFT * MINDIST;
907 } else if ((X >= 0.0) && (Z <= 0.0)) {
908 // field point outside the element
909 if (DebugISLES) printf("Case3=>\n");
910 X1 = X + SHIFT * MINDIST;
911 Z1 = Z - SHIFT * MINDIST;
912 } else if ((X <= 0.0) && (Z <= 0.0)) {
913 // field point outside the element
914 if (DebugISLES) printf("Case4=>\n");
915 X1 = X - SHIFT * MINDIST;
916 Z1 = Z - SHIFT * MINDIST;
917 }
918 double Pot1;
919 Vector3D Flux1;
920 ExactTriSurf(zMax, X1, Y, Z1, &Pot1, &Flux1);
921 *Potential = Pot1;
922 Flux->X = Flux1.X;
923 Flux->Y = Flux1.Y;
924 Flux->Z = Flux1.Z;
925 return 0;
926 }
927 if ((fabs(D21) <= MINDIST)) {
928 if (DebugISLES) printf("D21 <= MINDIST\n");
929 double X1 = X;
930 double Z1 = Z;
931 if ((X >= 1.0) && (Z >= 0.0)) {
932 // point outside the element
933 if (DebugISLES) printf("Case1=>\n");
934 X1 = X + SHIFT * MINDIST;
935 Z1 = Z + SHIFT * MINDIST;
936 } else if ((X <= 1.0) && (Z >= 0.0)) {
937 // field point on the element
938 // difficult to decide (chk figure) - chk whether Z is beyond Zhyp
939 if (DebugISLES) printf("Case2=>\n");
940 X1 = X - SHIFT * MINDIST;
941 Z1 = Z + SHIFT * MINDIST;
942 } else if ((X >= 1.0) && (Z <= 0.0)) {
943 // field point outside the element
944 if (DebugISLES) printf("Case3=>\n");
945 X1 = X + SHIFT * MINDIST;
946 Z1 = Z - SHIFT * MINDIST;
947 } else if ((X <= 1.0) && (Z <= 0.0)) {
948 // field point outside the element
949 if (DebugISLES) printf("Case4=>\n");
950 X1 = X - SHIFT * MINDIST;
951 Z1 = Z - SHIFT * MINDIST;
952 }
953 double Pot1;
954 Vector3D Flux1;
955 ExactTriSurf(zMax, X1, Y, Z1, &Pot1, &Flux1);
956 *Potential = Pot1;
957 Flux->X = Flux1.X;
958 Flux->Y = Flux1.Y;
959 Flux->Z = Flux1.Z;
960 return 0;
961 }
962 if ((fabs(D12) <= MINDIST)) {
963 if (DebugISLES) printf("D12 <= MINDIST\n");
964 double X1 = X;
965 double Z1 = Z;
966 if ((X >= 0.0) && (Z >= zMax)) {
967 // point outside the element
968 if (DebugISLES) printf("Case1=>\n");
969 X1 = X + SHIFT * MINDIST;
970 Z1 = Z + SHIFT * MINDIST;
971 } else if ((X <= 0.0) && (Z >= zMax)) {
972 // field point outside the element
973 if (DebugISLES) printf("Case2=>\n");
974 X1 = X - SHIFT * MINDIST;
975 Z1 = Z + SHIFT * MINDIST;
976 } else if ((X >= 0.0) && (Z <= zMax)) {
977 // field point on the element
978 // can create problem for small element - chk whether Z is beyond Zhyp
979 if (DebugISLES) printf("Case3=>\n");
980 X1 = X + SHIFT * MINDIST;
981 Z1 = Z - SHIFT * MINDIST;
982 } else if ((X <= 0.0) && (Z <= zMax)) {
983 // field point outside the element
984 if (DebugISLES) printf("Case4=>\n");
985 X1 = X - SHIFT * MINDIST;
986 Z1 = Z - SHIFT * MINDIST;
987 }
988 double Pot1;
989 Vector3D Flux1;
990 ExactTriSurf(zMax, X1, Y, Z1, &Pot1, &Flux1);
991 *Potential = Pot1;
992 Flux->X = Flux1.X;
993 Flux->Y = Flux1.Y;
994 Flux->Z = Flux1.Z;
995 return 0;
996 }
997 // Three Y lines at corners
998 if ((fabs(X) < MINDIST) && (fabs(Z) < MINDIST)) {
999 // Y line at (0,0,0) corner
1000 if (DebugISLES) printf("Y line at (0,0,0) corner\n");
1001 double X1 = X;
1002 double Z1 = Z;
1003 if ((X >= 0.0) && (Z >= 0.0)) {
1004 // point on the element
1005 if (DebugISLES) printf("Case1=>\n");
1006 X1 = X + SHIFT * MINDIST;
1007 Z1 = Z + SHIFT * MINDIST;
1008 } else if ((X <= 0.0) && (Z >= 0.0)) {
1009 // field point outside the element
1010 if (DebugISLES) printf("Case2=>\n");
1011 X1 = X - SHIFT * MINDIST;
1012 Z1 = Z + SHIFT * MINDIST;
1013 } else if ((X >= 0.0) && (Z <= 0.0)) {
1014 // field point outside the element
1015 if (DebugISLES) printf("Case3=>\n");
1016 X1 = X + SHIFT * MINDIST;
1017 Z1 = Z - SHIFT * MINDIST;
1018 } else if ((X <= 0.0) && (Z <= 0.0)) {
1019 // field point outside the element
1020 if (DebugISLES) printf("Case4=>\n");
1021 X1 = X - SHIFT * MINDIST;
1022 Z1 = Z - SHIFT * MINDIST;
1023 }
1024 double Pot1;
1025 Vector3D Flux1;
1026 ExactTriSurf(zMax, X1, Y, Z1, &Pot1, &Flux1);
1027 *Potential = Pot1;
1028 Flux->X = Flux1.X;
1029 Flux->Y = Flux1.Y;
1030 Flux->Z = Flux1.Z;
1031 return 0;
1032 }
1033 if ((fabs(X - 1.0) < MINDIST) && (fabs(Z) < MINDIST)) {
1034 // Y line at (1,0) corner
1035 if (DebugISLES) printf("Y line at (1,0,0) corner\n");
1036 double X1 = X;
1037 double Z1 = Z;
1038 if ((X >= 1.0) && (Z >= 0.0)) {
1039 // point on the element
1040 if (DebugISLES) printf("Case1=>\n");
1041 X1 = X + SHIFT * MINDIST;
1042 Z1 = Z + SHIFT * MINDIST;
1043 } else if ((X <= 1.0) && (Z >= 0.0)) {
1044 // field point outside the element
1045 if (DebugISLES) printf("Case2=>\n");
1046 X1 = X - SHIFT * MINDIST;
1047 Z1 = Z + SHIFT * MINDIST;
1048 } else if ((X >= 1.0) && (Z <= 0.0)) {
1049 // field point outside the element
1050 if (DebugISLES) printf("Case3=>\n");
1051 X1 = X + SHIFT * MINDIST;
1052 Z1 = Z - SHIFT * MINDIST;
1053 } else if ((X <= 1.0) && (Z <= 0.0)) {
1054 // field point outside the element
1055 if (DebugISLES) printf("Case4=>\n");
1056 X1 = X - SHIFT * MINDIST;
1057 Z1 = Z - SHIFT * MINDIST;
1058 }
1059 double Pot1;
1060 Vector3D Flux1;
1061 ExactTriSurf(zMax, X1, Y, Z1, &Pot1, &Flux1);
1062 *Potential = Pot1;
1063 Flux->X = Flux1.X;
1064 Flux->Y = Flux1.Y;
1065 Flux->Z = Flux1.Z;
1066 return 0;
1067 }
1068 if ((fabs(X) < MINDIST) && (fabs(Z - zMax) < MINDIST)) {
1069 // Y line at (0,zMax)
1070 if (DebugISLES) printf("Y line at (0,0,zMax) corner\n");
1071 double X1 = X;
1072 double Z1 = Z;
1073 if ((X >= 0.0) && (Z >= zMax)) {
1074 // point on the element
1075 if (DebugISLES) printf("Case1=>\n");
1076 X1 = X + SHIFT * MINDIST;
1077 Z1 = Z + SHIFT * MINDIST;
1078 } else if ((X <= 0.0) && (Z >= zMax)) {
1079 // field point outside the element
1080 if (DebugISLES) printf("Case2=>\n");
1081 X1 = X - SHIFT * MINDIST;
1082 Z1 = Z + SHIFT * MINDIST;
1083 } else if ((X >= 0.0) && (Z <= zMax)) {
1084 // field point outside the element
1085 if (DebugISLES) printf("Case3=>\n");
1086 X1 = X + SHIFT * MINDIST;
1087 Z1 = Z - SHIFT * MINDIST;
1088 } else if ((X <= 0.0) && (Z <= zMax)) {
1089 // field point outside the element
1090 if (DebugISLES) printf("Case4=>\n");
1091 X1 = X - SHIFT * MINDIST;
1092 Z1 = Z - SHIFT * MINDIST;
1093 }
1094 double Pot1;
1095 Vector3D Flux1;
1096 ExactTriSurf(zMax, X1, Y, Z1, &Pot1, &Flux1);
1097 *Potential = Pot1;
1098 Flux->X = Flux1.X;
1099 Flux->Y = Flux1.Y;
1100 Flux->Z = Flux1.Z;
1101 return 0;
1102 }
1103 // Three edges outside the extent of the real element.
1104 // Plus two edges delineating the virtual right-triangle that complements
1105 // the real one to create a rectangle.
1106 // Special X=1 condition which turns R1+-I2 terms in dire straits
1107 // Similar problem occurs for X=0, where R1+-I1 is NaN but this has been taken
1108 // care of earlier.
1109 // But this particular problem has a remedy - check note on complex tanh-1
1110 // below. There is the problem of dealing with gsl_complex_log_real(1.0)
1111 // though!
1112 if (fabs(X) < MINDIST) {
1113 // edge along X - axis
1114 if (DebugISLES) printf("edge along X-axis\n");
1115 double X1 = X, X2 = X;
1116 if (X >= 0.0) {
1117 // field point on +ve side of YZ plane
1118 if (DebugISLES) printf("Case1=>\n");
1119 X1 = X + SHIFT * MINDIST;
1120 X2 = X - SHIFT * MINDIST;
1121 } else if (X <= 0.0) {
1122 // field point on -ve side of YZ plane
1123 if (DebugISLES) printf("Case2=>\n");
1124 X1 = X - SHIFT * MINDIST;
1125 X2 = X + SHIFT * MINDIST;
1126 }
1127 double Pot1, Pot2;
1128 Vector3D Flux1, Flux2;
1129 ExactTriSurf(zMax, X1, Y, Z, &Pot1, &Flux1);
1130 ExactTriSurf(zMax, X2, Y, Z, &Pot2, &Flux2);
1131 *Potential = 0.5 * (Pot1 + Pot2);
1132 Flux->X = 0.5 * (Flux1.X + Flux2.X);
1133 Flux->Y = 0.5 * (Flux1.Y + Flux2.Y);
1134 Flux->Z = 0.5 * (Flux1.Z + Flux2.Z);
1135 return 0;
1136 }
1137 /* do not erase these blocks as yet
1138 if (fabs(Z) < MINDIST) {
1139 // edge along Z - axis
1140 if (DebugISLES) printf("edge along Z-axis\n");
1141 double Z1 = Z;
1142 if (Z >= 0.0) {
1143 // field point on +ve side of XY plane
1144 if (DebugISLES) printf("Case1=>\n");
1145 Z1 = Z + SHIFT*MINDIST;
1146 } else if (Z <= 0.0) {
1147 // field point on -ve side of XY plane
1148 if (DebugISLES) printf("Case2=>\n");
1149 Z1 = Z - SHIFT*MINDIST;
1150 }
1151 ExactTriSurf(zMax, X, Y, Z1, &Pot1, &Flux1);
1152 *Potential = Pot1;
1153 Flux->X = Flux1.X;
1154 Flux->Y = Flux1.Y;
1155 Flux->Z = Flux1.Z;
1156 return 0;
1157 }
1158 */
1159 if (fabs(X - 1.0) < MINDIST) {
1160 // edge along X=1.0
1161 if (DebugISLES) printf("edge along X = 1.\n");
1162 double X1 = X, X2 = X;
1163 if (X <= 1.0) {
1164 // field point on +ve side of YZ plane
1165 if (DebugISLES) printf("Case1=>\n");
1166 X1 = X - SHIFT * MINDIST;
1167 X2 = X + SHIFT * MINDIST;
1168 } else if (X >= 1.0) {
1169 // field point on -ve side of YZ plane
1170 if (DebugISLES) printf("Case2=>\n");
1171 X1 = X + SHIFT * MINDIST;
1172 X2 = X + SHIFT * MINDIST;
1173 }
1174 double Pot1, Pot2;
1175 Vector3D Flux1, Flux2;
1176 ExactTriSurf(zMax, X1, Y, Z, &Pot1, &Flux1);
1177 ExactTriSurf(zMax, X2, Y, Z, &Pot2, &Flux2);
1178 *Potential = 0.5 * (Pot1 + Pot2);
1179 Flux->X = 0.5 * (Flux1.X + Flux2.X);
1180 Flux->Y = 0.5 * (Flux1.Y + Flux2.Y);
1181 Flux->Z = 0.5 * (Flux1.Z + Flux2.Z);
1182 return 0;
1183 }
1184 if (fabs(Z - zMax) < MINDIST) {
1185 // edge along Z=zMax
1186 if (DebugISLES) printf("edge along Z = zMax\n");
1187 double Z1 = Z, Z2 = Z;
1188 if (Z >= zMax) {
1189 // field point on +ve side of XY plane
1190 if (DebugISLES) printf("Case1=>\n");
1191 Z1 = Z + SHIFT * MINDIST;
1192 Z2 = Z - SHIFT * MINDIST;
1193 } else if (Z <= zMax) {
1194 // field point on -ve side of XY plane
1195 if (DebugISLES) printf("Case2=>\n");
1196 Z1 = Z - SHIFT * MINDIST;
1197 Z2 = Z + SHIFT * MINDIST;
1198 }
1199 double Pot1, Pot2;
1200 Vector3D Flux1, Flux2;
1201 ExactTriSurf(zMax, X, Y, Z1, &Pot1, &Flux1);
1202 ExactTriSurf(zMax, X, Y, Z2, &Pot2, &Flux2);
1203 *Potential = 0.5 * (Pot1 + Pot2);
1204 Flux->X = 0.5 * (Flux1.X + Flux2.X);
1205 Flux->Y = 0.5 * (Flux1.Y + Flux2.Y);
1206 Flux->Z = 0.5 * (Flux1.Z + Flux2.Z);
1207 return 0;
1208 }
1209 if (fabs(Z - Zhyp) < MINDIST) {
1210 // edge along the hypotenuse
1211 if (DebugISLES) printf("edge along Hypotenuse\n");
1212 double X1 = X, Z1 = Z;
1213 if (Z <= Zhyp) {
1214 // towards element origin
1215 if (DebugISLES) printf("Case1=>\n");
1216 X1 = X - SHIFT * MINDIST;
1217 Z1 = Z - SHIFT * MINDIST;
1218 } else if (Z >= Zhyp) {
1219 // going further away from the element origin
1220 if (DebugISLES) printf("Case2=>\n");
1221 X1 = X + SHIFT * MINDIST;
1222 Z1 = Z + SHIFT * MINDIST;
1223 }
1224 double Pot1;
1225 Vector3D Flux1;
1226 ExactTriSurf(zMax, X1, Y, Z1, &Pot1, &Flux1);
1227 *Potential = Pot1;
1228 Flux->X = Flux1.X;
1229 Flux->Y = Flux1.Y;
1230 Flux->Z = Flux1.Z;
1231 return 0;
1232 }
1233
1234 // Related to logarithmic terms (LTerm1 and LTerm2)
1235 if ((fabs(G) <= MINDIST) && (modY <= MINDIST)) {
1236 ApproxFlag = 10;
1237 ++FailureCntr;
1238 --ExactCntr;
1239 fprintf(
1240 fIsles,
1241 "(fabs(G) <= MINDIST) && (modY <= MINDIST) ... approximating: %d.\n",
1242 ApproxFlag);
1243 return (ApproxTriSurf(zMax, X, Y, Z, XNSegApprox, ZNSegApprox, Potential,
1244 Flux));
1245 }
1246 /* do not erase these blocks as yet
1247 if ((fabs(X) <= MINDIST) && (modY <= MINDIST)) {
1248 // denominator zero for LP1 and LM1
1249 ApproxFlag = 11; ++FailureCntr; --ExactCntr;
1250 fprintf(fIsles, "denominator zero for LP1 and LM1 ... approximating: %d.\n",
1251 ApproxFlag);
1252 return (ApproxTriSurf(zMax, X, Y, Z, XNSegApprox, ZNSegApprox, Potential,
1253 Flux));
1254 }
1255 if ((fabs(H1 + D12 * G) <= MINDIST2) && (modY <= MINDIST)) {
1256 // numerator zero for LP1 and LM1
1257 ApproxFlag = 12; ++FailureCntr; --ExactCntr;
1258 fprintf(fIsles, "numerator zero for LP1 and LM1 ... approximating: %d.\n",
1259 ApproxFlag);
1260 return(ApproxTriSurf(zMax, X, Y, Z, XNSegApprox, ZNSegApprox, Potential,
1261 Flux));
1262 }
1263 if ((fabs(1.0 - X) <= MINDIST) && (modY <= MINDIST) ) {
1264 // denominator zero for LP2 and LM2
1265 ApproxFlag = 13; ++FailureCntr; --ExactCntr;
1266 fprintf(fIsles, "denominator zero for LP2 and LM2 ... approximating: %d.\n",
1267 ApproxFlag);
1268 return(ApproxTriSurf(zMax, X, Y, Z, XNSegApprox, ZNSegApprox, Potential,
1269 Flux));
1270 }
1271 if ((fabs(H2 + D21 * G) <= MINDIST2) && (modY <= MINDIST)) {
1272 // numerator zero for LP2 and LM2
1273 ApproxFlag = 14; ++FailureCntr; --ExactCntr;
1274 fprintf(fIsles, "numerator zero for LP2 and LM2 ... approximating: %d.\n",
1275 ApproxFlag);
1276 return (ApproxTriSurf(zMax, X, Y, Z, XNSegApprox, ZNSegApprox, Potential,
1277 Flux));
1278 }
1279 */
1280
1281 // Related to complex inverse tan hyperbolic terms - TanhTerm1 and TanhTerm2
1282 /* do not erase these blocks as yet
1283 if (D11 * fabs(Z) <= MINDIST2) {
1284 ApproxFlag = 15; ++FailureCntr; --ExactCntr;
1285 fprintf(fIsles, "D11*fabs(Z) zero for TanhTerms ... approximating: %d.\n",
1286 ApproxFlag);
1287 return(ApproxTriSurf(zMax, X, Y, Z, XNSegApprox, ZNSegApprox, Potential, Flux));
1288 }
1289 if (D21 * fabs(Z) <= MINDIST2) {
1290 ApproxFlag = 16; ++FailureCntr; --ExactCntr;
1291 fprintf(fIsles, "D21*fabs(Z) zero for TanhTerms ... approximating: %d.\n",
1292 ApproxFlag);
1293 return(ApproxTriSurf(zMax, X, Y, Z, XNSegApprox, ZNSegApprox, Potential, Flux));
1294 }
1295 */
1296
1297 // Exact computations begin here, at last!
1298 // DTerm1 and DTerm2
1299 double DblTmp1 = (Hypot * D12 - E1) / (Hypot * D21 - E2);
1300 if (DebugISLES) {
1301 printf("DblTmp1: %.16lg\n", DblTmp1);
1302 fflush(stdout);
1303 }
1304 if (DblTmp1 < MINDIST2) DblTmp1 = MINDIST2;
1305 double DTerm1 = log(DblTmp1);
1306 if (isnan(DTerm1) || isinf(DTerm1)) {
1307 ApproxFlag = 18;
1308 ++FailureCntr;
1309 --ExactCntr;
1310 fprintf(fIsles, "DTerm1 nan or inf ... approximating: %d.\n", ApproxFlag);
1311 return (ApproxTriSurf(zMax, X, Y, Z, XNSegApprox, ZNSegApprox, Potential,
1312 Flux));
1313 }
1314 DTerm1 /= Hypot;
1315
1316 double DblTmp2 = (D11 - X) / (D21 - X + 1.0);
1317 if (DebugISLES) {
1318 printf("DblTmp2: %.16lg\n", DblTmp2);
1319 fflush(stdout);
1320 }
1321 if (DblTmp2 < MINDIST2) DblTmp2 = MINDIST2;
1322 double DTerm2 = log(DblTmp2);
1323 if (isnan(DTerm2) || isinf(DTerm2)) {
1324 ApproxFlag = 18;
1325 ++FailureCntr;
1326 --ExactCntr;
1327 fprintf(fIsles, "DTerm2 nan or inf ... approximating: %d.\n", ApproxFlag);
1328 return (ApproxTriSurf(zMax, X, Y, Z, XNSegApprox, ZNSegApprox, Potential,
1329 Flux));
1330 }
1331
1332 if (DebugISLES) {
1333 printf("DTerm1: %.16lg, DTerm2: %.16lg\n", DTerm1, DTerm2);
1334 fflush(stdout);
1335 }
1336
1337 // Logarithmic terms
1338 double deltaLog = 0.;
1339 double deltaArg = 0.;
1340 if (modY > MINDIST) {
1341 const double s1 = 1. / (X * X + modY * modY);
1342 double Re1 = ((-X) * (H1 + G * D12) + modY * modY * (E1 - zMax * D12)) * s1;
1343 double Im1 = ((-X) * modY * (E1 - zMax * D12) - (H1 + G * D12) * modY) * s1;
1344 const double s2 = 1. / ((1. - X) * (1. - X) + modY * modY);
1345 double Re2 = ((1. - X) * (H2 + G * D21) + modY * modY * (E2 - zMax * D21)) * s2;
1346 double Im2 = ((1. - X) * modY * (E2 - zMax * D21) - (H2 + G * D21) * modY) * s2;
1347 deltaLog = 0.5 * log((Re1 * Re1 + Im1 * Im1) / (Re2 * Re2 + Im2 * Im2));
1348 deltaArg = atan2(Im1, Re1) - atan2(Im2, Re2);
1349 } else {
1350 double f1 = fabs((H1 + G * D12) / (-X));
1351 if (f1 < MINDIST) f1 = MINDIST;
1352 double f2 = fabs((H2 + G * D21) / (1. - X));
1353 if (f2 < MINDIST) f2 = MINDIST;
1354 deltaLog = log(f1 / f2);
1355 }
1356 const double c1 = 2. / (G * G + zMax * zMax * modY * modY);
1357 double LTerm1 = c1 * (G * deltaLog - zMax * modY * deltaArg);
1358 double LTerm2 = c1 * (G * deltaArg + zMax * modY * deltaLog);
1359 if (DebugISLES) {
1360 printf("LTerm1: %.16lg, LTerm2: %.16lg\n", LTerm1, LTerm2);
1361 fflush(stdout);
1362 }
1363
1364 // Computations for estimating TanhTerm1, TanhTerm2
1365 // Possible singularities - D11, D21 corners and Z=0 line / surface
1366 // Corners have been taken care of, as well as Z=0 (latter, through S1)
1367 // Incorrectly implemented - CHECK!!! CORRECTED!
1368 // If the argument of atanh is real and greater than 1.0, it is
1369 // perfectly computable. The value returned is a complex number,
1370 // however. This computation has to be carried out by invoking
1371 // gsl_complex_arctanh_real(double z).
1372 // The branch cut needs to be enforced only if the
1373 // argument is real and is equal to 1.0. Then the returned value
1374 // is undefined. But this happens only if Y=0 besides X=1.0!
1375 // Similar implementation issues are also there for X=0.
1376 // Coincides with X=1 since I2 turns out to be zero in such
1377 // cases. CmTmp3.dat[0] can be >= 1.0 in a variety of situations
1378 // since in the denominator D21, sqrt( (X-1)^2 + Y^2 + Z^2 ) is
1379 // sqrt(Y^2+Z^2) for X=1. This is multiplied by |Z| while on the
1380 // numerator we have Y^2+Z^2.
1381
1382 double TanhTerm1 = 0.;
1383 double TanhTerm2 = 0.;
1384 if (abs(S1) > 0) {
1385 if (modY < MINDIST) {
1386 const double f1 = R1 / (D11 * fabs(Z));
1387 const double f2 = R1 / (D21 * fabs(Z));
1388 TanhTerm1 = log1p(f1) - log1p(-f1) - log1p(f2) + log1p(-f2);
1389 } else {
1390 TanhTerm1 = 1.;
1391 if (fabs(R1 - D11 * fabs(Z)) > MINDIST2 || true) {
1392 const double rp = D11 * fabs(Z) + R1;
1393 const double rm = D11 * fabs(Z) - R1;
1394 TanhTerm1 *= (I1 * I1 + rp * rp) / (I1 * I1 + rm * rm);
1395 }
1396 if (fabs(R1 - D21 * fabs(Z)) > MINDIST2 || true) {
1397 const double rp = D21 * fabs(Z) + R1;
1398 const double rm = D21 * fabs(Z) - R1;
1399 TanhTerm1 *= (I2 * I2 + rm * rm) / (I2 * I2 + rp * rp);
1400 }
1401 TanhTerm1 = 0.5 * log(TanhTerm1);
1402 }
1403
1404 if (fabs(I1) > MINDIST2) {
1405 double a = D11 * D11 * Z * Z;
1406 double b = R1 * R1 + I1 * I1;
1407 double tmp = atan(2 * I1 * D11 * fabs(Z) / (a - b));
1408 if (b > a) {
1409 if (X > 0.) {
1410 tmp += ST_PI;
1411 } else {
1412 tmp -= ST_PI;
1413 }
1414 }
1415 TanhTerm2 += tmp;
1416 }
1417 if (fabs(I2) > MINDIST2) {
1418 double a = D21 * D21 * Z * Z;
1419 double b = R1 * R1 + I2 * I2;
1420 double tmp = atan(2 * I2 * D21 * fabs(Z) / (a - b));
1421 if (b > a) {
1422 if (X > 1.) {
1423 tmp += ST_PI;
1424 } else {
1425 tmp -= ST_PI;
1426 }
1427 }
1428 TanhTerm2 -= tmp;
1429 }
1430 }
1431
1432 if (DebugISLES) {
1433 printf("TanhTerm1: %.16lg, TanhTerm2: %.16lg\n", TanhTerm1, TanhTerm2);
1434 fflush(stdout);
1435 }
1436
1437 double Pot = (zMax * Y * Y - X * G) * LTerm1 -
1438 (zMax * X + G) * modY * LTerm2 +
1439 (S1 * X * TanhTerm1) + S1 * modY * TanhTerm2 +
1440 2. * (G * DTerm1 - Z * DTerm2);
1441 Pot *= 0.5;
1442 double Fx = G * LTerm1 + zMax * modY * LTerm2 - S1 * TanhTerm1 - 2.0 * zMax * DTerm1;
1443 Fx *= 0.5;
1444 double Fy = -zMax * Y * LTerm1 + G * SY * LTerm2 - S1 * SY * TanhTerm2;
1445 Fy *= 0.5;
1446
1447 double Fz = DTerm2 - DTerm1;
1448 if (DebugISLES) {
1449 printf("Pot: %.16lg, Fx: %.16lg, Fy: %.16lg, Fz: %.16lg\n", Pot, Fx, Fy,
1450 Fz);
1451 fflush(stdout);
1452 }
1453
1454 // Final adjustments
1455 // Constants (?) of integration - Carry out only one of the options
1456 // Depends criticially on > or >=; similarly < or <=
1457 // Needs further investigation
1458 // As far as Fy is concerned, conditions 1 & 3, and 2 & 4 can be combined.
1459 // So, instead of the triangular area, it seems, that the rectangular bound
1460 // is more important. In such an event, Zhyp will be redundant.
1461 if ((X >= 0.0) && (X <= 1.0)) {
1462 // Possibility of point within element bounds
1463 int ConstAdd = 0;
1464 if ((Z >= 0.0) && (Z <= Zhyp)) {
1465 // within the element
1466 ConstAdd = 1;
1467 Pot -= modY * ST_PI;
1468 if (fabs(Y) < MINDIST) // on the element surface
1469 Fy = ST_PI;
1470 else if (Y > 0.0) // above or below the element surface
1471 Fy += ST_PI;
1472 else if (Y < 0.0)
1473 Fy -= ST_PI;
1474 } else if ((Z <= 0.0) && (Z >= -Zhyp)) {
1475 // within -ve element
1476 ConstAdd = 2;
1477 Pot += modY * ST_PI;
1478 if (fabs(Y) < MINDIST)
1479 Fy = 0.0;
1480 else if (Y > 0.0)
1481 Fy -= ST_PI;
1482 else if (Y < 0.0)
1483 Fy += ST_PI;
1484 } else if (((Z > Zhyp) && (Z <= zMax))) {
1485 // within +rect bounds
1486 // merge with +ve shadow?
1487 ConstAdd = 3;
1488 Pot -= modY * ST_PI;
1489 if (fabs(Y) < MINDIST)
1490 Fy = 0.0;
1491 else if (Y > 0.0)
1492 Fy += ST_PI;
1493 else if (Y < 0.0)
1494 Fy -= ST_PI;
1495 } else if ((Z < -Zhyp) && (Z >= -zMax)) {
1496 // within -rect bounds
1497 // merge with -ve shadow?
1498 ConstAdd = 4;
1499 Pot += modY * ST_PI;
1500 if (fabs(Y) < MINDIST)
1501 Fy = 0.0;
1502 else if (Y > 0.0)
1503 Fy -= ST_PI;
1504 else if (Y < 0.0)
1505 Fy += ST_PI;
1506 } else if (Z > zMax) {
1507 // +ve shadow of the triangle - WHY?
1508 ConstAdd = 5;
1509 Pot -= modY * ST_PI;
1510 if (fabs(Y) < MINDIST)
1511 Fy = 0.0;
1512 else if (Y > 0.0)
1513 Fy += ST_PI;
1514 else if (Y < 0.0)
1515 Fy -= ST_PI;
1516 } else if (Z < -zMax) {
1517 // -ve shadow of the triangle - WHY?
1518 ConstAdd = 6;
1519 Pot += modY * ST_PI;
1520 if (fabs(Y) < MINDIST)
1521 Fy = 0.0;
1522 else if (Y > 0.0)
1523 Fy -= ST_PI;
1524 else if (Y < 0.0)
1525 Fy += ST_PI;
1526 }
1527
1528 /*
1529 if ((fabs(X - 1.0) < MINDIST) && (fabs(Z-zMax) < MINDIST)) {
1530 // (1,zMax) corner
1531 if (Y > 0.0) Fy += 0.5*ST_PI;
1532 if (Y < 0.0) Fy -= 0.5*ST_PI;
1533 } else if ((fabs(X - 1.0) < MINDIST) && (fabs(Z + zMax) < MINDIST)) {
1534 // (1,-zMax) corner
1535 if (Y > 0.0) Fy -= 0.5*ST_PI;
1536 if (Y < 0.0) Fy += 0.5*ST_PI;
1537 } else if ((fabs(X) < MINDIST) && (Z < 0.0)) {
1538 // X = 0, Z < 0 - off ele
1539 // along -ve Z axis
1540 if (Y > 0.0) Fy -= 0.5*ST_PI;
1541 if (Y < 0.0) Fy += 0.5*ST_PI;
1542 } else if ((fabs(X) < MINDIST) && (Z > 0.0)) {
1543 // X = 0, Z > 0 - on & off ele
1544 // along +ve Z axis
1545 if (Y > 0.0) Fy += 0.5*ST_PI;
1546 if (Y < 0.0) Fy -= 0.5*ST_PI;
1547 } else if (Z > 0.0) {
1548 // within rectangular bounds - changed for INPC
1549 if (Y > 0.0) Fy += ST_PI;
1550 if (Y < 0.0) Fy -= ST_PI;
1551 } else if (Z < 0.0) {
1552 // within -ve rectangular bounds - changed for INPC
1553 if (Y > 0.0) Fy -= ST_PI;
1554 if (Y < 0.0) Fy += ST_PI;
1555 }
1556 */
1557 // two other conditions, one for Z > zMax and another for Z < -zMax expected
1558
1559 if (DebugISLES) {
1560 printf("After constant addition %d\n", ConstAdd);
1561 printf("Pot: %.16lg, Fx: %.16lg, Fy: %.16lg, Fz: %.16lg\n", Pot, Fx, Fy,
1562 Fz);
1563 fflush(stdout);
1564 }
1565 } // point within element bounds
1566
1567 // Error situations handled before returning the potential value
1568 if ((Pot < 0) || isnan(Pot) || isinf(Pot)) {
1569 fprintf(fIsles, "\n---Approximation in ExactTriSurf---\n");
1570 fprintf(fIsles, "Problem with potential ... approximating!\n");
1571 fprintf(fIsles, "zMax: %.16lg, X: %.16lg, Y: %.16lg, Z: %.16lg\n", zMax, X, Y,
1572 Z);
1573 fprintf(fIsles, "D11: %.16lg, D21: %.16lg, D12: %.16lg, Hypot: %.16lg\n", D11,
1574 D21, D12, Hypot);
1575 fprintf(fIsles, "modY: %.16lg, G: %.16lg\n", modY, G);
1576 fprintf(fIsles, "E1: %.16lg, E2: %.16lg, H1: %.16lg, H2: %.16lg\n", E1, E2,
1577 H1, H2);
1578 fprintf(fIsles, "S1: %d, R1: %.16lg, I1: %.16lg, I2: %.16lg\n", S1, R1, I1,
1579 I2);
1580 fprintf(fIsles, "Pot: %.16lg\n", Pot);
1581 ApproxFlag = 23;
1582 ++FailureCntr;
1583 --ExactCntr;
1584 return (ApproxTriSurf(zMax, X, Y, Z, XNSegApprox, ZNSegApprox, Potential, Flux));
1585 }
1586
1587 if (isnan(Fx) || isinf(Fx)) {
1588 fprintf(fIsles, "\n---Approximation in ExactTriSurf---\n");
1589 fprintf(fIsles, "Problem with Fx ... approximating!\n");
1590 fprintf(fIsles, "zMax: %.16lg, X: %.16lg, Y: %.16lg, Z: %.16lg\n", zMax, X, Y,
1591 Z);
1592 fprintf(fIsles, "D11: %.16lg, D21: %.16lg, D12: %.16lg, Hypot: %.16lg\n", D11,
1593 D21, D12, Hypot);
1594 fprintf(fIsles, "modY: %.16lg, G: %.16lg\n", modY, G);
1595 fprintf(fIsles, "E1: %.16lg, E2: %.16lg, H1: %.16lg, H2: %.16lg\n", E1, E2,
1596 H1, H2);
1597 fprintf(fIsles, "S1: %d, R1: %.16lg, I1: %.16lg, I2: %.16lg\n", S1, R1, I1,
1598 I2);
1599 fprintf(fIsles, "Fx: %.16lg\n", Fx);
1600 ApproxFlag = 24;
1601 ++FailureCntr;
1602 --ExactCntr;
1603 return (ApproxTriSurf(zMax, X, Y, Z, XNSegApprox, ZNSegApprox, Potential, Flux));
1604 }
1605
1606 if (isnan(Fy) || isinf(Fy)) {
1607 fprintf(fIsles, "\n---Approximation in ExactTriSurf---\n");
1608 fprintf(fIsles, "Problem with Fy ... approximating!\n");
1609 fprintf(fIsles, "zMax: %.16lg, X: %.16lg, Y: %.16lg, Z: %.16lg\n", zMax, X, Y,
1610 Z);
1611 fprintf(fIsles, "D11: %.16lg, D21: %.16lg, D12: %.16lg, Hypot: %.16lg\n", D11,
1612 D21, D12, Hypot);
1613 fprintf(fIsles, "modY: %.16lg, G: %.16lg\n", modY, G);
1614 fprintf(fIsles, "E1: %.16lg, E2: %.16lg, H1: %.16lg, H2: %.16lg\n", E1, E2,
1615 H1, H2);
1616 fprintf(fIsles, "S1: %d, R1: %.16lg, I1: %.16lg, I2: %.16lg\n", S1, R1, I1,
1617 I2);
1618 fprintf(fIsles, "Fy: %.16lg\n", Fy);
1619 ApproxFlag = 25;
1620 ++FailureCntr;
1621 --ExactCntr;
1622 return (ApproxTriSurf(zMax, X, Y, Z, XNSegApprox, ZNSegApprox, Potential, Flux));
1623 }
1624
1625 if (isnan(Fz) || isinf(Fz)) {
1626 fprintf(fIsles, "\n---Approximation in ExactTriSurf---\n");
1627 fprintf(fIsles, "Problem with Fz ... approximating!\n");
1628 fprintf(fIsles, "zMax: %.16lg, X: %.16lg, Y: %.16lg, Z: %.16lg\n", zMax, X, Y,
1629 Z);
1630 fprintf(fIsles, "D11: %.16lg, D21: %.16lg, D12: %.16lg, Hypot: %.16lg\n", D11,
1631 D21, D12, Hypot);
1632 fprintf(fIsles, "modY: %.16lg\n", modY);
1633 fprintf(fIsles, "E1: %.16lg, E2: %.16lg\n", E1, E2);
1634 fprintf(fIsles, "Fz: %.16lg\n", Fz);
1635 ApproxFlag = 26;
1636 ++FailureCntr;
1637 --ExactCntr;
1638 return (ApproxTriSurf(zMax, X, Y, Z, XNSegApprox, ZNSegApprox, Potential, Flux));
1639 }
1640
1641 // Return the computed value of potential and flux
1642 *Potential = Pot;
1643 Flux->X = Fx;
1644 Flux->Y = Fy;
1645 Flux->Z = Fz;
1646
1647 if (DebugISLES) printf("Going out of ExactTriSurf ...\n\n");
1648 return 0;
1649} // ExactTriSurf ends
ApproxTriSurf
int ApproxTriSurf(double zMax, double X, double Y, double Z, int nbxseg, int nbzseg, double *Potential, Vector3D *Flux)
Definition: Isles.c:1660
ExactTriSurf
int ExactTriSurf(double zMax, double X, double Y, double Z, double *Potential, Vector3D *Flux)
Definition: Isles.c:784

Referenced by ExactTriSurf(), TriFlux(), TriPF(), TriPot(), and TriPrimPF().

◆ ImprovedFX_W()

double ImprovedFX_W ( double  rW,
double  lW,
double  X,
double  Y,
double  Z 
)

Definition at line 1946 of file Isles.c.

1946 {
1947 if (DebugISLES) printf("In ImprovedFX_W ...\n");
1948
1949 double dist = sqrt(X * X + Y * Y + Z * Z);
1950 if (dist < MINDIST) {
1951 // distance less than MINDIST
1952 return (0.0);
1953 } else if ((fabs(X) < MINDIST) && (fabs(Y) < MINDIST)) {
1954 // point on the axis of the wire element
1955 return (0.0);
1956 } else {
1957 // point far away from the wire element
1958 double C = (Z) - (lW / 2.0);
1959 double D = (Z) + (lW / 2.0);
1960 double A = sqrt(X * X + Y * Y + C * C);
1961 double B = sqrt(X * X + Y * Y + D * D);
1962
1963 double tmp1 = 2.0 * ST_PI * rW;
1964 double tmp2 = (X / (A * (B - D))) - ((A - C) * X / (B * (B - D) * (B - D)));
1965 double tmp3 = (B - D) / (A - C);
1966 return (-1.0 * tmp1 * tmp2 * tmp3);
1967 } // dist ... if ... else if ... else
1968} // ImprovedFX_W ends

◆ ImprovedFY_W()

double ImprovedFY_W ( double  rW,
double  lW,
double  X,
double  Y,
double  Z 
)

Definition at line 1970 of file Isles.c.

1970 {
1971 if (DebugISLES) printf("In ImprovedFY_W ...\n");
1972
1973 double dist = sqrt(X * X + Y * Y + Z * Z);
1974 if (dist < MINDIST) {
1975 // distance less than MINDIST
1976 return (0.0);
1977 } else if ((fabs(X) < MINDIST) && (fabs(Y) < MINDIST)) {
1978 // point on the axis of the wire element
1979 return (0.0);
1980 } else {
1981 // point far away from the wire element
1982 double C = (Z) - (lW / 2.0);
1983 double D = (Z) + (lW / 2.0);
1984 double A = sqrt(X * X + Y * Y + C * C);
1985 double B = sqrt(X * X + Y * Y + D * D);
1986
1987 double tmp1 = 2.0 * ST_PI * rW;
1988 double tmp2 = (Y / (A * (B - D))) - ((A - C) * Y / (B * (B - D) * (B - D)));
1989 double tmp3 = (B - D) / (A - C);
1990 return (-1.0 * tmp1 * tmp2 * tmp3);
1991 } // dist ... if ... else if ... else
1992} // ImprovedFY_W ends

◆ ImprovedFZ_W()

double ImprovedFZ_W ( double  rW,
double  lW,
double  X,
double  Y,
double  Z 
)

Definition at line 1994 of file Isles.c.

1994 {
1995 if (DebugISLES) printf("In ImprovedFZ_W ...\n");
1996
1997 double dist = sqrt(X * X + Y * Y + Z * Z);
1998 if (dist < MINDIST) {
1999 // distance less than MINDIST
2000 return (0.0);
2001 } else if ((fabs(X) < MINDIST) && (fabs(Y) < MINDIST)) {
2002 // point on the axis of the wire element
2003 double C = Z - (lW / 2.0);
2004 double D = Z + (lW / 2.0);
2005 double A = sqrt(X * X + Y * Y + C * C);
2006 double B = sqrt(X * X + Y * Y + D * D);
2007
2008 double tmp1 = 2.0 * ST_PI * rW;
2009 double tmp2 = (1.0 / B) - (1.0 / A);
2010 return (-1.0 * tmp1 * tmp2);
2011 } else {
2012 // point far away from the wire element
2013 double C = Z - (lW / 2.0);
2014 double D = Z + (lW / 2.0);
2015 double A = sqrt(X * X + Y * Y + C * C);
2016 double B = sqrt(X * X + Y * Y + D * D);
2017
2018 double tmp1 = 2.0 * ST_PI * rW;
2019 double tmp2 = (1.0 / B) - (1.0 / A);
2020 return (-1.0 * tmp1 * tmp2);
2021 } // dist ... if ... else if ... else
2022} // ImprovedFZ_W ends

◆ ImprovedP_W()

double ImprovedP_W ( double  rW,
double  lW,
double  X,
double  Y,
double  Z 
)

Definition at line 1930 of file Isles.c.

1930 {
1931 // Improved potential at far away points:
1932 // wire2.m applicable for thin wires
1933 if (DebugISLES) printf("In ImprovedP_W ...\n");
1934
1935 double dz = 0.5 * lW; // half length of the wire segment
1936 double dtmp1 = (X * X + Y * Y + (Z - dz) * (Z - dz));
1937 dtmp1 = sqrt(dtmp1);
1938 dtmp1 = dtmp1 - (Z - dz);
1939 double dtmp2 = (X * X + Y * Y + (Z + dz) * (Z + dz));
1940 dtmp2 = sqrt(dtmp2);
1941 dtmp2 = dtmp2 - (Z + dz);
1942 double dtmp3 = log(dtmp1 / dtmp2);
1943 return (2.0 * ST_PI * rW * dtmp3);
1944} // ImprovedP_W ends

◆ ImprovedWire()

int ImprovedWire ( double  rW,
double  lW,
double  X,
double  Y,
double  Z,
double *  potential,
Vector3D *  Flux 
)

Definition at line 2024 of file Isles.c.

2025 {
2026 if (DebugISLES) printf("In ImprovedWire ...\n");
2027
2028 double dz = 0.5 * lW; // half length of the wire segment
2029
2030 double dtmp1 = (X * X + Y * Y + (Z - dz) * (Z - dz));
2031 dtmp1 = sqrt(dtmp1);
2032 dtmp1 = dtmp1 - (Z - dz);
2033 double dtmp2 = (X * X + Y * Y + (Z + dz) * (Z + dz));
2034 dtmp2 = sqrt(dtmp2);
2035 dtmp2 = dtmp2 - (Z + dz);
2036 double dtmp3 = log(dtmp1 / dtmp2);
2037 *potential = 2.0 * ST_PI * rW * dtmp3;
2038
2039 double dist = sqrt(X * X + Y * Y + Z * Z);
2040 double Fx = 0.0, Fy = 0.0, Fz = 0.0;
2041
2042 if (dist < MINDIST) {
2043 // distance less than MINDIST from centroid
2044 Fx = 0.0;
2045 Fy = 0.0;
2046 Fz = 0.0;
2047 } else if ((fabs(X) < MINDIST) && (fabs(Y) < MINDIST)) {
2048 // point on the axis of the wire element
2049 Fx = 0.0;
2050 Fy = 0.0;
2051
2052 double C = Z - (lW / 2.0);
2053 double D = Z + (lW / 2.0);
2054 double A = sqrt(X * X + Y * Y + C * C);
2055 double B = sqrt(X * X + Y * Y + D * D);
2056
2057 double tmp1 = 2.0 * ST_PI * rW;
2058 double tmp2 = (1.0 / B) - (1.0 / A);
2059 Fz = -1.0 * tmp1 * tmp2;
2060 } else {
2061 // point far away from the wire element
2062 double C = (Z) - (lW / 2.0);
2063 double D = (Z) + (lW / 2.0);
2064 double A = sqrt(X * X + Y * Y + C * C);
2065 double B = sqrt(X * X + Y * Y + D * D);
2066
2067 double tmp1 = 2.0 * ST_PI * rW;
2068 double tmp2 = (X / (A * (B - D))) - ((A - C) * X / (B * (B - D) * (B - D)));
2069 double tmp3 = (B - D) / (A - C);
2070 Fx = -1.0 * tmp1 * tmp2 * tmp3;
2071
2072 tmp2 = (Y / (A * (B - D))) - ((A - C) * Y / (B * (B - D) * (B - D)));
2073 Fy = -1.0 * tmp1 * tmp2 * tmp3;
2074
2075 tmp2 = (1.0 / B) - (1.0 / A);
2076 Fz = -1.0 * tmp1 * tmp2;
2077 } // dist ... if ... else if ... else
2078
2079 Flux->X = Fx;
2080 Flux->Y = Fy;
2081 Flux->Z = Fz;
2082
2083 return 0;
2084}

◆ LineKnChPF()

double LineKnChPF ( Point3D  LineStart,
Point3D  LineStop,
double  radius,
Point3D  FieldPt,
Vector3D *  globalF 
)

Definition at line 2199 of file Isles.c.

2200 {
2201 double xvert[2], yvert[2], zvert[2];
2202 xvert[0] = LineStart.X;
2203 xvert[1] = LineStop.X;
2204 yvert[0] = LineStart.Y;
2205 yvert[1] = LineStop.Y;
2206 zvert[0] = LineStart.Z;
2207 zvert[1] = LineStop.Z;
2208
2209 double xfld = FieldPt.X;
2210 double yfld = FieldPt.Y;
2211 double zfld = FieldPt.Z;
2212
2213 double xorigin = 0.5 * (xvert[1] + xvert[0]);
2214 double yorigin = 0.5 * (yvert[1] + yvert[0]);
2215 double zorigin = 0.5 * (zvert[1] + zvert[0]);
2216 double LZ = GetDistancePoint3D(&LineStop, &LineStart);
2217
2218 // Create a local coordinate system for which the z axis is along the length
2219 // of the wire
2220 // FieldPt needs to be transformed to localPt
2221
2222 // Find direction cosines of the line charge to initiate the transformation.
2223 DirnCosn3D DirCos;
2224
2225 // Copied from the DiscretizeWire function of neBEM/src/PreProcess/ReTrim.c
2226 // Direction cosines along the wire - note difference from surface primitives!
2227 // The direction along the wire is considered to be the z axis of the LCS.
2228 // So, let us fix that axial vector first
2229 DirCos.ZUnit.X = (xvert[1] - xvert[0]) / LZ; // useful
2230 DirCos.ZUnit.Y = (yvert[1] - yvert[0]) / LZ;
2231 DirCos.ZUnit.Z = (zvert[1] - zvert[0]) / LZ; // useful
2232 // Now direction cosines for the X and Y axes.
2233 {
2234 Vector3D XUnit, YUnit, ZUnit;
2235 ZUnit.X = DirCos.ZUnit.X;
2236 ZUnit.Y = DirCos.ZUnit.Y;
2237 ZUnit.Z = DirCos.ZUnit.Z;
2238
2239 if (fabs(ZUnit.X) >= fabs(ZUnit.Z) && fabs(ZUnit.Y) >= fabs(ZUnit.Z)) {
2240 XUnit.X = -ZUnit.Y;
2241 XUnit.Y = ZUnit.X;
2242 XUnit.Z = 0.0;
2243 } else if (fabs(ZUnit.X) >= fabs(ZUnit.Y) &&
2244 fabs(ZUnit.Z) >= fabs(ZUnit.Y)) {
2245 XUnit.X = -ZUnit.Z;
2246 XUnit.Y = 0.0;
2247 XUnit.Z = ZUnit.X;
2248 } else {
2249 XUnit.X = 0.0;
2250 XUnit.Y = ZUnit.Z;
2251 XUnit.Z = -ZUnit.Y;
2252 }
2253 XUnit = UnitVector3D(&XUnit);
2254
2255 // y-Axis: vectorial product of axes 1 and 2.
2256 YUnit = Vector3DCrossProduct(&ZUnit, &XUnit);
2257 YUnit = UnitVector3D(&YUnit);
2258 // end of replacement
2259
2260 DirCos.XUnit.X = XUnit.X;
2261 DirCos.XUnit.Y = XUnit.Y;
2262 DirCos.XUnit.Z = XUnit.Z;
2263 DirCos.YUnit.X = YUnit.X;
2264 DirCos.YUnit.Y = YUnit.Y;
2265 DirCos.YUnit.Z = YUnit.Z;
2266 } // X and Y direction cosines computed
2267
2268 Point3D localPt;
2269 // Through InitialVector[], field point gets translated to ECS origin.
2270 // Axes direction are, however, still global which when rotated to ECS
2271 // system, yields FinalVector[].
2272 { // Rotate point3D from global to local system to get localPt.
2273 double InitialVector[4];
2274 double TransformationMatrix[4][4] = {{0.0, 0.0, 0.0, 0.0},
2275 {0.0, 0.0, 0.0, 0.0},
2276 {0.0, 0.0, 0.0, 0.0},
2277 {0.0, 0.0, 0.0, 1.0}};
2278 double FinalVector[4];
2279
2280 InitialVector[0] = xfld - xorigin;
2281 InitialVector[1] = yfld - yorigin;
2282 InitialVector[2] = zfld - zorigin;
2283 InitialVector[3] = 1.0;
2284
2285 TransformationMatrix[0][0] = DirCos.XUnit.X;
2286 TransformationMatrix[0][1] = DirCos.XUnit.Y;
2287 TransformationMatrix[0][2] = DirCos.XUnit.Z;
2288 TransformationMatrix[1][0] = DirCos.YUnit.X;
2289 TransformationMatrix[1][1] = DirCos.YUnit.Y;
2290 TransformationMatrix[1][2] = DirCos.YUnit.Z;
2291 TransformationMatrix[2][0] = DirCos.ZUnit.X;
2292 TransformationMatrix[2][1] = DirCos.ZUnit.Y;
2293 TransformationMatrix[2][2] = DirCos.ZUnit.Z;
2294
2295 for (int i = 0; i < 4; ++i) {
2296 FinalVector[i] = 0.0;
2297 for (int j = 0; j < 4; ++j) {
2298 FinalVector[i] += TransformationMatrix[i][j] * InitialVector[j];
2299 }
2300 }
2301
2302 localPt.X = FinalVector[0];
2303 localPt.Y = FinalVector[1];
2304 localPt.Z = FinalVector[2];
2305 } // Point3D rotated
2306
2307 double Pot;
2308 Vector3D localF;
2309 double xpt = localPt.X;
2310 double ypt = localPt.Y;
2311 double zpt = localPt.Z;
2312
2313 double rW = radius;
2314 double lW = LZ;
2315
2316 // field point from element centroid
2317 double dist = sqrt(xpt * xpt + ypt * ypt + zpt * zpt);
2318
2319 if (dist >= FarField * lW) // all are distances and, hence, +ve
2320 {
2321 double dA = 2.0 * ST_PI * rW * lW;
2322 Pot = dA / dist;
2323 double dist3 = dist * dist * dist;
2324 localF.X = dA * xpt / dist3;
2325 localF.Y = dA * ypt / dist3;
2326 localF.Z = dA * zpt / dist3;
2327 } else if ((fabs(xpt) < MINDIST) && (fabs(ypt) < MINDIST) &&
2328 (fabs(zpt) < MINDIST)) {
2329 Pot = ExactCentroidalP_W(rW, lW);
2330 localF.X = 0.0; // CHECK - these flux values need to be confirmed
2331 localF.Y = 0.0;
2332 localF.Z = 0.0;
2333 } else if ((fabs(xpt) < MINDIST) && (fabs(ypt) < MINDIST)) {
2334 Pot = ExactAxialP_W(rW, lW, localPt.Z);
2335 localF.X = localF.Y = 0.0;
2336 localF.Z = ExactThinFZ_W(rW, lW, xpt, ypt, zpt);
2337 } else {
2338 Pot = ExactThinP_W(rW, lW, xpt, ypt, zpt);
2339 localF.X = ExactThinFX_W(rW, lW, xpt, ypt, zpt);
2340 localF.Y = ExactThinFY_W(rW, lW, xpt, ypt, zpt);
2341 localF.Z = ExactThinFZ_W(rW, lW, xpt, ypt, zpt);
2342 }
2343
2344 (*globalF) = RotateVector3D(&localF, &DirCos, local2global);
2345 return (Pot);
2346} // LineKnChPF ends
ExactThinFX_W
double ExactThinFX_W(double rW, double lW, double X, double Y, double Z)
Definition: Isles.c:2103
ExactCentroidalP_W
double ExactCentroidalP_W(double rW, double lW)
Definition: Isles.c:1779
ExactThinFY_W
double ExactThinFY_W(double rW, double lW, double X, double Y, double Z)
Definition: Isles.c:2121
ExactAxialP_W
double ExactAxialP_W(double rW, double lW, double Z)
Definition: Isles.c:1789
ExactThinFZ_W
double ExactThinFZ_W(double rW, double lW, double X, double Y, double Z)
Definition: Isles.c:2139
ExactThinP_W
double ExactThinP_W(double rW, double lW, double X, double Y, double Z)
Definition: Isles.c:2088
GetDistancePoint3D
double GetDistancePoint3D(Point3D *a, Point3D *b)
Definition: Vector.c:192
UnitVector3D
Vector3D UnitVector3D(Vector3D *v)
Definition: Vector.c:227

Referenced by KnChPFAtPoint().

◆ PointKnChPF()

double PointKnChPF ( Point3D  SourcePt,
Point3D  FieldPt,
Vector3D *  globalF 
)

Definition at line 2177 of file Isles.c.

2177 {
2178 double dist = GetDistancePoint3D(&SourcePt, &FieldPt);
2179 double dist3 = pow(dist, 3.0);
2180
2181 if (dist3 < MINDIST3) {
2182 globalF->X = 0.0;
2183 globalF->Y = 0.0;
2184 globalF->Z = 0.0;
2185 } else {
2186 globalF->X = (FieldPt.X - SourcePt.X) / dist3;
2187 globalF->Y = (FieldPt.Y - SourcePt.Y) / dist3;
2188 globalF->Z = (FieldPt.Z - SourcePt.Z) / dist3;
2189 }
2190
2191 if (dist < MINDIST)
2192 return (0.0);
2193 else
2194 return (1.0 / dist);
2195} // PointKnChPF ends
MINDIST3
#define MINDIST3
Definition: Isles.h:19

Referenced by KnChPFAtPoint().

◆ Sign()

int Sign ( double  x)

Definition at line 2506 of file Isles.c.

2506 {
2507 if (fabs(x) < MINDIST)
2508 return (0);
2509 else
2510 return (x < 0 ? -1 : 1);
2511} // Sign ends

Referenced by ExactRecSurf(), and ExactTriSurf().

◆ VolumeKnChPF()

double VolumeKnChPF ( int  NbPts,
Point3D *  Vertices,
Point3D  FieldPt,
Vector3D *  globalF 
)

Definition at line 2491 of file Isles.c.

2492 {
2493 printf("VolumeKnChPF not implemented yet ... returning zero flux\n");
2494 {
2495 globalF->X = 0.0;
2496 globalF->Y = 0.0;
2497 globalF->Z = 0.0;
2498 }
2499
2500 printf("VolumeKnChPF not implemented yet ... returning 0.0\n");
2501 return (0.0);
2502} // VolumeKnChPF ends

Referenced by KnChPFAtPoint().