You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 

170 lines
4.5 KiB

import matplotlib.pyplot as plt
import numpy as np
from itertools import cycle
import argparse
import pickle
import yaml
from matplotlib import rc
rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
rc('text', usetex=True)
fig1, ax1 = plt.subplots(1,1,figsize=(8, 5))
lwidth = 2
font_size = 28
################ Flow Parameters
Rd = 2.5
Rt = 0.5
GradP = 4
mu = 0.5
fac = 1
nr = 50
VENC = 0.6
gamma = 267.513e6 # rad/Tesla/sec Gyromagnetic ratio for H nuclei
Bo = 1.5 # Tesla Magnetic Field Strenght
TE = 5e-3 # Echo-time
M = np.ones(nr) # Magnetization
phi0 = gamma*Bo*TE # Reference phase
phi02 = phi0%3.14
M1 = np.pi/(gamma*VENC)
ff = np.pi/(1000*gamma*M1)
uv = np.arange(-4*VENC,4*VENC,ff)
r = np.linspace(-Rd, Rd, nr)
dr = r[2]-r[1]
vmax = 1
v = vmax/Rt**2*( Rt**2 - r**2 )*(np.abs(r)<Rt); # Poiseuille Formula
ai = v/vmax
theta = np.linspace(-4,5,2000)
vtest = np.linspace(-5,5,2000)
JF = 0*theta
jv = 0*theta
JV = 0*theta
Mjv = np.zeros([len(theta),len(ai)])
jv0 = 0*theta
JV0 = 0*theta
Mjv0 = np.zeros([len(theta),len(ai)])
#################################### MAGNETIZACION FROM V
phiv = phi02 + v*np.pi/VENC
modv = np.ones(phiv.shape)
M1 = modv*np.cos(phi02) + 1j*modv*np.sin(phi02)
M2 = modv*np.cos(phiv) + 1j*modv*np.sin(phiv)
################################### FFT to COMPLEX M
S1 = np.fft.fft(M1)
S2 = np.fft.fft(M2)
################################### SubSampling
a1 = 0
a2 = 1
##### FILLED WITH ZEROS
US1 = S1
US2 = 0*S2
US2[a1::a2] = S2[a1::a2]
MR1 = np.fft.ifft(US1)
MR2 = np.fft.ifft(US2)
vrec1 = (np.angle(MR2)-phi02)*VENC/(np.pi)
for k in range(len(ai)):
jv0 = 1-np.cos(np.pi*(vrec1[k]-vtest)/VENC)
Mjv0[:,k] = jv0[:]
JV0 = JV0 + jv0
for k in range(len(ai)):
jv = 1-np.cos(np.pi*(vrec1[k]-theta*ai[k])/VENC)
Mjv[:,k] = jv[:]
JV = JV + jv
NJV1 = JV*100/np.max(JV)
MV = Mjv0
V =NJV1
left, bottom, width, height = [0.2, 0.2, 0.1, 0.1]
fig = plt.figure(figsize=(12, 6), dpi=100)
ax1 = plt.subplot(1,2,1)
ch1 = 20
ch2 = 23
ax0 = fig.add_axes([left, bottom, width, height])
ax0.plot(r,v,'b-')
ax0.plot([r[ch1]],[v[ch1]],color='xkcd:coral',marker='o')
ax0.plot([r[ch2]],[v[ch2]],color='xkcd:azure',marker='o')
ax0.set_xlim((-1.5,1.5))
#for k in range(22,39):
# if k!=ch1 and k!=ch2 and np.sum(MV[:,k])!=0:
# ax1.plot(vtest, MV[:,k],color='xkcd:beige',alpha=0.8)
ax1.plot(vtest, MV[:,ch1],color='xkcd:coral',label='$v_1$')
ax1.plot(vtest, MV[:,ch2],color='xkcd:azure',label='$v_2$')
m1x = vtest[np.where( np.abs(MV[:,ch1] - np.min(MV[:,ch1]))<0.001 )]
m1y = np.min(MV[:,ch1])
m2x = vtest[np.where( np.abs(MV[:,ch2] - np.min(MV[:,ch2]))<0.001 )]
#m2x = vtest[np.where(MV[:,ch2]==np.min(MV[:,ch2]))]
m2y = np.min(MV[:,ch2])
ax1.plot([m1x],[m1y],color='xkcd:coral',marker='o')
ax1.plot([m2x],[m2y],color='xkcd:azure',marker='o')
ax1.axvline(x=v[ch1], color='xkcd:coral', linestyle='--',label='$v_{1,true}$')
ax1.axvline(x=v[ch2], color='xkcd:azure', linestyle='--',label='$v_{2,true}$')
ax1.set_xlabel(r'$u$',fontsize=20)
ax1.set_ylabel(r'$J_i(u)$',fontsize=20)
#ax1.legend(loc='upper right', bbox_to_anchor=(0.5, 1.05),ncol=2, fancybox=True, shadow=True,fontsize=15)
ax1.set_yticks([])
ax1.set_xticks([])
ax1.set_xlim((-3.5,3.5))
ax1.set_ylim((-1.0,2.4))
ax2 = plt.subplot(1,2,2)
ax2.plot(theta,V,'b-')
ax2.axvline(x=1, color='k', linestyle='--')
ax2.set_xlabel(r'$\theta$',fontsize=20)
ax2.set_ylabel(r'$J_T(\theta)$',fontsize=20)
plt.yticks([])
ax2.set_xticks([])
plt.title(r'$\theta_{true}=1$' + '\n' +'$venc < v_{max}$',fontsize=15)
plt.xlim((-2,3))
plt.show()
#ax1.plot(u, J1, color = 'orangered', label = '$venc = 0.9 u_{true}$', linestyle='-',linewidth=lwidth)
#ax1.plot(u, J2, color = 'dodgerblue', label = '$venc = 0.6 u_{true}$', linestyle='-',linewidth=lwidth)
#ax1.axvline(x=1,color = 'black',linewidth = lwidth , label = '$u_{true}$')
ax1.legend(fontsize=20, loc= 'upper right')
ax1.tick_params(axis='both', which='major', labelsize=22)
ax1.set_yticks([])
ax1.set_xlabel('$u$',fontsize=font_size)
plt.show()
#fig1.savefig('functionals.png', dpi=500, bbox_inches='tight')