Garfield++ 5.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 "Isles.h"
#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_sf.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)
 
ISLESGLOBAL int ExactRingPF (double a, Point3D localPt, double *potential, Vector3D *Flux)
 
double PointKnChPF (Point3D SourcePt, Point3D FieldPt, Vector3D *globalF)
 
double LineKnChPF (Point3D LineStart, Point3D LineStop, Point3D FieldPt, Vector3D *globalF)
 
double WireKnChPF (Point3D WireStart, Point3D WireStop, double radius, Point3D FieldPt, Vector3D *globalF)
 
double AreaKnChPF (int NbVertices, Point3D *Vertex, Point3D FieldPt, Vector3D *globalF)
 
double ApproxVolumeKnChPF (int, Point3D *, Point3D, Vector3D *globalF)
 
double VolumeKnChPF (int, Point3D *, Point3D, Vector3D *globalF)
 
int Sign (double x)
 

Macro Definition Documentation

◆ ARMAX

#define ARMAX   10000.0

Definition at line 17 of file Isles.c.

Referenced by ExactRecSurf(), and ExactTriSurf().

◆ ARMIN

#define ARMIN   0.0001

Definition at line 18 of file Isles.c.

Referenced by ExactRecSurf(), and ExactTriSurf().

◆ DEFINE_ISLESGLOBAL

#define DEFINE_ISLESGLOBAL

Definition at line 5 of file Isles.c.

◆ SHIFT

#define SHIFT   2.0

Definition at line 16 of file Isles.c.

Referenced by ExactRecSurf(), and ExactTriSurf().

◆ XNSegApprox

#define XNSegApprox   100

Definition at line 22 of file Isles.c.

Referenced by ExactRecSurf(), and ExactTriSurf().

◆ ZNSegApprox

#define ZNSegApprox   100

Definition at line 23 of file Isles.c.

Referenced by ExactRecSurf(), and ExactTriSurf().

Function Documentation

◆ ApproxFX_W()

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

Definition at line 1837 of file Isles.c.

1838 {
1839 if (DebugISLES) printf("In ApproxFX_W ...\n");
1840 ++ApproxCntr;
1841
1842 double dz = lW / zseg;
1843 double area = 2.0 * ST_PI * rW * dz;
1844
1845 double Fx = 0.0;
1846 double z0 = -0.5 * lW + 0.5 * dz;
1847 for (int k = 1; k <= zseg; ++k) {
1848 double zk = z0 + (k - 1) * dz;
1849 double dist = sqrt(X * X + Y * Y + (Z - zk) * (Z - zk));
1850 double dist3 = pow(dist, 3.0);
1851 if (fabs(dist) >= MINDIST) {
1852 Fx += (area * X / dist3);
1853 }
1854 } // zseg
1855
1856 return (Fx);
1857} // ApproxFX_W ends
ApproxCntr
ISLESGLOBAL int ApproxCntr
Definition Isles.h:35
DebugISLES
ISLESGLOBAL int DebugISLES
Definition Isles.h:36
ST_PI
#define ST_PI
Definition Isles.h:25
MINDIST
#define MINDIST
Definition Isles.h:18
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 1859 of file Isles.c.

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

◆ ApproxFZ_W()

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

Definition at line 1881 of file Isles.c.

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

◆ ApproxP_W()

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

Definition at line 1816 of file Isles.c.

1816 {
1817 // Approximate potential due to a wire segment
1818 if (DebugISLES) printf("In ApproxP_W ...\n");
1819 ++ApproxCntr;
1820
1821 double dz = lW / zseg;
1822 double area = 2.0 * ST_PI * rW * dz;
1823
1824 double Pot = 0.0;
1825 double z0 = -0.5 * lW + 0.5 * dz;
1826 for (int k = 1; k <= zseg; ++k) {
1827 double zk = z0 + (k - 1) * dz;
1828 double dist = sqrt(X * X + Y * Y + (Z - zk) * (Z - zk));
1829 if (fabs(dist) >= MINDIST) {
1830 Pot += area / dist;
1831 }
1832 } // zseg
1833
1834 return (Pot);
1835} // 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 if (DebugISLES) {
710 printf("In ApproxRecSurf ...\n");
711 }
712
713 ++ApproxCntr;
714
715 double dx = (xhi - xlo) / xseg;
716 double dz = (zhi - zlo) / zseg;
717 double xel = (xhi - xlo) / xseg;
718 double zel = (zhi - zlo) / zseg;
719 double diag = sqrt(dx * dx + dz * dz);
720 double area = xel * zel;
721
722 double Pot = 0., XFlux = 0., YFlux = 0., ZFlux = 0.;
723
724 if (area > MINDIST2) { // else not necessary
725 for (int i = 1; i <= xseg; ++i) {
726 double xi = xlo + (dx / 2.0) + (i - 1) * dx;
727 for (int k = 1; k <= zseg; ++k) {
728 double zk = zlo + (dz / 2.0) + (k - 1) * dz;
729
730 double dist = sqrt((X - xi) * (X - xi) + Y * Y + (Z - zk) * (Z - zk));
731 if (DebugISLES) printf("dist: %lg\n", dist);
732 if (dist >= diag) {
733 Pot += area / dist;
734 } else if (dist <= MINDIST) {
735 // Self influence
736 Pot += 2.0 * (xel * log((zel + sqrt(xel * xel + zel * zel)) / xel) +
737 zel * log((xel + sqrt(xel * xel + zel * zel)) / zel));
738 } else {
739 // in the intermediate region where diag > dist > MINDIST
740 Pot += area / diag; // replace by expression of self-influence
741 if (DebugISLES) printf("Special Pot: %lg\n", area / diag);
742 }
743
744 if (dist >= diag) {
745 double f = area / (dist * dist * dist);
746 XFlux += f * (X - xi);
747 YFlux += f * Y;
748 ZFlux += f * (Z - zk);
749 } else {
750 double f = area / (diag * diag * diag);
751 XFlux += f * (X - xi);
752 YFlux += f * Y;
753 ZFlux += f * (Z - zk);
754 if (DebugISLES) {
755 printf("Special XFlux: %lg, YFlux: %lg, ZFlux: %lg\n", f * (X - xi),
756 f * Y, f * (Z - zk));
757 }
758 } // else dist >= diag
759 } // zseg
760 } // xseg
761 } // if area > MINDIST2
762
763 *Potential = Pot;
764 Flux->X = XFlux;
765 Flux->Y = YFlux;
766 Flux->Z = ZFlux;
767
768 return 0;
769} // ApproxRecSurf ends
MINDIST2
#define MINDIST2
Definition Isles.h:19
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 1666 of file Isles.c.

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

Referenced by ExactTriSurf().

◆ ApproxVolumeKnChPF()

double ApproxVolumeKnChPF ( int NbPts,
Point3D * SourcePt,
Point3D FieldPt,
Vector3D * globalF )

Definition at line 2827 of file Isles.c.

2829 {
2830 printf("ApproxVolumeKnChPF not implemented yet ... returning zero flux\n");
2831 globalF->X = 0.0;
2832 globalF->Y = 0.0;
2833 globalF->Z = 0.0;
2834
2835 printf("ApproxVolumeKnChPF not implemented yet ... returning 0.0\n");
2836 return 0.0;
2837} // ApproxVolumeKnChPF ends

◆ ApproxWire()

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

Definition at line 1903 of file Isles.c.

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

◆ AreaKnChPF()

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

Definition at line 2694 of file Isles.c.

2695 {
2696 if (NbVertices != 4) {
2697 printf(
2698 "Only rectangular areas with known charges allowed at present ...\n");
2699 exit(-1);
2700 }
2701
2702 double xvert[4], yvert[4], zvert[4];
2703 xvert[0] = Vertex[0].X;
2704 xvert[1] = Vertex[1].X;
2705 xvert[2] = Vertex[2].X;
2706 xvert[3] = Vertex[3].X;
2707 yvert[0] = Vertex[0].Y;
2708 yvert[1] = Vertex[1].Y;
2709 yvert[2] = Vertex[2].Y;
2710 yvert[3] = Vertex[3].Y;
2711 zvert[0] = Vertex[0].Z;
2712 zvert[1] = Vertex[1].Z;
2713 zvert[2] = Vertex[2].Z;
2714 zvert[3] = Vertex[3].Z;
2715
2716 double xfld = FieldPt.X;
2717 double yfld = FieldPt.Y;
2718 double zfld = FieldPt.Z;
2719
2720 double xorigin = 0.25 * (xvert[3] + xvert[2] + xvert[1] + xvert[0]);
2721 double yorigin = 0.25 * (yvert[3] + yvert[2] + yvert[1] + yvert[0]);
2722 double zorigin = 0.25 * (zvert[3] + zvert[2] + zvert[1] + zvert[0]);
2723
2724 // lengths of the sides
2725 double SurfLX = sqrt((xvert[1] - xvert[0]) * (xvert[1] - xvert[0]) +
2726 (yvert[1] - yvert[0]) * (yvert[1] - yvert[0]) +
2727 (zvert[1] - zvert[0]) * (zvert[1] - zvert[0]));
2728 double SurfLZ = sqrt((xvert[2] - xvert[1]) * (xvert[2] - xvert[1]) +
2729 (yvert[2] - yvert[1]) * (yvert[2] - yvert[1]) +
2730 (zvert[2] - zvert[1]) * (zvert[2] - zvert[1]));
2731
2732 // Direction cosines - CHECK!
2733 //
2734 // 3 2
2735 // ________
2736 // | |
2737 // | |
2738 // | |
2739 // -------- -> +ve x-axis
2740 // 0 1
2741 // |
2742 // |
2743 // V +ve z-axis
2744 // The resulting +ve y-axis is out of the paper towards the reader.
2745 DirnCosn3D DirCos;
2746 DirCos.XUnit.X = (xvert[1] - xvert[0]) / SurfLX;
2747 DirCos.XUnit.Y = (yvert[1] - yvert[0]) / SurfLX;
2748 DirCos.XUnit.Z = (zvert[1] - zvert[0]) / SurfLX;
2749 DirCos.ZUnit.X = (xvert[0] - xvert[3]) / SurfLZ;
2750 DirCos.ZUnit.Y = (yvert[0] - yvert[3]) / SurfLZ;
2751 DirCos.ZUnit.Z = (zvert[0] - zvert[3]) / SurfLZ;
2752 DirCos.YUnit = Vector3DCrossProduct(&DirCos.ZUnit, &DirCos.XUnit);
2753
2754 Point3D localPt;
2755 // Through InitialVector[], field point gets translated to ECS origin.
2756 // Axes direction are, however, still global which when rotated to ECS
2757 // system, yields FinalVector[].
2758 { // Rotate point3D from global to local system to get localPt.
2759 double InitialVector[3] = {xfld - xorigin, yfld - yorigin, zfld - zorigin};
2760 double TransformationMatrix[3][3] = {{0.0, 0.0, 0.0},
2761 {0.0, 0.0, 0.0},
2762 {0.0, 0.0, 0.0}};
2763 TransformationMatrix[0][0] = DirCos.XUnit.X;
2764 TransformationMatrix[0][1] = DirCos.XUnit.Y;
2765 TransformationMatrix[0][2] = DirCos.XUnit.Z;
2766 TransformationMatrix[1][0] = DirCos.YUnit.X;
2767 TransformationMatrix[1][1] = DirCos.YUnit.Y;
2768 TransformationMatrix[1][2] = DirCos.YUnit.Z;
2769 TransformationMatrix[2][0] = DirCos.ZUnit.X;
2770 TransformationMatrix[2][1] = DirCos.ZUnit.Y;
2771 TransformationMatrix[2][2] = DirCos.ZUnit.Z;
2772 double FinalVector[3];
2773
2774 for (int i = 0; i < 3; ++i) {
2775 FinalVector[i] = 0.0;
2776 for (int j = 0; j < 3; ++j) {
2777 FinalVector[i] += TransformationMatrix[i][j] * InitialVector[j];
2778 }
2779 }
2780
2781 localPt.X = FinalVector[0];
2782 localPt.Y = FinalVector[1];
2783 localPt.Z = FinalVector[2];
2784 } // Point3D rotated
2785
2786 double Pot;
2787 Vector3D localF;
2788 double xpt = localPt.X;
2789 double ypt = localPt.Y;
2790 double zpt = localPt.Z;
2791
2792 double a = SurfLX;
2793 double b = SurfLZ;
2794 double diag = sqrt(a * a + b * b); // diagonal of the element
2795
2796 // distance of field point from element centroid
2797 double dist = sqrt(xpt * xpt + ypt * ypt + zpt * zpt);
2798 double dist3 = dist * dist * dist;
2799
2800 if (dist >= FarField * diag) // all are distances and, hence, +ve
2801 {
2802 double dA = a * b;
2803 Pot = dA / dist;
2804 localF.X = xpt * dA / dist3;
2805 localF.Y = ypt * dA / dist3;
2806 localF.Z = zpt * dA / dist3;
2807 } else {
2808 // normalize distances by `a' while sending - likely to improve accuracy
2809 int fstatus =
2810 ExactRecSurf(xpt / a, ypt / a, zpt / a, -1.0 / 2.0, -(b / a) / 2.0,
2811 1.0 / 2.0, (b / a) / 2.0, &Pot, &localF);
2812 if (fstatus) { // non-zero
2813 printf("problem in computing Potential of rectangular element ... \n");
2814 printf("a: %lg, b: %lg, X: %lg, Y: %lg, Z: %lg\n", a, b, xpt, ypt, zpt);
2815 // printf("returning ...\n");
2816 // return -1; void function at present
2817 }
2818 Pot *= a; // rescale Potential - cannot be done outside because of if(dist)
2819 }
2820
2821 (*globalF) = RotateVector3D(&localF, &DirCos, local2global);
2822 return (Pot);
2823} // 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:42
FarField
#define FarField
Definition Isles.h:22
RotateVector3D
Vector3D RotateVector3D(Vector3D *A, DirnCosn3D *DC, int Sense)
Definition Vector.c:395
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:77
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::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 ContinuityChUp(), ContinuityKnCh(), KnChPFAtPoint(), ValueChUp(), and ValueKnCh().

◆ ExactAxialFZ_W()

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

Definition at line 1805 of file Isles.c.

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

◆ ExactAxialP_W()

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

Definition at line 1792 of file Isles.c.

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

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

◆ ExactCentroidalP_W()

double ExactCentroidalP_W ( double rW,
double lW )

Definition at line 1782 of file Isles.c.

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

Referenced by WireKnChPF(), 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 42 of file Isles.c.

43 {
44 if (DebugISLES) printf("In ExactRecSurf ...\n");
45
46 ++IslesCntr;
47 ++ExactCntr;
48
49 ApproxFlag = 0;
50
51 if ((fabs(xhi - xlo) < 3.0 * MINDIST) || (fabs(zhi - zlo) < 3.0 * MINDIST)) {
52 fprintf(stdout, "Element size too small! ... returning ...\n");
53 return -1;
54 }
55
56 double ar = fabs((zhi - zlo) / (xhi - xlo));
57 if (ar > ARMAX || ar < ARMIN) {
58 fprintf(stdout, "Element too thin! ... returning ...\n");
59 return -2;
60 }
61
62 if (fabs(X) < MINDIST) X = 0.0;
63 if (fabs(Y) < MINDIST) Y = 0.0;
64 if (fabs(Z) < MINDIST) Z = 0.0;
65
66 double dxlo = X - xlo; // zero at the X=xlo edge
67 if (fabs(dxlo) < MINDIST) dxlo = 0.0;
68 double dxhi = X - xhi; // zero at the X=xhi edge
69 if (fabs(dxhi) < MINDIST) dxhi = 0.0;
70 double dzlo = Z - zlo; // zero at the Z=zlo edge
71 if (fabs(dzlo) < MINDIST) dzlo = 0.0;
72 double dzhi = Z - zhi; // zero at the Z=zhi edge
73 if (fabs(dzhi) < MINDIST) dzhi = 0.0;
74
75 // These four parameters can never be zero except at the four corners where
76 // one of them can become zero. For example, at X=xlo, Y=0, Z=zlo, D11
77 // is zero but the others are nonzero.
78 double D11 = sqrt(dxlo * dxlo + Y * Y + dzlo * dzlo);
79 if (fabs(D11) < MINDIST) D11 = 0.0;
80 double D21 = sqrt(dxhi * dxhi + Y * Y + dzlo * dzlo);
81 if (fabs(D21) < MINDIST) D21 = 0.0;
82 double D12 = sqrt(dxlo * dxlo + Y * Y + dzhi * dzhi);
83 if (fabs(D12) < MINDIST) D12 = 0.0;
84 double D22 = sqrt(dxhi * dxhi + Y * Y + dzhi * dzhi);
85 if (fabs(D22) < MINDIST) D22 = 0.0;
86
87 // Parameters related to the Y terms
88 int S1 = Sign(dzlo);
89 int S2 = Sign(dzhi);
90 double modY = fabs(Y);
91 int SY = Sign(Y);
92 double I1 = dxlo * modY;
93 double I2 = dxhi * modY;
94 double R1 = Y * Y + dzlo * dzlo;
95 double R2 = Y * Y + dzhi * dzhi;
96 if (fabs(I1) < MINDIST2) I1 = 0.0;
97 if (fabs(I2) < MINDIST2) I2 = 0.0;
98 if (fabs(R1) < MINDIST2) R1 = 0.0;
99 if (fabs(R2) < MINDIST2) R2 = 0.0;
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:18
ZNSegApprox
#define ZNSegApprox
Definition Isles.c:23
ARMAX
#define ARMAX
Definition Isles.c:17
Sign
int Sign(double x)
Definition Isles.c:2854
SHIFT
#define SHIFT
Definition Isles.c:16
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:22
IslesCntr
ISLESGLOBAL int IslesCntr
Definition Isles.h:35
fIsles
ISLESGLOBAL FILE * fIsles
Definition Isles.h:37
FailureCntr
ISLESGLOBAL int FailureCntr
Definition Isles.h:35
ApproxFlag
ISLESGLOBAL int ApproxFlag
Definition Isles.h:36
ExactCntr
ISLESGLOBAL int ExactCntr
Definition Isles.h:35

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

◆ ExactRingPF()

ISLESGLOBAL int ExactRingPF ( double a,
Point3D localPt,
double * potential,
Vector3D * Flux )

Definition at line 2193 of file Isles.c.

2194 {
2195 int dbgFn = 0;
2196
2197 // Coordinate transformation as performed for lines is necessary
2198 // Origin of the system is the center of the ring
2199 // The plane of the ring will need to be identified
2200 // An axis vertical to the plane will be the Y-axis
2201 // An axis orthogonal to the Y-axis will be X-axis (as done for lines)
2202 // Cross between the X and Y units will yield unit Z
2203 // Input parameter will need to be added - we now have only the radius
2204 // A way of identifying the plane of the ring is needed. It is quite likel
2205 // that the normal to the ring plane will be provided as an input. That
2206 // will immediately give us unit Y. The rest can be easily deduced.
2207
2208 double roe = localPt.X * localPt.X + localPt.Y * localPt.Y;
2209 roe = sqrt(roe);
2210 double z = localPt.Z;
2211 double r1dot = ((a + roe) * (a + roe)) + z * z;
2212 r1dot = sqrt(r1dot);
2213 double u = 2.0 * sqrt(a * roe) / r1dot;
2214
2215 // K1 and K2 are complete elliptic integrals
2216 // K implies first kind, according to GSL (and Wikipedia) convention
2217 gsl_mode_t mode = GSL_PREC_DOUBLE;
2218 double K1, K2;
2219 // Following are the routines with debugging information
2220 if (dbgFn) {
2221 gsl_sf_result rK1, rK2;
2222 gsl_sf_ellint_Kcomp_e(u, mode, &rK1);
2223 gsl_sf_ellint_Ecomp_e(u, mode, &rK2);
2224 K1 = rK1.val;
2225 K2 = rK2.val;
2226 } else { // If debugging information is not necessary
2227 K1 = gsl_sf_ellint_Kcomp(u, mode);
2228 K2 = gsl_sf_ellint_Ecomp(u, mode);
2229 }
2230
2231 double Vring = (a / ST_PI) * K1 / r1dot;
2232 // field in the radial direction
2233 double Eringroe = (a / ST_PI) * ((1.0 / (2.0 * r1dot * roe)) *
2234 (K1 - (((a * a - roe * roe + z * z) * K2) /
2235 (r1dot * r1dot * (1.0 - u * u)))));
2236 // field in the z direction
2237 double Eringz = (a / ST_PI) *
2238 ((z * K2) / ((r1dot * r1dot * r1dot) * (1.0 - u * u)));
2239
2240 *potential = Vring;
2241 Flux->X = Eringroe * localPt.X / roe;
2242 Flux->Y = Eringroe * localPt.Y / roe;
2243 Flux->Z = Eringz;
2244
2245 return 0;
2246} // ExactRingPF ends

◆ ExactThinFX_W()

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

Definition at line 2108 of file Isles.c.

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

Referenced by WireFlux(), and WireKnChPF().

◆ ExactThinFY_W()

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

Definition at line 2128 of file Isles.c.

2128 {
2129 if (DebugISLES) {
2130 printf("In ExactThinFY_W ...\n");
2131 printf("rW: %lg, lW: %lg, X: %lg, Y: %lg, Z: %lg\n", rW, lW, X, Y, Z);
2132 }
2133 // Expressions from PFXYZ_thin.m
2134 double h = 0.5 * lW;
2135 double Fy = 2.0 * Y *
2136 (sqrt(X * X + Y * Y + Z * Z + h * h + 2.0 * Z * h) * h -
2137 sqrt(X * X + Y * Y + Z * Z + h * h + 2.0 * Z * h) * Z +
2138 sqrt(X * X + Y * Y + Z * Z - 2.0 * Z * h + h * h) * h +
2139 sqrt(X * X + Y * Y + Z * Z - 2.0 * Z * h + h * h) * Z) /
2140 (X * X + Y * Y) /
2141 sqrt(X * X + Y * Y + Z * Z - 2.0 * Z * h + h * h) /
2142 sqrt(X * X + Y * Y + Z * Z + h * h + 2.0 * Z * h) * ST_PI;
2143 return (rW * Fy);
2144} // ExactThinFY_W ends

Referenced by WireFlux(), and WireKnChPF().

◆ ExactThinFZ_W()

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

Definition at line 2148 of file Isles.c.

2148 {
2149 if (DebugISLES) {
2150 printf("In ExactThinFZ_W ...\n");
2151 printf("rW: %lg, lW: %lg, X: %lg, Y: %lg, Z: %lg\n", rW, lW, X, Y, Z);
2152 }
2153
2154 // Expressions from PFXYZ_thin.m
2155 double h = 0.5 * lW;
2156 double Fz = 2.0 *
2157 (sqrt(X * X + Y * Y + Z * Z + h * h + 2.0 * Z * h) -
2158 sqrt(X * X + Y * Y + Z * Z - 2.0 * Z * h + h * h)) /
2159 sqrt(X * X + Y * Y + Z * Z - 2.0 * Z * h + h * h) /
2160 sqrt(X * X + Y * Y + Z * Z + h * h + 2.0 * Z * h) * ST_PI;
2161 return (rW * Fz);
2162} // ExactThinFZ_W ends

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

◆ ExactThinP_W()

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

Definition at line 2092 of file Isles.c.

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

Referenced by WireKnChPF(), and WirePot().

◆ ExactThinWire()

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

Definition at line 2164 of file Isles.c.

2165 {
2166 if (DebugISLES) {
2167 printf("In ExactThinWire_W ...\n");
2168 printf("rW: %lg, lW: %lg, X: %lg, Y: %lg, Z: %lg\n", rW, lW, X, Y, Z);
2169 }
2170 const double h = 0.5 * lW;
2171 const double r2 = X * X + Y * Y;
2172 const double a = h + Z;
2173 const double b = h - Z;
2174 const double c = sqrt(r2 + a * a);
2175 const double d = sqrt(r2 + b * b);
2176 const double s = 2. * ST_PI * rW;
2177 *potential = s * log((a + c) * (b + d) / r2);
2178 double f = s * (c * b + d * a) / (r2 * c * d);
2179 Flux->X = f * X;
2180 Flux->Y = f * Y;
2181 Flux->Z = s * (c - d) / (c * d);
2182 return 0;
2183}

Referenced by WirePF(), and WirePrimPF().

◆ ExactTriSurf()

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

Definition at line 783 of file Isles.c.

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

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 1950 of file Isles.c.

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

◆ ImprovedFY_W()

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

Definition at line 1974 of file Isles.c.

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

◆ ImprovedFZ_W()

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

Definition at line 1998 of file Isles.c.

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

◆ ImprovedP_W()

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

Definition at line 1934 of file Isles.c.

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

◆ ImprovedWire()

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

Definition at line 2028 of file Isles.c.

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

◆ LineKnChPF()

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

Definition at line 2271 of file Isles.c.

2272 {
2273 int debugFn = 0;
2274
2275 if (debugFn) {
2276 printf("In LineKnChPF:\n");
2277 printf("LineStart: %lg %lg %lg\n", LineStart.X, LineStart.Y, LineStart.Z);
2278 printf("LineStop: %lg %lg %lg\n", LineStop.X, LineStop.Y, LineStop.Z);
2279 printf("FieldPt: %lg %lg %lg\n", FieldPt.X, FieldPt.Y, FieldPt.Z);
2280 }
2281
2282 double xvert[2], yvert[2], zvert[2];
2283 xvert[0] = LineStart.X;
2284 xvert[1] = LineStop.X;
2285 yvert[0] = LineStart.Y;
2286 yvert[1] = LineStop.Y;
2287 zvert[0] = LineStart.Z;
2288 zvert[1] = LineStop.Z;
2289
2290 double xfld = FieldPt.X;
2291 double yfld = FieldPt.Y;
2292 double zfld = FieldPt.Z;
2293
2294 double xorigin = 0.5 * (xvert[1] + xvert[0]);
2295 double yorigin = 0.5 * (yvert[1] + yvert[0]);
2296 double zorigin = 0.5 * (zvert[1] + zvert[0]);
2297 double LZ = GetDistancePoint3D(&LineStop, &LineStart);
2298 if (debugFn) {
2299 printf("xorigin: %lg, yorigin: %lg, zorigin: %lg\n", xorigin, yorigin,
2300 zorigin);
2301 printf("LZ: %lg\n", LZ);
2302 }
2303
2304 // Create a local coordinate system for which the z axis is along the length
2305 // of the line
2306 // FieldPt needs to be transformed to localPt
2307
2308 // Find direction cosines of the line charge to initiate the transformation.
2309 DirnCosn3D DirCos;
2310
2311 // Copied from the DiscretizeWire function of neBEM/src/PreProcess/ReTrim.c
2312 // Direction cosines along the wire - note difference from surface primitives!
2313 // The direction along the wire is considered to be the z axis of the LCS.
2314 // So, let us fix that axial vector first
2315 DirCos.ZUnit.X = (xvert[1] - xvert[0]) / LZ; // useful
2316 DirCos.ZUnit.Y = (yvert[1] - yvert[0]) / LZ;
2317 DirCos.ZUnit.Z = (zvert[1] - zvert[0]) / LZ; // useful
2318 // Now direction cosines for the X and Y axes.
2319 {
2320 Vector3D XUnit, YUnit, ZUnit;
2321 ZUnit.X = DirCos.ZUnit.X;
2322 ZUnit.Y = DirCos.ZUnit.Y;
2323 ZUnit.Z = DirCos.ZUnit.Z;
2324
2325 // Find any vector that is orthogonal to the Z unit vector, and make that
2326 // the X-axis
2327 int caseComp = 1;
2328 double maxComp = fabs(ZUnit.X);
2329 if (fabs(ZUnit.Y) > maxComp) {
2330 caseComp = 2;
2331 maxComp = fabs(ZUnit.Y);
2332 }
2333 if (fabs(ZUnit.Z) > maxComp) {
2334 caseComp = 3;
2335 maxComp = fabs(ZUnit.Z);
2336 }
2337
2338 switch (caseComp) {
2339 case 1: // ZUnit.X is larger than the other direction cosines
2340 XUnit.Y = 1.0;
2341 XUnit.Z = 1.0;
2342 XUnit.X = (ZUnit.Y + ZUnit.Z) / ZUnit.X;
2343 break;
2344 case 2: // ZUnit.Y is larger than the other direction cosines
2345 XUnit.X = 1.0;
2346 XUnit.Z = 1.0;
2347 XUnit.Y = (ZUnit.X + ZUnit.Z) / ZUnit.Y;
2348 break;
2349 case 3: // ZUnit.Z is larger than the other direction cosines
2350 XUnit.X = 1.0;
2351 XUnit.Y = 1.0;
2352 XUnit.Z = (ZUnit.X + ZUnit.Y) / ZUnit.Z;
2353 break;
2354 default:;
2355 } // switch ends
2356 XUnit = UnitVector3D(&XUnit); // direction cosines for X is created
2357
2358 // y-Axis: vectorial product of axes 1 and 2.
2359 YUnit = Vector3DCrossProduct(&ZUnit, &XUnit);
2360 YUnit = UnitVector3D(&YUnit);
2361 // end of replacement
2362
2363 DirCos.XUnit.X = XUnit.X;
2364 DirCos.XUnit.Y = XUnit.Y;
2365 DirCos.XUnit.Z = XUnit.Z;
2366 DirCos.YUnit.X = YUnit.X;
2367 DirCos.YUnit.Y = YUnit.Y;
2368 DirCos.YUnit.Z = YUnit.Z;
2369 } // X and Y direction cosines computed
2370
2371 Point3D localPt;
2372 // Through InitialVector[], field point gets translated to ECS origin.
2373 // Axes direction are, however, still global which when rotated to ECS
2374 // system, yields FinalVector[].
2375 { // Rotate point3D from global to local system to get localPt.
2376 double InitialVector[3] = {xfld - xorigin, yfld - yorigin, zfld - zorigin};
2377 double TransformationMatrix[3][3] = {{0.0, 0.0, 0.0},
2378 {0.0, 0.0, 0.0},
2379 {0.0, 0.0, 0.0}};
2380 TransformationMatrix[0][0] = DirCos.XUnit.X;
2381 TransformationMatrix[0][1] = DirCos.XUnit.Y;
2382 TransformationMatrix[0][2] = DirCos.XUnit.Z;
2383 TransformationMatrix[1][0] = DirCos.YUnit.X;
2384 TransformationMatrix[1][1] = DirCos.YUnit.Y;
2385 TransformationMatrix[1][2] = DirCos.YUnit.Z;
2386 TransformationMatrix[2][0] = DirCos.ZUnit.X;
2387 TransformationMatrix[2][1] = DirCos.ZUnit.Y;
2388 TransformationMatrix[2][2] = DirCos.ZUnit.Z;
2389 double FinalVector[3];
2390
2391 for (int i = 0; i < 3; ++i) {
2392 FinalVector[i] = 0.0;
2393 for (int j = 0; j < 3; ++j) {
2394 FinalVector[i] += TransformationMatrix[i][j] * InitialVector[j];
2395 }
2396 }
2397
2398 localPt.X = FinalVector[0];
2399 localPt.Y = FinalVector[1];
2400 localPt.Z = FinalVector[2];
2401 } // Point3D rotated
2402
2403 double Pot;
2404 Vector3D localF;
2405 double xpt = localPt.X;
2406 double ypt = localPt.Y;
2407 double zpt = localPt.Z;
2408
2409 if (debugFn) {
2410 // printf("LineStart: %lg %lg %lg\n", LineStart.X, LineStart.Y,
2411 // LineStart.Z); printf("LineStop: %lg %lg %lg\n", LineStop.X, LineStop.Y,
2412 // LineStop.Z);
2413 printf("localPt: %lg %lg %lg\n", localPt.X, localPt.Y, localPt.Z);
2414 printf("LZ: %lg\n", LZ);
2415 }
2416
2417 // field point from element centroid
2418 double zptplus = zpt + 0.5 * LZ;
2419 double zptminus = zpt - 0.5 * LZ;
2420
2421 // NOTE: MatLab15, Maple Tanay's expressions are the same, only rearranged in
2422 // different ways. In MatLab, there are several possibilities indicated, but
2423 // Maple produces only the reasonable ones.
2424
2425 /*
2426 // Maple expressions: DefLineElementTry1.mw
2427 gsl_complex tmpcmp;
2428 tmpcmp.dat[0] = zptplus / sqrt(xpt*xpt + ypt*ypt); tmpcmp.dat[1] = 0.0;
2429 gsl_complex Pot1 = gsl_complex_arcsinh(tmpcmp);
2430 tmpcmp.dat[0] = zptminus / sqrt(xpt*xpt + ypt*ypt); tmpcmp.dat[1] = 0.0;
2431 gsl_complex Pot2 = gsl_complex_arcsinh(tmpcmp);
2432 Pot = Pot1.dat[0] - Pot2.dat[0];
2433 // The following expressions need to be rewritten after consulting the Maple
2434 // print-out.
2435 localF.X = xpt *((zptplus*sqrt(xpt*xpt + ypt*ypt + zptminus*zptminus) -
2436 zptminus*sqrt(xpt*xpt + ypt*ypt + zptplus*zptplus)) /
2437 (sqrt(xpt*xpt + ypt*ypt + zptminus*zptminus) *
2438 sqrt(xpt*xpt + ypt*ypt + zptplus*zptplus)) *
2439 (1.0/sqrt(xpt*xpt + ypt*ypt));
2440 localF.Y = ypt *((zptplus*sqrt(xpt*xpt + ypt*ypt + zptminus*zptminus) -
2441 zptminus*sqrt(xpt*xpt + ypt*ypt + zptplus*zptplus)) /
2442 (sqrt(xpt*xpt + ypt*ypt + zptminus*zptminus) *
2443 sqrt(xpt*xpt + ypt*ypt + zptplus*zptplus) *
2444 sqrt(xpt*xpt + ypt*ypt));
2445 localF.Z = (sqrt(xpt*xpt + ypt*ypt + zptplus*zptplus) -
2446 sqrt(xpt*xpt + ypt*ypt + zptminus*zptminus)) /
2447 (sqrt(xpt*xpt + ypt*ypt + zptplus*zptplus) *
2448 sqrt(xpt*xpt + ypt*ypt + zptminus*zptminus));
2449
2450 if (debugFn) {
2451 printf("Using Maple expressions =>\n");
2452 printf("Pot1: %lg + i %lg\n", Pot1.dat[0], Pot1.dat[1]);
2453 printf("Pot2: %lg + i %lg\n", Pot2.dat[0], Pot2.dat[1]);
2454 printf("Pot: %lg\n", Pot);
2455 printf("localF: %lg %lg %lg\n", localF.X, localF.Y, localF.Z);
2456 }
2457 */
2458
2459 // MatLab expressions: PFXYZ_MatLab15_4.out
2460 double tmpxy = xpt * xpt + ypt * ypt;
2461 double tmpzplus = sqrt(zptplus * zptplus + xpt * xpt + ypt * ypt);
2462 double tmpzminus = sqrt(xpt * xpt + ypt * ypt + zptminus * zptminus);
2463 Pot =
2464 log(zptplus + tmpzplus) - log(zptminus + tmpzminus);
2465 localF.X = xpt * ((zptplus / tmpxy / tmpzplus) -
2466 (zptminus / tmpxy / tmpzminus));
2467 localF.Y = ypt * ((zptplus / tmpxy / tmpzplus) -
2468 (zptminus / tmpxy / tmpzminus));
2469 localF.Z = (1.0 / tmpzminus) - (1.0 / tmpzplus);
2470 if (debugFn) {
2471 printf("Using MatLab expressions =>\n");
2472 printf("Pot: %lg\n", Pot);
2473 printf("localF: %lg %lg %lg\n", localF.X, localF.Y, localF.Z);
2474 }
2475
2476 /*
2477 // Tanay expressions - here the line is parallel to Z-axis, rather than Y.
2478 // According to Tanay's paper, the lines are parallel to Y.
2479 Pot = 0.0; // Tanay does not have an expression for the potential.
2480 // borrow from Matlab
2481 // Pot = log(zptplus + sqrt(zptplus*zptplus+xpt*xpt+ypt*ypt)) -
2482 // log(zptminus + sqrt(xpt*xpt+ypt*ypt+zptminus*zptminus));
2483 localF.X = (xpt / (xpt*xpt + ypt*ypt)) *
2484 ((zptplus / sqrt(zptplus*zptplus + xpt*xpt + ypt*ypt)) -
2485 (zptminus / sqrt(zptminus*zptminus + xpt*xpt + ypt*ypt)));
2486 localF.Y = (ypt / (xpt*xpt + ypt*ypt)) *
2487 ((zptplus / sqrt(zptplus*zptplus + xpt*xpt + ypt*ypt)) -
2488 (zptminus / sqrt(zptminus*zptminus + xpt*xpt + ypt*ypt)));
2489 localF.Z = - ((1.0 / sqrt(zptplus*zptplus + xpt*xpt + ypt*ypt)) -
2490 (1.0 / sqrt(zptminus*zptminus + xpt*xpt + ypt*ypt)));
2491 if (debugFn) {
2492 printf("Using Tanay expressions =>\n");
2493 printf("Pot: %lg\n", Pot);
2494 printf("localF: %lg %lg %lg\n", localF.X, localF.Y, localF.Z);
2495 }
2496 */
2497
2498 if (debugFn) {
2499 printf("Pot: %lg\n", Pot);
2500 printf("localF: %lg %lg %lg\n", localF.X, localF.Y, localF.Z);
2501 }
2502
2503 (*globalF) = RotateVector3D(&localF, &DirCos, local2global);
2504 if (debugFn) {
2505 printf("globalF: %lg %lg %lg\n", globalF->X, globalF->Y, globalF->Z);
2506 exit(-1);
2507 }
2508
2509 return (Pot);
2510} // LineKnChPF ends
GetDistancePoint3D
double GetDistancePoint3D(Point3D *a, Point3D *b)
Definition Vector.c:192
UnitVector3D
Vector3D UnitVector3D(Vector3D *v)
Definition Vector.c:227

Referenced by ContinuityChUp(), ContinuityKnCh(), KnChPFAtPoint(), ValueChUp(), and ValueKnCh().

◆ PointKnChPF()

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

Definition at line 2249 of file Isles.c.

2249 {
2250 double dist = GetDistancePoint3D(&SourcePt, &FieldPt);
2251 double dist3 = pow(dist, 3.0);
2252
2253 if (dist3 < MINDIST3) {
2254 globalF->X = 0.0;
2255 globalF->Y = 0.0;
2256 globalF->Z = 0.0;
2257 } else {
2258 globalF->X = (FieldPt.X - SourcePt.X) / dist3;
2259 globalF->Y = (FieldPt.Y - SourcePt.Y) / dist3;
2260 globalF->Z = (FieldPt.Z - SourcePt.Z) / dist3;
2261 }
2262
2263 if (dist < MINDIST)
2264 return (0.0);
2265 else
2266 return (1.0 / dist);
2267} // PointKnChPF ends
MINDIST3
#define MINDIST3
Definition Isles.h:20

Referenced by ContinuityChUp(), ContinuityKnCh(), KnChPFAtPoint(), ValueChUp(), and ValueKnCh().

◆ Sign()

int Sign ( double x)

Definition at line 2854 of file Isles.c.

2854 {
2855 if (fabs(x) < MINDIST)
2856 return (0);
2857 else
2858 return (x < 0 ? -1 : 1);
2859} // Sign ends

Referenced by ExactRecSurf(), and ExactTriSurf().

◆ VolumeKnChPF()

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

Definition at line 2841 of file Isles.c.

2842 {
2843 printf("VolumeKnChPF not implemented yet ... returning zero flux\n");
2844 globalF->X = 0.0;
2845 globalF->Y = 0.0;
2846 globalF->Z = 0.0;
2847
2848 printf("VolumeKnChPF not implemented yet ... returning 0.0\n");
2849 return 0.0;
2850} // VolumeKnChPF ends

Referenced by ContinuityChUp(), ContinuityKnCh(), KnChPFAtPoint(), ValueChUp(), and ValueKnCh().

◆ WireKnChPF()

double WireKnChPF ( Point3D WireStart,
Point3D WireStop,
double radius,
Point3D FieldPt,
Vector3D * globalF )

Definition at line 2514 of file Isles.c.

2515 {
2516 int debugFn = 1;
2517
2518 if (debugFn) {
2519 printf("In WireKnChPF:\n");
2520 printf("WireStart: %lg %lg %lg\n", WireStart.X, WireStart.Y, WireStart.Z);
2521 printf("WireStop: %lg %lg %lg\n", WireStop.X, WireStop.Y, WireStop.Z);
2522 printf("radius: %lg\n", radius);
2523 printf("FieldPt: %lg %lg %lg\n", FieldPt.X, FieldPt.Y, FieldPt.Z);
2524 }
2525
2526 double xvert[2], yvert[2], zvert[2];
2527 xvert[0] = WireStart.X;
2528 xvert[1] = WireStop.X;
2529 yvert[0] = WireStart.Y;
2530 yvert[1] = WireStop.Y;
2531 zvert[0] = WireStart.Z;
2532 zvert[1] = WireStop.Z;
2533
2534 double xfld = FieldPt.X;
2535 double yfld = FieldPt.Y;
2536 double zfld = FieldPt.Z;
2537
2538 double xorigin = 0.5 * (xvert[1] + xvert[0]);
2539 double yorigin = 0.5 * (yvert[1] + yvert[0]);
2540 double zorigin = 0.5 * (zvert[1] + zvert[0]);
2541 double LZ = GetDistancePoint3D(&WireStop, &WireStart);
2542 if (debugFn) {
2543 printf("xorigin: %lg, yorigin: %lg, zorigin: %lg\n", xorigin, yorigin,
2544 zorigin);
2545 printf("LZ: %lg\n", LZ);
2546 }
2547
2548 // Create a local coordinate system for which the z axis is along the length
2549 // of the wire
2550 // FieldPt needs to be transformed to localPt
2551
2552 // Find direction cosines of the line charge to initiate the transformation.
2553 DirnCosn3D DirCos;
2554
2555 // Copied from the DiscretizeWire function of neBEM/src/PreProcess/ReTrim.c
2556 // Direction cosines along the wire - note difference from surface primitives!
2557 // The direction along the wire is considered to be the z axis of the LCS.
2558 // So, let us first fix that as the axial vector.
2559 // Check whether this method is appropriate. If found not ok, we may need
2560 // to change ReTrim.c as well.
2561 // Finding the unit X and unit Y could work as in line elements.
2562 DirCos.ZUnit.X = (xvert[1] - xvert[0]) / LZ; // useful
2563 DirCos.ZUnit.Y = (yvert[1] - yvert[0]) / LZ;
2564 DirCos.ZUnit.Z = (zvert[1] - zvert[0]) / LZ; // useful
2565 // Now direction cosines for the X and Y axes.
2566 {
2567 Vector3D XUnit, YUnit, ZUnit;
2568 ZUnit.X = DirCos.ZUnit.X;
2569 ZUnit.Y = DirCos.ZUnit.Y;
2570 ZUnit.Z = DirCos.ZUnit.Z;
2571
2572 if (fabs(ZUnit.X) >= fabs(ZUnit.Z) && fabs(ZUnit.Y) >= fabs(ZUnit.Z)) {
2573 XUnit.X = -ZUnit.Y;
2574 XUnit.Y = ZUnit.X;
2575 XUnit.Z = 0.0;
2576 } else if (fabs(ZUnit.X) >= fabs(ZUnit.Y) &&
2577 fabs(ZUnit.Z) >= fabs(ZUnit.Y)) {
2578 XUnit.X = -ZUnit.Z;
2579 XUnit.Y = 0.0;
2580 XUnit.Z = ZUnit.X;
2581 } else {
2582 XUnit.X = 0.0;
2583 XUnit.Y = ZUnit.Z;
2584 XUnit.Z = -ZUnit.Y;
2585 }
2586 XUnit = UnitVector3D(&XUnit);
2587
2588 // y-Axis: vectorial product of axes 1 and 2.
2589 YUnit = Vector3DCrossProduct(&ZUnit, &XUnit);
2590 YUnit = UnitVector3D(&YUnit);
2591 // end of replacement
2592
2593 DirCos.XUnit.X = XUnit.X;
2594 DirCos.XUnit.Y = XUnit.Y;
2595 DirCos.XUnit.Z = XUnit.Z;
2596 DirCos.YUnit.X = YUnit.X;
2597 DirCos.YUnit.Y = YUnit.Y;
2598 DirCos.YUnit.Z = YUnit.Z;
2599 } // X and Y direction cosines computed
2600
2601 Point3D localPt;
2602 // Through InitialVector[], field point gets translated to ECS origin.
2603 // Axes direction are, however, still global which when rotated to ECS
2604 // system, yields FinalVector[].
2605 { // Rotate point3D from global to local system to get localPt.
2606 double InitialVector[3] = {xfld - xorigin, yfld - yorigin, zfld - zorigin};
2607 double TransformationMatrix[3][3] = {{0.0, 0.0, 0.0},
2608 {0.0, 0.0, 0.0},
2609 {0.0, 0.0, 0.0}};
2610 TransformationMatrix[0][0] = DirCos.XUnit.X;
2611 TransformationMatrix[0][1] = DirCos.XUnit.Y;
2612 TransformationMatrix[0][2] = DirCos.XUnit.Z;
2613 TransformationMatrix[1][0] = DirCos.YUnit.X;
2614 TransformationMatrix[1][1] = DirCos.YUnit.Y;
2615 TransformationMatrix[1][2] = DirCos.YUnit.Z;
2616 TransformationMatrix[2][0] = DirCos.ZUnit.X;
2617 TransformationMatrix[2][1] = DirCos.ZUnit.Y;
2618 TransformationMatrix[2][2] = DirCos.ZUnit.Z;
2619 double FinalVector[3];
2620
2621 for (int i = 0; i < 3; ++i) {
2622 FinalVector[i] = 0.0;
2623 for (int j = 0; j < 3; ++j) {
2624 FinalVector[i] += TransformationMatrix[i][j] * InitialVector[j];
2625 }
2626 }
2627
2628 localPt.X = FinalVector[0];
2629 localPt.Y = FinalVector[1];
2630 localPt.Z = FinalVector[2];
2631 } // Point3D rotated
2632
2633 double Pot;
2634 Vector3D localF;
2635 double xpt = localPt.X;
2636 double ypt = localPt.Y;
2637 double zpt = localPt.Z;
2638
2639 if (radius < MINDIST) radius = MINDIST;
2640 double rW = radius;
2641 double lW = LZ;
2642
2643 if (debugFn) {
2644 // printf("WireStart: %lg %lg %lg\n", WireStart.X, WireStart.Y,
2645 // WireStart.Z); printf("WireStop: %lg %lg %lg\n", WireStop.X, WireStop.Y,
2646 // WireStop.Z);
2647 printf("localPt: %lg %lg %lg\n", localPt.X, localPt.Y, localPt.Z);
2648 printf("rW: %lg, lW: %lg\n", rW, lW);
2649 }
2650
2651 // field point from element centroid
2652 double dist = sqrt(xpt * xpt + ypt * ypt + zpt * zpt);
2653
2654 if (dist >= FarField * lW) {
2655 double dA = 2.0 * ST_PI * rW * lW;
2656 Pot = dA / dist;
2657 double dist3 = dist * dist * dist;
2658 localF.X = dA * xpt / dist3;
2659 localF.Y = dA * ypt / dist3;
2660 localF.Z = dA * zpt / dist3;
2661 } else if ((fabs(xpt) < MINDIST) && (fabs(ypt) < MINDIST) &&
2662 (fabs(zpt) < MINDIST)) {
2663 Pot = ExactCentroidalP_W(rW, lW);
2664 localF.X = 0.0; // CHECK - these flux values need to be confirmed
2665 localF.Y = 0.0;
2666 localF.Z = 0.0;
2667 } else if ((fabs(xpt) < MINDIST) && (fabs(ypt) < MINDIST)) {
2668 Pot = ExactAxialP_W(rW, lW, zpt);
2669 localF.X = localF.Y = 0.0;
2670 localF.Z = ExactThinFZ_W(rW, lW, xpt, ypt, zpt);
2671 } else {
2672 Pot = ExactThinP_W(rW, lW, xpt, ypt, zpt);
2673 localF.X = ExactThinFX_W(rW, lW, xpt, ypt, zpt);
2674 localF.Y = ExactThinFY_W(rW, lW, xpt, ypt, zpt);
2675 localF.Z = ExactThinFZ_W(rW, lW, xpt, ypt, zpt);
2676 }
2677
2678 if (debugFn) {
2679 printf("Pot: %lg\n", Pot);
2680 printf("localF: %lg %lg %lg\n", localF.X, localF.Y, localF.Z);
2681 }
2682
2683 (*globalF) = RotateVector3D(&localF, &DirCos, local2global);
2684 if (debugFn) {
2685 printf("globalF: %lg %lg %lg\n", globalF->X, globalF->Y, globalF->Z);
2686 exit(-1);
2687 }
2688
2689 return (Pot);
2690} // WireKnChPF ends
ExactThinFX_W
double ExactThinFX_W(double rW, double lW, double X, double Y, double Z)
Definition Isles.c:2108
ExactCentroidalP_W
double ExactCentroidalP_W(double rW, double lW)
Definition Isles.c:1782
ExactThinFY_W
double ExactThinFY_W(double rW, double lW, double X, double Y, double Z)
Definition Isles.c:2128
ExactAxialP_W
double ExactAxialP_W(double rW, double lW, double Z)
Definition Isles.c:1792
ExactThinFZ_W
double ExactThinFZ_W(double rW, double lW, double X, double Y, double Z)
Definition Isles.c:2148
ExactThinP_W
double ExactThinP_W(double rW, double lW, double X, double Y, double Z)
Definition Isles.c:2092