\documentclass[11pt,leqno]{article}
\usepackage{amssymb, amsmath} 
\usepackage{latexsym} 
\usepackage{graphics} 
%\usepackage{hyperref} 
%\hypersetup{breaklinks=true}
\usepackage{url} 
%\Urlmuskip=0mu plus 2mu
\usepackage{color} 


%\usepackage[]{showkeys}
%\usepackage[notcite,notref,color]{showkeys}


%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 	SPANISH
%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[spanish]{babel}
\usepackage[latin1]{inputenc}
\decimalpoint
%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 	FORMAT
%%%%%%%%%%%%%%%%%%%%%%%%%%
%\setlength{\textheight}{23.0cm}
\setlength{\textwidth}{15.75cm}
\setlength{\oddsidemargin}{0.5cm}
\setlength{\evensidemargin}{0.5cm}
%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 	LISTINGS
%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{listings}

\definecolor{dgreencolor}{rgb}{0.2,0.4,0.0}

\newcommand{\code}[1]{\text{\small\tt #1}}

\lstset{
	aboveskip = 12pt,
	belowskip = 10pt,
	basicstyle = {\fontsize{7pt}{7pt}\rmfamily},
	commentstyle = {\ttfamily\color{dgreencolor}},
	stringstyle = {\ttfamily\color{red}\bfseries},
	keywordstyle = {\ttfamily\color{blue}\bfseries},
	inputencoding = latin1,	
	showstringspaces = true, %only spaces in strings
	numbers = left,
	numberstyle = {\fontsize{7pt}{7pt}\ttfamily},
	numbersep = 8pt, %separation numbers - listing
	breaklines = true,
	breakatwhitespace = true, %breaks only at white spaces
	postbreak = \raisebox{0ex}[0ex][0ex]{\ensuremath{\color{red}\hookrightarrow\space}},
	breakindent = 40pt, %indentation after line break
	fontadjust = yes,
	xleftmargin = 2.5cm
}


%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 	COUNTERS
%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcounter{ejeno}
\setcounter{ejeno}{1}


%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 	MACROS
%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\R}{{\mathbb R}} 
\newcommand{\C}{{\mathbb C}} 
\newcommand{\N}{{\mathbb N}} 
\newcommand{\Z}{{\mathbb Z}} 
\newcommand{\Q}{{\mathbb Q}} 

\newcommand{\pvec}[1]{\vec{#1}\mkern2mu\vphantom{#1}'}

\newcommand{\ejer}{\ \par\vspace{-7pt}{\bf\number\value{ejeno})\ }\addtocounter{ejeno}{1}}



\begin{document}
\pagestyle{myheadings}
\markboth{}{\rm\textsf{\jobname}\qquad\hfill}

Esta hoja es la última de tu trabajo. Consta de dos partes bien diferenciadas. La primera,  consiste en algunos experimentos numéricos con la ecuación de Bloch. Para ganar tiempo, el programa te lo daré en su mayoría yo. La segunda, se centra sobre todo en uno de los métodos matemáticos para lograr la localización de los protones que están resonando.  

\

Considera el siguiente programa {\tt sage}. Lo puedes copiar a partir de la fuente de esta página. 

\begin{lstlisting}[language=python]
# Si no tienes sage instalado en tu ordenador, lo puedes
# ejecutar en https://sagecell.sagemath.org/

# Suponemos \gamma=B0=1 con las unidades elegidas,
# con lo que la frecuencia de resonancia es -1.

from sage.calculus.desolvers import desolve_odeint

# Ajusta aquí alfa, beta, la frecuencia y el L inicial
alp = 0.07
bet = 0.0
ome = -1
Linic = [0.0,0.0,1.0]

# B1 debe ser mucho más pequeño que B0
B1 = 0.05
# Se entiende que frecuencia cero, significa
# que hemos desconectado el campo externo
if ome == 0: B1 = 0

# Orientación de B cuando no hay giro
L30 = 1

# Cuando mayor sea N, más es la precisión y más lento es el programa
N = 100
# Se estudia la ecuación entre 0 y 200 unidades de tiempo
tn = srange(0.0,200.0,1.0/N)
L1,L2,L3,t=var('L1,L2,L3,t')

# Este es el segundo miembro de la ecuación 
f = [L2-B1*sin(ome*t)*L3-bet*L1, -L1+B1*cos(ome*t)*L3-bet*L2 ,B1*sin(ome*t)*L1-B1*cos(ome*t)*L2 -alp*(L3-L30), 1.0+0.0*t]

# Con esto se resuelve la ecuación,
# en L1n, L2n, L3n y tn,
# están respectivamente las listas de valores numéricos de 
# L1, L2, L3 y t
sol = desolve_odeint(f,Linic+[0.0],tn,[L1,L2,L3,t])
L1n, L2n, L3n = list(sol[:,0]), list(sol[:,1]), list(sol[:,2])

# Por ejemplo, lo siguiente dibuja L3 en función del tiempo
#show( line(zip(tn, L3n)) )
\end{lstlisting}


El propósito del primer ejercicio es sobre todo que haya una doble comprobación de que no me he equivocado al teclear la fórmula. El segundo ejercicio es para asegurarse de que no se ha colado ningún error al copiar el código en tu ordenador o en {\tt sage cell}.

\ejer
Comprueba que las tres primeras coordenadas de \code{f} en la línea 31 son justamente el segundo miembro de la ecuación de Bloch (con $\gamma=B_0=1$). 

\ejer
Quita el \code{\#} de comentario en la línea 41 y ejecuta el programa. Esta es la gráfica de~$L_3(t)$ y debería salir una figura en la que desde el valor inicial $L_3(0)=1$ acaba tendiendo a un valor casi nulo, porque con los parámetros tal como están, hay resonancia y la precesión ocurre prácticamente en el plano~$XY$. 

\medskip

El programa genera las listas \code{L1n}, \code{L2n} y \code{L3n} que dan las aproximaciones numéricas de las coordenadas de $\vec{L}(t)$ con $t\in [0,200]$, con una discretización consistente en \code{200N} puntos.  El reto informático es el siguiente:

\ejer
Representa la aproximación numérica de $\vec{L}(t)/\|\vec{L}(t)\|$ en la superficie de la esfera unidad con el ecuador destacado y las partes positivas de los tres ejes. Te mandaré un ejemplo por correo electrónico. Si se me olvida, recuérdamelo. 
\medskip

Para lo anterior, puedes completar el programa {\tt sage} anterior (los comandos \code{line3d} y \code{sphere} te serán de utilidad) o, si {\tt matlab} u otro programa te es más familiar, puedes descargar las listas \code{L1n}, \code{L2n} y \code{L3n} a ficheros y leerlas desde allí.
Si tus conocimientos informáticos impiden cualquier posibilidad, te daré la solución pero antes, por favor, inténtalo. 

\ejer
Con el programa que has hecho para representar $\vec{L}(t)/\|\vec{L}(t)\|$, obtén las figuras correspondientes a las elecciones de parámetros en las líneas 10--13 indicadas en las 8 situaciones siguientes, donde $\vec{e}_j$ son los vectores de la base canónica. 
\begin{enumerate}\setlength{\itemsep=0pt}
 \item 
 $\alpha=0.07$, 
 $\beta=0.0$, 
 $\omega=-1$,
 $\vec{L}(0)=\vec{e}_3$.
 Como vimos en la hoja anterior, se tiende a un plano cercano al ecuador.
 
 \item 
 $\alpha=0.05$, 
 $\beta=0.015$, 
 $\omega=-1$,
 $\vec{L}(0)=\vec{e}_3$.
 Ejemplo con $\beta\ne0$ de que todavía se produce la localización a la larga cercana a un plano.
 
 \item 
 $\alpha=0.05$, 
 $\beta=0.015$, 
 $\omega=-1$,
 $\vec{L}(0)=\vec{e}_2$.
 Lo mismo que lo anterior pero ahora partiendo del ecuador.

 \item 
 $\alpha=0.0$, 
 $\beta=0.0$, 
 $\omega=-1$,
 $\vec{L}(0)=\vec{e}_3$.
 El caso de resonancia pura sin el ``rozamiento'' representado por $\alpha$ y $\beta$. Como vimos, se combinan dos rotaciones.

 \item 
 $\alpha=0.5$, 
 $\beta=0.01$, 
 $\omega=0$,
 $\vec{L}(0)=\vec{e}_2$.
 Si desconectamos el campo variable, $\alpha$ atrae a $\vec{L}$ a la vertical pero si $\beta$ es muy pequeña, la rotación se mantiene.

 \item 
 $\alpha=0.01$, 
 $\beta=0.5$, 
 $\omega=0$,
 $\vec{L}(0)=\vec{e}_2$.
 La situación recíproca a la anterior: cuando desconectamos el campo variable con  $\beta$ significativo, solo se aprecian unos pocos giros.

 \item 
 $\alpha=0.0$, 
 $\beta=0.0$, 
 $\omega=-2$,
 $\vec{L}(0)=\vec{e}_3$.
 En el caso en que no hay resonancia $\omega\ne -\omega_0$, la amplitud es pequeña.

 \item 
 $\alpha=0.0$, 
 $\beta=0.0$, 
 $\omega=-1$,
 $\vec{L}(0)=\vec{e}_1$.
 Para el caso de resonancia pura con $\alpha=0$ y $\beta=0$
 se tiene la solución $\vec{L}(t)=\big(\cos(\omega t), \sen(\omega t), 0\big)$ que parte de $\vec{e}_1$. 
\end{enumerate}

\

Hasta ahora sabemos de las oscilaciones colectivas de los protones con una frecuencia característica que depende del campo magnético principal pero nosotros queremos distinguir densidades, zonas con huecos de otras que no los tienen. Centrémonos en un plano $XY$ e imaginemos que añadimos un campo magnético  secundario (un \emph{gradiente de campo} \cite[III.A]{HiLe}) que crece a lo largo del eje~$OY$, lo cual se consigue en la práctica con electroimanes móviles. Entonces los protones que están más arriba ($y$ mayor) oscilarán más rápido y los de abajo más lentos. Desde fuera detectaremos una superposición de estas ocilaciones en forma de radiofrecuencia. Como habrás visto en más de un curso, el análisis de Fourier, concretamente la transformada de Fourier, da el contenido que tiene una señal de cada frecuencia por tanto podemos distinguir cuántos protones hay en cada horizontal. Si $\rho(x,y)$ es la densidad de protones en la ``rodaja'' que estamos examinando, sabremos $\int_{-\infty}^\infty \rho(x,y)\; dx$ para cada $y$. En  
\cite{chamizo_oc} intenté dar una versión para el gran público de esta explicación. Te pasaré un par de páginas de una versión antigua (de nuevo, si se me olvida, recuérdamelo). También puedes ver un resumen rápido de la idea en 6:00--16:05 de la lección 25 del curso:
\[
 \text{\footnotesize\url{https://www.youtube.com/watch?v=35gfOtjRcic&list=PLH_k6d9j2lGUms4QZzokrkXGEoP6oxkYH}}
\]
En realidad todo esto es una simplificación, centrarnos en una rodaja requiere también un gradiente. Para más detalle te remito a las lecciones 20-26 del curso anterior. Son más de cuatro horas y por tanto solo me atrevo a recomendártelo si tienes mucho interés. 

\ejer
Relee el párrafo anterior hasta que lo entiendas, y con las dos páginas que te pasaré y los 10 minutos de la lección 25, reproduce la idea con tus palabras en unas líneas. 
\medskip

Si la variación del campo magnético se produce en otras direcciones\footnote{Los ruidos que se oyen cuando nos hacen una resonancia vienen del movimiento de los electroimanes secundarios para ajustar la dirección.}, tendremos información sobre otras integrales de la densidad. El modelo matemático es que los datos obtenidos (tratados con la transformada de Fourier) nos dan la función 
\[
P_{\theta}(t)=\int_{r_{\theta,t}}\rho
\qquad\text{con }\quad
r_{\theta,t}\equiv x\cos\theta+y\sen\theta=t.
\]
Esto es lo 
que se llama, a veces con una normalización diferente, \emph{transformada de Radon}.
Nota que $r_{\theta,t}$ es la recta a distancia~$t$ del origen con vector normal formando un ángulo~$\theta$ con el eje~$OX$. Tomar como ángulo el de la normal en vez del de la recta, es solo un convenio habitual. 

\ejer
Prueba que $r_{\theta,t}$ tiene estas propiedades. 
\medskip

El problema matemático que se plantea en la tomografía es obtener la función densidad~$\rho$ a partir de la transformada de Radon.  Lo resolvió J.~Radon en 1917 mucho antes de que nadie soñara con la posibilidad del uso de la tomografía en imágenes médicas (en \cite{radon} hay una reimpresión). 
Este mismo problema aparece también en las TAC (\emph{Tomografía Axial Computarizada}), que son el pariente barato de las resonancias y allí el modelo es mucho más sencillo porque no hay resonancias de espines sino simples radiografías. 



\ejer
Lee lo que escribí para un curso en 
\[
 \text{\small\url{https://matematicas.uam.es/~fernando.chamizo/asignaturas/model1415/fourier.pdf}}
\]
sobre tomografías y transformada de Radon (pp.\,12--13) y complétalo  con 
\cite[pp.\,75--77]{chamizo_mo}. Si tienes acceso a \cite[7.5]{StSh}, allí hay también información, que seguramente esté mejor escrita que la de mis documentos. 

\ejer
Motiva la transformada de Radon en relación con la tomografía (como en las líneas anteriores) y enuncia y prueba la fórmula de reconstrucción, que obtiene $\rho$ a partir de $P_\theta(t)$, sin preocuparte por temas de regularidad y convergencia. 
\medskip


Hay un caso de la reconstrucción tomográfica que es  particularmente sencillo. Aunque no es relevante en las imágenes médicas, sí tiene otras aplicaciones. Corresponde a la situación en la que la muestra tiene simetría radial, matemáticamente 
$\rho(x,y)=f\big(\sqrt{x^2+y^2}\big)$,
lo que implica que $P_\theta(t)$ no depende del ángulo y se tiene (si no lo ves claro, hazlo como ejercicio)
\[
P_{\theta}(t)=\int_{-\infty}^\infty
f\big(\sqrt{t^2+u^2}\big)\; du.
\]
A esta función $F(t)=P_{\theta}(t)$ se le llama \emph{transformada de Abel} de $f$. Resulta que en este caso la fórmula para la reconstrucción se reduce a 
\[
 f(r)=
 -\frac{1}{\pi}
 \int_r^\infty
 \frac{F'(t)}{\sqrt{t^2-r^2}}dt
\]
y su prueba es relativamente simple. 
Realmente el tratamiento original de la transformada de Radon dado por el propio Radon \cite{radon}, está más cercano al que ahora hacemos de la transformada de Abel. 

\ejer
Busca la prueba de la fórmula anterior (está por ejemplo en \cite{wiki_abelt}, lee el comentario que añado a continuación) y escríbela con todas las explicaciones que consideres convenientes, sin entrar en temas de regularidad y convergencia. 
\medskip

Un comentario al margen es que las integrales 
\[
 I_k
 =
 \int_a^1
 \frac{t^{2k+1}\; dt}{\sqrt{(1-t^2)(t^2-a^2)}}
 \qquad\text{con}\quad k\in\mathbb{N},\quad a\in\R^+
\]
se pueden calcular explícitamente con variable compleja. En particular, se tiene  $I_0=\pi/2$, $I_1=\pi (a^2+1)/4$. La primera integral se usa, con un cambio de variable, en la explicación de \cite{wiki_abelt}. Si no sabes cómo hallarla, puedes apelar a $\text{WolframAlpha}^{\text{\textregistered}}$ en tu redacción.
\medskip

El siguiente ejercicio yo creo que es ilustrativo pero si estás agobiado o te resulta difícil, te lo puedes saltar. A mi juicio, la única dificultad es que te asusten los cálculos al comprobar la reconstrucción. 

\ejer
Imagina un tubo de sección circular de radio~$1$ cuya densidad crece cuadráticamente desde~$0$  en el centro a $3$ en el borde. Escribe las fórmulas para $\rho$, $f$ y $F$. Comprueba que se verifica la fórmula de reconstrucción empleando los valores de $I_0$ e $I_1$ anteriores, que puedes dar por supuestos.



\begin{center}
\noindent\rule{0.95\textwidth}{1.1pt}
\end{center}
\nopagebreak


\paragraph{Tarea a entregar.}
Como en anteriores ocasiones, debes escribir un documento que combine los ejercicios anteriores. 
La primera parte, la de los cálculos numéricos, debería ir en el capítulo correspondiente a la hoja anterior (la ecuación de Bloch) del trabajo final, mientras que el resto debe estar separado en uno nuevo. En el proyecto original lo llamé ``Tomografía y la transformada de Radon'' pero puedes poner otro nombre si te parece más adecuado. 

En \url{http://matematicas.uam.es/~fernando.chamizo/supervision/TFG/enlaces.html} tienes una plantilla para incluir figuras y código en {\LaTeX}, como ya te indiqué. Si la extensión de tu trabajo  excediera mucho los límites recomendados, una posibilidad es relegar total o parcialmente los programas y las figuras a un apéndice, que el tribunal no tiene obligación de consultar. 

Te dejo libertad con respecto a la extensión. Si ves que queda muy largo, puedes minimizar o incluso eliminar lo que se refiere a la transformada de Abel. 


\begin{center}
\noindent\rule{0.95\textwidth}{1.1pt}
\end{center}

\begin{thebibliography}{{Wik}20}

\bibitem[Cha04]{chamizo_mo}
F.~Chamizo.
\newblock Modelización {II} (un pase de modelos).
\newblock
  {\small\url{http://matematicas.uam.es/~fernando.chamizo/libreria/libreria.html}},
  2004.

\bibitem[Cha20]{chamizo_oc}
F.~Chamizo.
\newblock {\em La Matemática y la Ciencia oculta}.
\newblock UAM Ediciones, 2020.

\bibitem[HL83]{HiLe}
W.~S. Hinshaw and A.~H. Lent.
\newblock {An introduction to NMR imaging: From the Bloch equation to the
  imaging equation}.
\newblock {\em Proceedings of the IEEE}, 71(3):338--350, 1983.

\bibitem[Rad86]{radon}
J.~Radon.
\newblock On the determination of functions from their integral values along
  certain manifolds.
\newblock {\em IEEE Transactions on Medical Imaging}, 5:170--176, 1986.

\bibitem[SS03]{StSh}
E.~M. Stein and R.~Shakarchi.
\newblock {\em Fourier analysis}, volume~1 of {\em Princeton Lectures in
  Analysis}.
\newblock Princeton University Press, Princeton, NJ, 2003.
\newblock An introduction.

\bibitem[{Wik}20]{wiki_abelt}
{Wikipedia contributors}.
\newblock Abel transform --- {Wikipedia}{,} the free encyclopedia, 2020.
\newblock [Online; accessed 14-April-2021].

\end{thebibliography}

%\bibliography{bibtfg20}
%\bibliographystyle{alpha}

\end{document}

