Background models mathematical fixes

deerlab/bg_models.py

@@ -142,7 +142,7 @@ def _hom3dphase(t,conc,lam):
     conc = conc*1e-6*1e3*Nav # umol/L -> mol/L -> mol/m^3 -> spins/m^3
     # Compute background function
     ξ = 8*pi/9/np.sqrt(3)*(np.sqrt(3)+np.log(2-np.sqrt(3)))*D
-    B = np.exp(-1j*ξ*lam*conc*(t*1e-6))
+    B = np.exp(1j*ξ*lam*conc*(t*1e-6))
 
     return B
 # Create model
@@ -246,7 +246,7 @@ def _hom3dex_phase(t,conc,rex,lam):
 
     # Background function
     C_k = - Ic - np.squeeze(Vex*(dipolarkernel(t,rex,integralop=False,complex=True)).imag)
-    B = np.exp(-1j*lam*conc*C_k)
+    B = np.exp(1j*lam*conc*C_k)
 
     return B
 # Create model
@@ -284,7 +284,7 @@ def _homfractal(t,fconc,fdim,lam):
     elif d==4.5:
         κd = 5.35506  # d->4.5 limit of general expression
     else:
-        κd = 2/9*(-1)**(-d/3)*pi*np.cos(d*pi/6)*gamma(-d/3)*(
+        κd = 2/9*(-1)**(-d/3+1)*pi*np.cos(d*pi/6)*gamma(-d/3)*(
                 (-1 + (-1)**(d/3))*np.sqrt(3*np.pi)*gamma(1+d/3)/gamma(3/2+d/3)
                 + 6*hyp2f1_repro(1/2, -d/3, 3/2, 3)
               )
@@ -330,13 +330,13 @@ def _homfractal_phase(t,fconc,fdim,lam):
     elif d==4.5:
         ξd = 1j*np.inf # Limit of d->4.5 of equation below 
     else:
-        ξd = D**(d/3)*pi**(3/2)/9/gamma(3/2 + d/3) * (
+        ξd = 2*D**(d/3)*pi**(3/2)/9/gamma(3/2 + d/3) * (
                 np.sqrt(3)*pi*np.cos(d*pi/6)/np.cos(d*pi/3) 
                 - 2**(2+2*d/3)*3**(1 + d/3)*gamma(-1-2*d/3)*np.sin(d*pi/6)*gamma(3/2+d/3)*hyp2f1((-3-2*d)/6, -d/3, (3-2*d)/6, 1/3)/gamma((3-2*d)/6) 
             ) 
 
     # Compute background function
-    B = np.exp(-1j*ξd*fconc*lam*np.sign(t)*abs(t*1e-6)**(d/3))
+    B = np.exp(1j*ξd*fconc*lam*np.sign(t)*abs(t*1e-6)**(d/3))
     return B
 # ======================================================================
 # Create model

test/test_bgmodels.py

@@ -62,38 +62,56 @@ def test_bg_hom3d_value():
     conc = 143
     lam = 0.47
     B = dl.bg_hom3d(t, conc, lam)
-    Bref = 0.9227092
-    assert abs(B - Bref) < 1e-7
+    Bref = 0.92269778443975   # Value from numerical integration
+    assert abs(B - Bref) < 1e-4
 
 
-def test_bg_hom3d_phase():
-    Bref = 0.9998641918107823 - 0.01648022859583373j
-    assert_bgmodel_value(dl.bg_hom3d_phase, Bref)
-
+def test_bg_hom3d_phase_value():
+    t = 1.2
+    conc = 143
+    lam = 0.47
+    B = dl.bg_hom3d_phase(t, conc, lam)
+    Bref = 0.99994353+0.01062727j   # Value from numerical integration
+    assert np.abs(B - Bref) < 1e-4
 
-def test_bg_homfractal3d_value():
-    conc = 100
-    dim = 3
-    lam = 0.423
+def test_bg_homfractal_value():
+    conc = 1e-3
+    dim = 2.2
+    lam = 0.47
     t = 0.1
     B = dl.bg_homfractal(t, conc, dim, lam)
-    Bref = 0.99578995 
-    assert abs(B - Bref) < 1e-6
+    Bref = 0.9443678720378192    # Value from numerical integration
+    assert abs(B - Bref) < 1e-4
 
 
 def test_bg_homfractal_phase():
-    Bref = 0.999999988356966 - 0.00015259773164288322j
-    assert_bgmodel_value(dl.bg_homfractal_phase, Bref)
+    conc = 1e-3
+    dim = 2.2
+    lam = 0.47
+    t = 0.1
+    B = dl.bg_homfractal_phase(t, conc, dim, lam)
+    Bref = 0.9999771819749+0.006755407429425621j    # Value from numerical integration
+    assert abs(B - Bref) < 1e-4
 
 
 def test_bg_hom3dex():
-    Bref = 0.8828975362813715
-    assert_bgmodel_value(dl.bg_hom3dex, Bref)
+    t = 1.2
+    Rex = 3
+    conc = 143
+    lam = 0.47
+    B = dl.bg_hom3dex(t, conc, Rex, lam)
+    Bref = 0.926971275363862   # Value from numerical integration
+    assert abs(B - Bref) < 1e-4
 
 
 def test_bg_hom3dex_phase():
-    Bref = 0.9998641899814665 - 0.016480339581022234j
-    assert_bgmodel_value(dl.bg_hom3dex_phase, Bref)
+    t = 1.2
+    Rex = 3
+    conc = 143
+    lam = 0.47
+    B = dl.bg_hom3dex_phase(t, conc, Rex, lam)
+    Bref = 0.999943779807736+0.010603642007255394j   # Value from numerical integration
+    assert abs(B - Bref) < 1e-4
 
 
 def test_bg_exp_value():