update doc

README.md

@@ -8,8 +8,10 @@
 
 This package is inspired by [SYCOMORE](https://sycomore.readthedocs.io/) with the intention to use Automatic differentiation (with ForwardDiff.jl)
 
+Take a look at the [documentation](https://atrotier.github.io/EPGsim.jl/dev/).
+
 ## TO DO list :
-- Delete states with amplitude < defined threshold
--  discrete / discrete 3D
+- Delete states with amplitude < defined threshold or number of states (implemented in epgDephasing ??)
+- discrete / discrete 3D
 - Diffusion
-- Magnetization transfer
+- Magnetization Transfer

docs/src/AD.md

@@ -111,6 +111,36 @@ lines!(ax,TE_vec,df)
 f
 ```
 
+## Differentiation along multiple variables
+If we want to obtain the derivation along T1 and T2 we need to change the EPG_MESE function. The function should take as input a vector containing T2 and T1 (here noted T2/T1) :
+```@example AD
+function MESE_EPG2(T2T1,TE,ETL,delta)
+  T2,T1 = T2T1
+  T = complex(eltype(T2))
+  E = EPGStates([T(0.0)],[T(0.0)],[T(1.0)])
+  echo_vec = Vector{Complex{eltype(T2)}}()
+
+  E = epgRotation(E,pi/2*delta, pi/2)
+  ##loop over refocusing-pulses
+  for i = 1:ETL
+    E = epgDephasing(E,1)
+    E = epgRelaxation(E,TE,T1,T2)
+    E  = epgRotation(E,pi*delta,0.0)
+    E  = epgDephasing(E,1)
+    push!(echo_vec,E.Fp[1])
+  end
+
+  return abs.(echo_vec)
+end
+
+j2 = ForwardDiff.jacobian(x -> MESE_EPG2(x,TE,ETL,deltaB1),[T2,T1])
+```
+Here we can see that the second column corresponding to T1 is equal to 0 which is expected for a MESE sequence and the derivative along T2 gives the same results :
+
+```@example AD
+j2[:,1] ≈ vec(j)
+```
+
 ## Reproducibility
 
 This page was generated with the following version of Julia:

docs/src/index.md

@@ -14,7 +14,7 @@ The principal aspect of this package was to make it compatible with Automatic Di
 !!! note
     EPGsim.jl is work in progress and in some parts not entirely optimized. The interface is susceptible to change between version
 
-## physics concept
+## EPG concept 
 Introduction to the physics concepts behing EPG as well as their usage can be found on the rad229 youtube channels by Brian Hargreaves and Daniel Ennis :
 - [Lecture-04A: Definition of the Extended Phase Graph Basis](https://www.youtube.com/watch?v=bskhnaoJVNY)
 - [Lecture-04B: Sequence Operations in the Extended Phase Graph Domain](https://www.youtube.com/watch?v=kToL-9ZTzCs)

docs/src/regular.md

@@ -54,6 +54,10 @@ E = epgDephasing(E,1)
 E = epgRotation(E,deg2rad(60),deg2rad(117))
 ```
 
+!!! note
+    Currently, all the EPGstates are stored and used for calculation. 
+    The states equal or really close to zero are not deleted
+
 # Accessing states
 States can seen directly as a vector :
 ```@example Regular

src/RegularEPG.jl

@@ -1,6 +1,23 @@
 export epgDephasing,epgRelaxation,epgRotation,rfRotation
 export EPGStates, getStates
 
+"""
+    EPGStates{T <: Real} 
+
+Stores the EPG states in 3 vectors Fp,Fn and Z.
+
+# Constructors :
+    EPGStates(Fp::Vector{Complex{S}},Fn::Vector{Complex{S}},Z::Vector{Complex{S}}) where {S <: Real}
+    EPGStates(Fp::T=0,Fn::T=0,Z::T=1) where T <: Number
+
+# Fields
+- `Fp::Vector{Complex{T}}`
+- `Fn::Vector{Complex{T}}`
+- `Z::Vector{Complex{T}}`
+
+# Related functions
+- `getStates(E::EPGStates)` : extract EPG states as matrix 3xN
+"""
 mutable struct EPGStates{T <: Real} 
   Fp::Vector{Complex{T}}
   Fn::Vector{Complex{T}}
@@ -18,7 +35,11 @@ mutable struct EPGStates{T <: Real}
   end
 end
 
+"""
+    getStates(E::EPGStates)
 
+Extract EPG states as matrix 3xN
+"""
 function EPGStates(Fp::T=0,Fn::T=0,Z::T=1) where T <: Number
   T2 = ComplexF64
 
@@ -33,8 +54,9 @@ end
     epgDephasing(E::EPGStates, n=1) where T
   
 shifts the transverse dephasing states `F` corresponding to n dephasing-cycles.
+n can be any integer
 """
-function epgDephasing(E::EPGStates, n=1)
+function epgDephasing(E::EPGStates, n::Int=1)
 
   if(abs(n)>1)
     for i in 1:abs(n)
@@ -88,9 +110,13 @@ function epgRelaxation(E::EPGStates,t,T1, T2)
 end
 
 """
-    rfRotation_AT(alpha, phi=0.)
+    rfRotation(alpha, phi=0.)
 
 returns the rotation matrix for a pulse with flip angle `alpha` and phase `phi`.
+
+  # Arguments
+* `alpha`  - flip angle (radian)
+* `phi=0.` - phase of the flip angle (radian)
 """
 function rfRotation(alpha, phi=0.)
   R = [ cos(alpha/2.)^2   exp(2*im*phi)*sin(alpha/2.)^2   -im*exp(im*phi)*sin(alpha);