|
@@ -15,23 +15,23 @@ U_points = voxel_num * 3000
|
|
|
noise_ratio = mpf(0.0*1e-38) #1e8
|
|
|
|
|
|
total_periods=10#DO NOT CHANGE to fit T2s_scale autoscale
|
|
|
-rf_samples_per_period = 40
|
|
|
+rf_samples_per_period = 80
|
|
|
# max polynomial order equals rf_samples_per_period * total_periods
|
|
|
|
|
|
# B0=1.5T freq=64Mhz, period = 15.6 ns
|
|
|
period = mpf(10) #ms
|
|
|
omega = 2.0*mp.pi/period
|
|
|
#T2s_scale = 0.01 #ms # need to be 10ms
|
|
|
-T2s_scale = total_periods*period/rf_samples_per_period/voxel_num*5 #ms # need to be 10ms
|
|
|
-T2s_min = T2s_scale/10.0
|
|
|
+T2s_scale = total_periods*period/rf_samples_per_period/voxel_num*8 #ms # need to be 10ms
|
|
|
+T2s_min = T2s_scale/20.0
|
|
|
#print(period)
|
|
|
#ms
|
|
|
time_steps = np.array(mp.linspace(mpf(0), mpf(period*total_periods), rf_samples_per_period*total_periods))
|
|
|
-tmp = [mp.rand() for n in range(voxel_num)]
|
|
|
+tmp = [mp.rand()*0.0+mpf(1.0) for n in range(voxel_num)]
|
|
|
voxel_amplitudes = np.array(tmp)
|
|
|
tmp = [mp.rand() for n in range(voxel_num)]
|
|
|
-voxel_T2s_decay = np.array(tmp)*(T2s_scale-2*T2s_min) + T2s_min
|
|
|
-print(voxel_T2s_decay)
|
|
|
+voxel_T2s_decay = np.sort(np.array(tmp))*(T2s_scale-2*T2s_min) + T2s_min
|
|
|
+
|
|
|
voxel_all = np.append(voxel_amplitudes,voxel_T2s_decay/T2s_scale)
|
|
|
if voxel_num == 5:
|
|
|
# voxel_all = np.array([mpf(0.2),mpf(0.6),mpf(0.02),mpf(0.1)])
|
|
@@ -168,10 +168,19 @@ lambdas = GetLambda(mag_sin)
|
|
|
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)
|
|
|
x = np.linspace(0,1, U_points)
|
|
|
mag_x = np.linspace(0,1, len(mag_sin))
|
|
|
+import matplotlib as matplotlib
|
|
|
+matplotlib.rcParams.update({'font.size': 28})
|
|
|
+plt.figure(figsize=(30, 20))
|
|
|
plt.plot(x,U, lw=0.2)
|
|
|
plt.title(r"$\sum^M \lambda_i L_i(x)$", y=1.02)
|
|
|
+plt.xlim(0,1)
|
|
|
+plt.xlabel(sign, fontsize=28)
|
|
|
plt.savefig("plt.pdf")
|
|
|
plt.clf()
|
|
|
|