correct criteria Gnk

phase-encoding-method-of-moments.py

@@ -7,12 +7,12 @@ import matplotlib.pyplot as plt
 from mpmath import mp, mpf
 mp.dps = 2000
 
-voxel_num = 15
-phase_range = np.pi/2
-phase_init = np.pi/4 #(0.0)
+voxel_num = 5
+phase_range = np.pi/2 *0.0
+phase_init = np.pi/2  #(0.0)
 U_points = voxel_num * 3000
 
-noise_ratio = (0.0*1e-38) #1e8
+noise_ratio = (1e-60) #1e8
 
 total_periods=10#DO NOT CHANGE to fit T2s_scale autoscale
 rf_samples_per_period = 10
@@ -66,13 +66,15 @@ def gen_rf_signal(time):
             amp = voxel_amplitudes[i] * (
                 mp.exp(-time[t]/voxel_T2s_decay[i])
                 ) + ( 0.0 
-                          # + np.random.rand()*noise_ratio
+                     #     + np.random.rand()*noise_ratio
                         )
-            if t == 0:
+            # if t == 0:
                 #print("a_{:d}".format(i),float(voxel_amplitudes[i]* np.sin(voxel_phases[i] + phase_init)))
-                print("d_{:d}".format(i),float( np.exp(-(period/rf_samples_per_period)/voxel_T2s_decay[i]) ))
+                #print("d_{:d}".format(i),float( np.exp(-(period/rf_samples_per_period)/voxel_T2s_decay[i]) ))
             mag_sin[t] += amp * np.sin((voxel_phases[i] + phase_init))
             mag_cos[t] += amp * np.cos((voxel_phases[i] + phase_init))
+        mag_sin[t] += np.random.rand()*noise_ratio
+        #mag_cos[t] += np.random.rand()*noise_ratio
     return mag_sin, mag_cos
 
 def factorial(n):
@@ -101,7 +103,7 @@ def GetU (lambdas):
         P[n+1] = ((2*n+1)/(n+1)) * (2*x-1) * P[n]  - (n/(n+1))*P[n-1]
     U = np.zeros(U_points)
     for i in range (len(lambdas)):
-        if i%10 == 0: print(i)
+        # if i%10 == 0: print(i)
         U += lambdas[i]*P[i]*np.sqrt(2*i+1)
         #polyL = L[i] #shiftedLegendre(i)
         #U += lambdas[i]*polyL(x)
@@ -122,8 +124,24 @@ def GetLambda(mag_rf):
     return all_lambda
 
 
+def GetGnk(c,n,k):
+    Gnk = 0
+    for m in range(k+1):
+        Gnk+=(-1)**m * binom(k,m) * c[n+m]
+    return Gnk
+
+    
+
 mag_sin, mag_cos = gen_rf_signal(time_steps)
 
+for i in range(3):
+    print("c",i," = ",float(mag_sin[i]))
+
+for n in range(3):
+    for k in range (3):
+        print("G",n,k,float(GetGnk(mag_sin, n, k)))
+        
+
 sign = ""
 # for i in range(voxel_num):
 #     if i%5 == 0:
@@ -152,16 +170,16 @@ x = np.linspace(0,1, U_points)
 # plt.title("Legendre polynom of order "+str(len(L)))
 # plt.savefig("polyL.pdf")
 # plt.clf()
-print("Before lambda.")
+#print("Before lambda.")
 
 lambdas = GetLambda(mag_sin)
-print(lambdas)
+#print(lambdas)
 U = GetU(lambdas)
 
 sign ="" 
 for i in range(voxel_num):
     sign+="{:4.3g}, ".format(float( mp.exp(-(period/rf_samples_per_period)/voxel_T2s_decay[i]) ))
-print(sign)
+#print(sign)
 x = np.linspace(0,1, U_points)
 mag_x = np.linspace(0,1, len(mag_sin))
 import matplotlib as matplotlib