Try to fix point dipole analytical solution to generate beneath the dipole position #42

mdovgialo · 2024-11-25T12:27:56Z

Zacznijmy od ostatniego. Proponuję w miejscu obliczeń tego potencjału zdefiniować
r_smaller = min(r, R)
r_bigger = max (r, R)

i zamiast bezpośrednio wstawiać r oraz R zamienić je na r_smaller i r_bigger tak, żeby dla r>R było jak obecnie a dla r<R odwrotnie.
Sprawdź i daj znać, czy działa

mdovgialo · 2024-12-03T16:42:33Z

Analytical solution - assumes there is no current flowing through skin outer layer, calculates potential in reference to infinity.

FEM - grounding in infinity and in the "neck" - electrode in the bottom of the 4 spheres. Otherwise, numerical precision doesn't allow to calculate difference between two point sources which generate VERY high amplitude potentials - due to our point dipole technical implementation as difference of two point sources.

Analytical solution has no air, FEM does.

After some experiments trying to fix, the best I've got:

Radial dipole:

Crossection through the dipole through Z axis:

the potential under the dipole R does not seem to follow FEM solution at all, but starts agreeng for r>R

Tangential dipole:

Seems to match FEM:

And on the crossection:

Except of the areas close to dipole, where analytical solution explodes...

mdovgialo · 2024-12-05T13:48:02Z

3 sphere model:
brain, skill, skin, air

Boundary conditions
FEM: grounding in infinity and in the "neck" of 3 sphere model
a) air domain is 100 meters in diameter
b) air domain is 0.5 meters in diameter

analytical: 0 current at outer sphere, 100 meter sphere

signal_simulator_kesi /home/mdovgialo/projects/halje_data_analysis/kESI/extras/data/generated/meshes/four_spheres_with_deep_electrode_in_air_boundary_bigger_0.005.msh -o /home/mdovgialo/test_dipole_below/simulations/sim_fem_x0_3_spheres -fs 1000 -l 0.100 -dx 2 -s /home/mdovgialo/projects/nencki_csd_tools/ncsd_tools/signal_simulator/example_data/point_dipole_source_test_z.csv -e /home/mdovgialo/projects/halje_data_analysis/kESI/extras/data/bundled/electrode_locations/10_20/10_20_FOUR_SPHERES_DEPTH_L2_0.088_neck.csv -c 0.33 0.33 0.0165 0.33 1e-10


signal_simulator_kesi /home/mdovgialo/projects/halje_data_analysis/kESI/extras/data/generated/meshes/four_spheres_with_deep_electrode_in_air_boundary0.005.msh -o /home/mdovgialo/test_dipole_below/simulations/sim_fem_x0_3_spheres_smaller -fs 1000 -l 0.100 -dx 2 -s /home/mdovgialo/projects/nencki_csd_tools/ncsd_tools/signal_simulator/example_data/point_dipole_source_test_z.csv -e /home/mdovgialo/projects/halje_data_analysis/kESI/extras/data/bundled/electrode_locations/10_20/10_20_FOUR_SPHERES_DEPTH_L2_0.088_neck.csv -c 0.33 0.33 0.0165 0.33 1e-10


signal_simulator_analytical -o /home/mdovgialo/test_dipole_below/simulations/sim_analytical_x0_test_rad_r_fix_v3_3_spheres -fs 1000 -l 0.100 -dx 2 -s /home/mdovgialo/projects/nencki_csd_tools/ncsd_tools/signal_simulator/example_data/point_dipole_source_test_z.csv -e /home/mdovgialo/projects/halje_data_analysis/kESI/extras/data/bundled/electrode_locations/10_20/10_20_FOUR_SPHERES_DEPTH_L2_0.088_neck.csv -r 0.082 0.086 0.09 100.0 -c 0.33 0.0165 0.33 1e-10

Result cross sections

python extras/draw_cross_section_nifty.py  /home/mdovgialo/test_dipole_below/simulations/sim_analytical_x0_test_rad_r_fix_v3_3_spheres/potential_profiles.nii.gz /home/mdovgialo/test_dipole_below/simulations/sim_fem_x0_3_spheres/potential_profiles.nii.gz /home/mdovgialo/test_dipole_below/simulations/sim_fem_x0_3_spheres_smaller/potential_profiles.nii.gz -f 0 -x=0.0 -y 0 -g 0.085

Analytical doesn't set potential in infinity to 0, so the "charges" can "accumulate"???

Let's try increasing analytical outer shell size to 1000000 meters (megameter):

  warnings.warn("no tangential dipole",
/home/mdovgialo/projects/halje_data_analysis/kESI/src/kesi/models/fourspheremodel.py:62: RuntimeWarning: overflow encountered in power
  Factor = ((self.r34 ** n - self.r43 ** (n + 1))
/home/mdovgialo/projects/halje_data_analysis/kESI/src/kesi/models/fourspheremodel.py:63: RuntimeWarning: overflow encountered in power
  / (k * self.r34 ** n + self.r43 ** (n + 1)))
/home/mdovgialo/projects/halje_data_analysis/kESI/src/kesi/models/fourspheremodel.py:62: RuntimeWarning: invalid value encountered in divide
  Factor = ((self.r34 ** n - self.r43 ** (n + 1))
/home/mdovgialo/projects/halje_data_analysis/kESI/src/kesi/models/fourspheremodel.py:327: RuntimeWarning: overflow encountered in power

Can't...

even 1000 meters causes an overflow...

(venv) mdovgialo@neuroinf011:~/projects/nencki_csd_tools$ signal_simulator_analytical -o /home/mdovgialo/test_dipole_below/simulations/sim_analytical_x0_test_rad_r_fix_v3_3_spheres_1000meters -fs 1000 -l 0.100 -dx 2 -s /home/mdovgialo/projects/nencki_csd_tools/ncsd_tools/signal_simulator/example_data/point_dipole_source_test_z.csv -e /home/mdovgialo/projects/halje_data_analysis/kESI/extras/data/bundled/electrode_locations/10_20/10_20_FOUR_SPHERES_DEPTH_L2_0.088_neck.csv -r 0.082 0.086 0.09 1000.0 -c 0.33 0.0165 0.33 1e-10

/home/mdovgialo/projects/halje_data_analysis/kESI/src/kesi/models/fourspheremodel.py:194: RuntimeWarning: no tangential dipole
  warnings.warn("no tangential dipole",
/home/mdovgialo/projects/halje_data_analysis/kESI/src/kesi/models/fourspheremodel.py:62: RuntimeWarning: overflow encountered in power
  Factor = ((self.r34 ** n - self.r43 ** (n + 1))
/home/mdovgialo/projects/halje_data_analysis/kESI/src/kesi/models/fourspheremodel.py:63: RuntimeWarning: overflow encountered in power
  / (k * self.r34 ** n + self.r43 ** (n + 1)))
/home/mdovgialo/projects/halje_data_analysis/kESI/src/kesi/models/fourspheremodel.py:62: RuntimeWarning: invalid value encountered in divide
  Factor = ((self.r34 ** n - self.r43 ** (n + 1))
^C^CTraceback (most recent call last):

But we can increase precision of the calculations from 64 bit floats to 128 bits...

Blue and orange are overlapping, this is weird...

let's add 0.2 and 0.5 meter air sphere for analytical:

Green and red are overlapping,

As we can see analytical solution in air pretty quickly settles, 0.5 meter radius air sphere with boundary condition of no current passing through the shell, is almost the same as 100 meters, and 1 megameters shell is basically equal to 100.

FEM solution which have boundary condition of 0 potential on the outside shell decay faster and all are consistently faster decaying in air than analytical...

mdovgialo · 2024-12-05T16:23:40Z

Adjusting n terms of the expansion in analytical solution:

n=100 takes a few minutes
n=500 takes 2.5 hours

for full sphere solution at 2 mm resolution.

We usually use n=100, but the problem is quadratic with n it seems:
With dipole being at r=0.03

we still see that at all tested n values there is something weird happening at r=0. But the n=500 doesn't have an instability around -0.03

At higher distances, at sphere surface n >= 10 is already quite accurate.

Let's look at n=500 in 3D:

The instabilities are way lower, compared to n=100 analytical solutions:

mdovgialo · 2024-12-12T17:01:23Z

Orange - tak jak było do poprawek
Niebieskie - próba reimplementacji w/g wzoru (6) dla wewnętrznej sfery
zielone - próba zmiany cos(theta) na takie theta które jest katem pomiedzy kierunkiem dipola a pozycją elektory, gdy układ współrzędnych na dipolu.

mdovgialo · 2024-12-13T16:30:17Z

analytyczne próby w nieskończonej przestrzeni:

Dipol skierowany do góry, w pozycji [0, 0, 0.030], o mocy p=100 / 1e6, kierunek dipola [0, 0, p]

Wzór klasyczny:

Wzór 1

Monopol?

Wzór 4:

Wzór na kropki:

Kod:

kESI/extras/solve_dipole_on_a_nifty.py

Line 56 in 1c7c039

    
           potential_flat = mag / (4 * np.pi * sigma) * (r - a * cos_theta) / np.cbrt(r**2 + a**2 - 2*a*r*cos_theta)

TODO:

it's not cube root! it's a square root and then cubed!!!!
just treat the classic formula as T2 and T1 should be the same as outside of the dipole in inner sphere

danek8317 · 2024-12-16T12:30:55Z

Regarding formula 4 above:

I lost b - the distance between the monopoles - should be q -> qb=p
I should have probably differentiate over a not r. This changes the formula to the same as above, that is
$\frac{p}{4 \pi \epsilon_0} \frac{a - r \cos \theta}{(r^2 + a^2 - 2ar \cos\theta)^{3/2}}$

Although I think I am still missing a minus :/

mdovgialo · 2024-12-19T16:46:33Z

Version 5, using classic free space solution + T1 is still slightly wrong:

mdovgialo self-assigned this Nov 25, 2024

mdovgialo mentioned this issue Jan 31, 2025

Draft: Dev/42 trying to fix dipole beneath, experiments #53

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Try to fix point dipole analytical solution to generate beneath the dipole position #42

Try to fix point dipole analytical solution to generate beneath the dipole position #42

mdovgialo commented Nov 25, 2024

mdovgialo commented Dec 3, 2024 •

edited

Loading

mdovgialo commented Dec 5, 2024 •

edited

Loading

mdovgialo commented Dec 5, 2024 •

edited

Loading

mdovgialo commented Dec 12, 2024

mdovgialo commented Dec 13, 2024 •

edited

Loading

danek8317 commented Dec 16, 2024 •

edited

Loading

mdovgialo commented Dec 19, 2024

Try to fix point dipole analytical solution to generate beneath the dipole position #42

Try to fix point dipole analytical solution to generate beneath the dipole position #42

Comments

mdovgialo commented Nov 25, 2024

mdovgialo commented Dec 3, 2024 • edited Loading

Radial dipole:

Tangential dipole:

mdovgialo commented Dec 5, 2024 • edited Loading

mdovgialo commented Dec 5, 2024 • edited Loading

mdovgialo commented Dec 12, 2024

mdovgialo commented Dec 13, 2024 • edited Loading

analytyczne próby w nieskończonej przestrzeni:

danek8317 commented Dec 16, 2024 • edited Loading

mdovgialo commented Dec 19, 2024

mdovgialo commented Dec 3, 2024 •

edited

Loading

mdovgialo commented Dec 5, 2024 •

edited

Loading

mdovgialo commented Dec 5, 2024 •

edited

Loading

mdovgialo commented Dec 13, 2024 •

edited

Loading

danek8317 commented Dec 16, 2024 •

edited

Loading