diff --git a/AUTHORS.md b/AUTHORS.md new file mode 100644 index 0000000..545bd28 --- /dev/null +++ b/AUTHORS.md @@ -0,0 +1,3 @@ +**BLOCHUS** is developed by + +* Thomas Hiller (since 2019), thomas.hiller@bgr.de \ No newline at end of file diff --git a/BLOCHUS/BLOCHUS.m b/BLOCHUS/BLOCHUS.m index 85c84d3..c34c9d1 100644 --- a/BLOCHUS/BLOCHUS.m +++ b/BLOCHUS/BLOCHUS.m @@ -33,10 +33,10 @@ %------------- BEGIN CODE -------------- %% GUI 'header' info and default GUI settings -myui.version = '0.1.4'; -myui.date = '25.09.2020'; +myui.version = '0.1.5'; +myui.date = '30.06.2022'; myui.author = 'Thomas Hiller'; -myui.email = 'thomas.hiller[at]leibniz-liag.de'; +myui.email = 'thomas.hiller[at]bgr.de'; myui.fontsize = 9; myui.axfontsize = 11; diff --git a/BLOCHUS/BLOCHUS_loadDefaults.m b/BLOCHUS/BLOCHUS_loadDefaults.m index d7e96b0..29d460b 100644 --- a/BLOCHUS/BLOCHUS_loadDefaults.m +++ b/BLOCHUS/BLOCHUS_loadDefaults.m @@ -118,7 +118,7 @@ % pre-polarization B-field amplitude in units of [B0] init.PrePolFactor = [100 1e-6 1e6]; % polar angle of pre-polarization field [deg] -init.PrePolTheta = [90 0 360]; +init.PrePolTheta = [90 1e-3 360]; % azimuthal angle of pre-polarization field [deg] init.PrePolPhi = [0 0 360]; % pre-polarization B-field switch amplitude in units of [B0] diff --git a/CHANGELOG.md b/CHANGELOG.md index 99495ae..9f29fdf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,16 @@ # Changelog +## [0.1.5] - 2022-06-30 + +### Added +- custom pre-polarization current switch-off ramps are now supported via an interpolation routine (so far only via script and in MATLABTM) + +### Changed +- now multiple angle-axis-pairs are allowed as input for function `getRotationMatrixFromAngleandAxis` + +### Fixed +- fixed an issue in `getRotationMatrixFromVectors` for very small angles between vectors A and B (also edge case like parallel and anti-parallel vectors are now considered) + ## [0.1.4] - 2020-09-25 ### Added @@ -43,6 +54,7 @@ Initial Version +[0.1.5]: https://github.com/ThoHiller/nmr-blochus/compare/v0.1.4...v0.1.5 [0.1.4]: https://github.com/ThoHiller/nmr-blochus/compare/v0.1.3...v0.1.4 [0.1.3]: https://github.com/ThoHiller/nmr-blochus/compare/v0.1.2...v0.1.3 [0.1.2]: https://github.com/ThoHiller/nmr-blochus/compare/v.0.1.1...v0.1.2 diff --git a/README.md b/README.md index 2f909b4..0f3618c 100644 --- a/README.md +++ b/README.md @@ -48,7 +48,7 @@ In order to work properly you need to meet the following requirements: #### Operating System -I tested it successfully under Windows 7 (64bit) and 10 (64bit) with Matlab R2017a and newer. Always with the latest version of the GUI Layout Toolbox (current version is v2.3.4) +I tested it successfully under Windows 7 (64bit) and 10 (64bit) with Matlab R2017a and newer. Always with the latest version of the GUI Layout Toolbox (current version is v2.3.5) **NOTE:** So far I did not test anything on Linux or a Mac. If you get it to work on either of the two systems (which it basically should I guess) please let me know. @@ -63,7 +63,7 @@ I tested it successfully under Windows 7 (64bit) and 10 (64bit) with Matlab R201 ### Usage -1. By executing the start scripts (see above) +1. By executing the start script (see above) 2. Simply type `BLOCHUS` on the MATLABTM prompt (make sure the `blochus` folder is on the MATLABTM path) 3. Check the example scripts for the usage of the core functions without the GUI (inside the `scripts` folder) @@ -87,14 +87,15 @@ In no particular order and without guarantee that it will ever happen :-) : ### Cite as If you use BLOCHUS for your research, please cite it as: -Thomas Hiller. (2020, September 25). ThoHiller/nmr-blochus: v0.1.4 (Version v0.1.4). Zenodo. [https://doi.org/10.5281/zenodo.3829464] +Thomas Hiller. (2020, September 25). ThoHiller/nmr-blochus: v0.1.5 (Version v0.1.5). Zenodo. [https://doi.org/10.5281/zenodo.3829464] Note: Even though the version number might change due to updates, this DOI is permanent (represents all versions) and always links to the latest version. ### References - -1. Hiller, T., Dlugosch, R. and Müller-Petke, M., "Utilizing pre-polarization to enhance SNMR signals - effect of imperfect switch-off", Geophysical Journal International **222**(2), p.815-826, 2020, [DOI](https://doi.org/10.1093/gji/ggaa216) +1. Hiller, T., Costabel, S., Dlugosch, R. and Müller-Petke, M. "First Measurements of Surface Nuclear Magnetic Resonance Signals Without an Oscillating Excitation Pulse – Exploiting Non-Adiabatic Prepolarization Switch-Off", Geophysical Research Letters, **48**(23), e2021GL095371, 2021, [DOI](https://doi.org/10.1029/2021GL095371) +1. Hiller, T., Costabel, S., Radic, T., Dlugosch, R. and Müller-Petke, M. "Feasibility study on prepolarized surface nuclear magnetic resonance for soil moisture measurements", Vadose Zone Journal, **20**(5), e20138, 2021, [DOI](https://doi.org/10.1002/vzj2.20138) +2. Hiller, T., Dlugosch, R. and Müller-Petke, M., "Utilizing pre-polarization to enhance SNMR signals - effect of imperfect switch-off", Geophysical Journal International **222**(2), p.815-826, 2020, [DOI](https://doi.org/10.1093/gji/ggaa216) - - -

MATLAB is a registered trademark of The Mathworks, Inc.

diff --git a/doc/blochus/BLOCHUS/BLOCHUS.html b/doc/blochus/BLOCHUS/BLOCHUS.html index d568020..ace6cdc 100644 --- a/doc/blochus/BLOCHUS/BLOCHUS.html +++ b/doc/blochus/BLOCHUS/BLOCHUS.html @@ -105,10 +105,10 @@

SOURCE CODE ^%------------- BEGIN CODE -------------- 0034 0035 %% GUI 'header' info and default GUI settings -0036 myui.version = '0.1.4'; -0037 myui.date = '25.09.2020'; +0036 myui.version = '0.1.5'; +0037 myui.date = '30.06.2022'; 0038 myui.author = 'Thomas Hiller'; -0039 myui.email = 'thomas.hiller[at]leibniz-liag.de'; +0039 myui.email = 'thomas.hiller[at]bgr.de'; 0040 0041 myui.fontsize = 9; 0042 myui.axfontsize = 11; diff --git a/doc/blochus/BLOCHUS/BLOCHUS_loadDefaults.html b/doc/blochus/BLOCHUS/BLOCHUS_loadDefaults.html index 5ee46e9..260cef1 100644 --- a/doc/blochus/BLOCHUS/BLOCHUS_loadDefaults.html +++ b/doc/blochus/BLOCHUS/BLOCHUS_loadDefaults.html @@ -170,15 +170,15 @@

SOURCE CODE ^% initial magnetization [A/m] 0101 init.Minit = [1 0 0]; -0102 % Earth's magnetic field [T] default: 48µT min: 1fT max: 1T +0102 % Earth's magnetic field [T] default: 48µT min: 1fT max: 1T 0103 init.B0 = [48/1e6 1e-12 1]; 0104 % corresponding Larmor frequencies [Hz] 0105 init.Omega0 = [getOmega0(init.gamma,init.B0(1))/2/pi -1e9 1e9]; -0106 % T1 relaxation time [ms] default: 100ms min: 1µs max: 1000s +0106 % T1 relaxation time [ms] default: 100ms min: 1µs max: 1000s 0107 init.T1relax = [100 1e-3 1e6]; -0108 % T2 relaxation time [ms] default: 50ms min: 1µs max: 1000s +0108 % T2 relaxation time [ms] default: 50ms min: 1µs max: 1000s 0109 init.T2relax = [50 1e-3 1e6]; -0110 % simulation time [ms] default: 50ms min: 1µs max: 1000s +0110 % simulation time [ms] default: 50ms min: 1µs max: 1000s 0111 init.Tsim = [50 1e-3 1e6]; 0112 0113 % --- PRE-POLARIZATION --- @@ -189,14 +189,14 @@

SOURCE CODE ^% pre-polarization B-field amplitude in units of [B0] 0119 init.PrePolFactor = [100 1e-6 1e6]; 0120 % polar angle of pre-polarization field [deg] -0121 init.PrePolTheta = [90 0 360]; +0121 init.PrePolTheta = [90 1e-3 360]; 0122 % azimuthal angle of pre-polarization field [deg] 0123 init.PrePolPhi = [0 0 360]; 0124 % pre-polarization B-field switch amplitude in units of [B0] 0125 init.PrePolSwitchFactor = [1 1e-6 1e6]; -0126 % switch-off ramp time [ms] default: 1ms min: 1µs max: 1000s +0126 % switch-off ramp time [ms] default: 1ms min: 1µs max: 1000s 0127 init.PrePolTramp = [1 1e-6 1e6]; -0128 % switch-off ramp slope time [ms] default: 0.1ms min: 1µs max: 1000s +0128 % switch-off ramp slope time [ms] default: 0.1ms min: 1µs max: 1000s 0129 init.PrePolTslope = [0.1 1e-6 1e6]; 0130 0131 % --- PULSE --- @@ -208,7 +208,7 @@

SOURCE CODE ^'circular'; 0138 % relaxation during pulse [0/1] 0139 init.PulseRDP = 0; -0140 % pulse length [ms] default: 20ms min: 1µs max: 1000s +0140 % pulse length [ms] default: 20ms min: 1µs max: 1000s 0141 init.PulseTtau = [20 1e-6 1e6]; 0142 % pulse amplitude in units of [B0] default: 0.0061163 min: 1e-12 max: 1e6 0143 init.PulseB1Factor = [0.0061163 1e-12 1e6]; diff --git a/doc/blochus/functions/blochsim/fcn_BLOCHUS_ode.html b/doc/blochus/functions/blochsim/fcn_BLOCHUS_ode.html index 20ecc7b..b295981 100644 --- a/doc/blochus/functions/blochsim/fcn_BLOCHUS_ode.html +++ b/doc/blochus/functions/blochsim/fcn_BLOCHUS_ode.html @@ -140,7 +140,7 @@

SOURCE CODE ^switch param.type 0059 case 'std' -0060 B = B0*zunit; +0060 B = B0*zunit; 0061 % dM/dt 0062 dM = gamma*cross(m,B) - ( (m(1)*xunit+m(2)*yunit) / T2 ) - ( ((m(3)-M0(3))*zunit) / T1 ); 0063 diff --git a/doc/blochus/functions/blochsim/getMIDI_Tx.html b/doc/blochus/functions/blochsim/getMIDI_Tx.html index fdf113d..2f98cdc 100644 --- a/doc/blochus/functions/blochsim/getMIDI_Tx.html +++ b/doc/blochus/functions/blochsim/getMIDI_Tx.html @@ -317,7 +317,7 @@

SOURCE CODE ^% now create an artificial y-component by shifting the original signal -0249 % (x-component) 90° (quarter of a full period) +0249 % (x-component) 90° (quarter of a full period) 0250 % get the length of a quarter period 0251 t90 = dtL/4; 0252 % cut out these first samples and put them to the end of the original signal @@ -330,7 +330,7 @@

SOURCE CODE ^% pulse this is already taken care of during assembling) 0260 if strcmp(param.Tx,'MIDI_OR') && isfield(param,'PulseAxis') 0261 % depending on the frequency, there is a different amount of time -0262 % steps for the "90° phase shift" +0262 % steps for the "90° phase shift" 0263 shiftind = sum(t<t90); 0264 % now shift both components by the necessary amount of samples 0265 switch param.PulseAxis diff --git a/doc/blochus/functions/blochsim/getPulseTimeSeries.html b/doc/blochus/functions/blochsim/getPulseTimeSeries.html index 81ffcac..dfc8ef1 100644 --- a/doc/blochus/functions/blochsim/getPulseTimeSeries.html +++ b/doc/blochus/functions/blochsim/getPulseTimeSeries.html @@ -214,50 +214,53 @@

SOURCE CODE ^if Imod.useQ && Imod.Q > 0 0125 % get line (band) width -> f_L/Q (simple bandwidth for bandpass) 0126 Lwidth = abs(param.omega0/2/pi) / Imod.Q; -0127 % apply "Breit-Wigner" formula (here already normalized to 1 by -0128 % multiplying with pi*Lwidth) -0129 L = 1 ./ ( ((df+Imod.Qdf).^2 ./ Lwidth.^2) + 1 ); -0130 I = I.*L; -0131 end -0132 -0133 % get instantaneous phase due to frequency modulation -0134 theta = getPulsePhase(t,df,param,0); -0135 -0136 switch PulsePolarization -0137 case 'circular' -0138 Bout = I.*Amp.*[cos(theta + phi + phi_ax + phi_ref) ... -0139 sin(theta + phi + phi_ax + phi_ref)]; -0140 case 'linear' -0141 Bout = 2.*I.*Amp.*[cos(theta + phi + phi_ax + phi_ref) ... -0142 0.*t]; -0143 end -0144 end -0145 % if B-field values are NaN, set them to zero -0146 % this can happen due to the interpolation at the end of the MIDI-pulses -0147 Bout(isnan(Bout)) = 0; -0148 -0149 end -0150 -0151 %------------- END OF CODE -------------- -0152 -0153 %% License: -0154 % GNU GPLv3 -0155 % -0156 % BLOCHUS -0157 % Copyright (C) 2019 Thomas Hiller +0127 % apply Cauchy-Lorentz type formula (here already normalized to +0128 % 1 by multiplying with pi*Lwidth) +0129 % this is basically a Cauchy distribution PDF of the form: +0130 % PDF = 1/pi * ( bw / ((f-f0)^2 + bw^2) ) +0131 % tweaked with some algebra +0132 L = 1 ./ ( ((df+Imod.Qdf).^2 ./ Lwidth.^2) + 1 ); +0133 I = I.*L; +0134 end +0135 +0136 % get instantaneous phase due to frequency modulation +0137 theta = getPulsePhase(t,df,param,0); +0138 +0139 switch PulsePolarization +0140 case 'circular' +0141 Bout = I.*Amp.*[cos(theta + phi + phi_ax + phi_ref) ... +0142 sin(theta + phi + phi_ax + phi_ref)]; +0143 case 'linear' +0144 Bout = 2.*I.*Amp.*[cos(theta + phi + phi_ax + phi_ref) ... +0145 0.*t]; +0146 end +0147 end +0148 % if B-field values are NaN, set them to zero +0149 % this can happen due to the interpolation at the end of the MIDI-pulses +0150 Bout(isnan(Bout)) = 0; +0151 +0152 end +0153 +0154 %------------- END OF CODE -------------- +0155 +0156 %% License: +0157 % GNU GPLv3 0158 % -0159 % This program is free software: you can redistribute it and/or modify -0160 % it under the terms of the GNU General Public License as published by -0161 % the Free Software Foundation, either version 3 of the License, or -0162 % (at your option) any later version. -0163 % -0164 % This program is distributed in the hope that it will be useful, -0165 % but WITHOUT ANY WARRANTY; without even the implied warranty of -0166 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -0167 % GNU General Public License for more details. -0168 % -0169 % You should have received a copy of the GNU General Public License -0170 % along with this program. If not, see <https://www.gnu.org/licenses/>. +0159 % BLOCHUS +0160 % Copyright (C) 2019 Thomas Hiller +0161 % +0162 % This program is free software: you can redistribute it and/or modify +0163 % it under the terms of the GNU General Public License as published by +0164 % the Free Software Foundation, either version 3 of the License, or +0165 % (at your option) any later version. +0166 % +0167 % This program is distributed in the hope that it will be useful, +0168 % but WITHOUT ANY WARRANTY; without even the implied warranty of +0169 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +0170 % GNU General Public License for more details. +0171 % +0172 % You should have received a copy of the GNU General Public License +0173 % along with this program. If not, see <https://www.gnu.org/licenses/>.
Generated by m2html © 2005
\ No newline at end of file diff --git a/doc/blochus/functions/blochsim/getRampAmplitude.html b/doc/blochus/functions/blochsim/getRampAmplitude.html index c337653..f629ec8 100644 --- a/doc/blochus/functions/blochsim/getRampAmplitude.html +++ b/doc/blochus/functions/blochsim/getRampAmplitude.html @@ -140,77 +140,80 @@

SOURCE CODE ^case 'lin' % linear 0060 Bp = Bmax.*(1-t./Tramp); -0061 end -0062 -0063 return -0064 % --- +0061 case 'custom' +0062 idata = param.interp; +0063 Bp = Bmax.*interp1(idata.t,idata.I,t); +0064 end 0065 -0066 % --- -0067 function Bp = getLinExpAmp(Bmax,Bstar,T,t) -0068 % linear + exponential ramp after: -0069 % Conradi et al., 2017, Journal of Magnetic Resonance 281, p.241-245 -0070 % https://doi.org/10.1016/j.jmr.2017.06.001 -0071 -0072 if numel(t)>1 -0073 % linear part -0074 Bplin = (-Bmax/T)*t + Bmax; -0075 % exponential part -0076 Bpexp = exp(-t /(Bstar*T/Bmax)); -0077 % find change -0078 index = find(abs(Bplin-Bstar)==min(abs(Bplin-Bstar)),1,'first'); -0079 % merge the lin- and exp-part and scale the amplitude of the exp-part -0080 % to that of the lin-part at the switch-over time t(index) -0081 scale_point = Bplin(index)/Bpexp(index); -0082 % in case something goes south due to very small numbers set the -0083 % amplitude to 0 -0084 if isinf(scale_point) || isnan(scale_point) -0085 scale_point = 0; -0086 end -0087 % the final amplitude vector -0088 Bp = [Bplin(1:index-1); scale_point * Bpexp(index:end)]; -0089 else -0090 % linear part -0091 Bplin = (-Bmax/T)*t + Bmax; -0092 % exponential part -0093 Bpexp = exp(-t /(Bstar*T/Bmax)); -0094 % Bstar time tstar -0095 tstar = (Bstar-Bmax)/(-Bmax/T); -0096 % amplitude at tstar for scaling -0097 Btstar = exp(-tstar /(Bstar*T/Bmax)); -0098 % apply -0099 if t<tstar -0100 Bp = Bplin; -0101 else -0102 Bp = (Bstar/Btstar) * Bpexp; -0103 end -0104 end -0105 -0106 % if due to division by "0" the value is NaN ... set it to 0 -0107 Bp(isnan(Bp)) = 0; +0066 return +0067 % --- +0068 +0069 % --- +0070 function Bp = getLinExpAmp(Bmax,Bstar,T,t) +0071 % linear + exponential ramp after: +0072 % Conradi et al., 2017, Journal of Magnetic Resonance 281, p.241-245 +0073 % https://doi.org/10.1016/j.jmr.2017.06.001 +0074 +0075 if numel(t)>1 +0076 % linear part +0077 Bplin = (-Bmax/T)*t + Bmax; +0078 % exponential part +0079 Bpexp = exp(-t /(Bstar*T/Bmax)); +0080 % find change +0081 index = find(abs(Bplin-Bstar)==min(abs(Bplin-Bstar)),1,'first'); +0082 % merge the lin- and exp-part and scale the amplitude of the exp-part +0083 % to that of the lin-part at the switch-over time t(index) +0084 scale_point = Bplin(index)/Bpexp(index); +0085 % in case something goes south due to very small numbers set the +0086 % amplitude to 0 +0087 if isinf(scale_point) || isnan(scale_point) +0088 scale_point = 0; +0089 end +0090 % the final amplitude vector +0091 Bp = [Bplin(1:index-1); scale_point * Bpexp(index:end)]; +0092 else +0093 % linear part +0094 Bplin = (-Bmax/T)*t + Bmax; +0095 % exponential part +0096 Bpexp = exp(-t /(Bstar*T/Bmax)); +0097 % Bstar time tstar +0098 tstar = (Bstar-Bmax)/(-Bmax/T); +0099 % amplitude at tstar for scaling +0100 Btstar = exp(-tstar /(Bstar*T/Bmax)); +0101 % apply +0102 if t<tstar +0103 Bp = Bplin; +0104 else +0105 Bp = (Bstar/Btstar) * Bpexp; +0106 end +0107 end 0108 -0109 -0110 return +0109 % if due to division by "0" the value is NaN ... set it to 0 +0110 Bp(isnan(Bp)) = 0; 0111 -0112 %------------- END OF CODE -------------- -0113 -0114 %% License: -0115 % GNU GPLv3 -0116 % -0117 % BLOCHUS -0118 % Copyright (C) 2019 Thomas Hiller +0112 +0113 return +0114 +0115 %------------- END OF CODE -------------- +0116 +0117 %% License: +0118 % GNU GPLv3 0119 % -0120 % This program is free software: you can redistribute it and/or modify -0121 % it under the terms of the GNU General Public License as published by -0122 % the Free Software Foundation, either version 3 of the License, or -0123 % (at your option) any later version. -0124 % -0125 % This program is distributed in the hope that it will be useful, -0126 % but WITHOUT ANY WARRANTY; without even the implied warranty of -0127 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -0128 % GNU General Public License for more details. -0129 % -0130 % You should have received a copy of the GNU General Public License -0131 % along with this program. If not, see <https://www.gnu.org/licenses/>. +0120 % BLOCHUS +0121 % Copyright (C) 2019 Thomas Hiller +0122 % +0123 % This program is free software: you can redistribute it and/or modify +0124 % it under the terms of the GNU General Public License as published by +0125 % the Free Software Foundation, either version 3 of the License, or +0126 % (at your option) any later version. +0127 % +0128 % This program is distributed in the hope that it will be useful, +0129 % but WITHOUT ANY WARRANTY; without even the implied warranty of +0130 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +0131 % GNU General Public License for more details. +0132 % +0133 % You should have received a copy of the GNU General Public License +0134 % along with this program. If not, see <https://www.gnu.org/licenses/>.
Generated by m2html © 2005
\ No newline at end of file diff --git a/doc/blochus/functions/blochsim/getRotationMatrixFromAngleandAxis.html b/doc/blochus/functions/blochsim/getRotationMatrixFromAngleandAxis.html index bb6b44e..518e565 100644 --- a/doc/blochus/functions/blochsim/getRotationMatrixFromAngleandAxis.html +++ b/doc/blochus/functions/blochsim/getRotationMatrixFromAngleandAxis.html @@ -33,11 +33,11 @@

DESCRIPTION ^SOURCE CODE ^% getRotationMatrixFromAngleandAxis(phi,n) 0007 % 0008 % Inputs: -0009 % phi - rotation angle [rad] -0010 % n - rotation axis vector [x y z] +0009 % phi - rotation angle [rad]; size Nx1 +0010 % n - rotation axis vector [x y z]; size Nx3 0011 % 0012 % Outputs: -0013 % R - 3x3 rotation matrix +0013 % R - 3x3xN rotation matrix 0014 % 0015 % Example: 0016 % R = getRotationMatrixFromAngleandAxis(pi,[0 0 1]') @@ -110,44 +110,80 @@

SOURCE CODE ^%------------- BEGIN CODE -------------- 0037 -0038 % make "n" a unit vector -0039 n = n./norm(n); -0040 % get the individual components -0041 nx = n(1); -0042 ny = n(2); -0043 nz = n(3); -0044 % matrix terms needed -0045 omcos = 1-cos(phi); -0046 cosp = cos(phi); -0047 sinp = sin(phi); -0048 -0049 % assemble rotation matrix R -0050 R = [nx*nx*omcos + cosp nx*ny*omcos - nz*sinp nx*nz*omcos + ny*sinp; ... -0051 ny*nx*omcos + nz*sinp ny*ny*omcos + cosp ny*nz*omcos - nx*sinp; ... -0052 nz*nx*omcos - ny*sinp nz*ny*omcos + nx*sinp nz*nz*omcos + cosp]; -0053 -0054 return -0055 -0056 %------------- END OF CODE -------------- -0057 -0058 %% License: -0059 % GNU GPLv3 -0060 % -0061 % BLOCHUS -0062 % Copyright (C) 2019 Thomas Hiller -0063 % -0064 % This program is free software: you can redistribute it and/or modify -0065 % it under the terms of the GNU General Public License as published by -0066 % the Free Software Foundation, either version 3 of the License, or -0067 % (at your option) any later version. -0068 % -0069 % This program is distributed in the hope that it will be useful, -0070 % but WITHOUT ANY WARRANTY; without even the implied warranty of -0071 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -0072 % GNU General Public License for more details. -0073 % -0074 % You should have received a copy of the GNU General Public License -0075 % along with this program. If not, see <https://www.gnu.org/licenses/>. +0038 % for only one axis-angle pair +0039 if numel(phi) == 1 +0040 % make "n" a unit vector +0041 n = n./norm(n); +0042 % get the individual components +0043 nx = n(1); +0044 ny = n(2); +0045 nz = n(3); +0046 % matrix terms needed +0047 omcos = 1-cos(phi); +0048 cosp = cos(phi); +0049 sinp = sin(phi); +0050 +0051 % assemble rotation matrix R +0052 R(1,1) = nx*nx*omcos + cosp; +0053 R(1,2) = nx*ny*omcos - nz*sinp; +0054 R(1,3) = nx*nz*omcos + ny*sinp; +0055 +0056 R(2,1) = ny*nx*omcos + nz*sinp; +0057 R(2,2) = ny*ny*omcos + cosp; +0058 R(2,3) = ny*nz*omcos - nx*sinp; +0059 +0060 R(3,1) = nz*nx*omcos - ny*sinp; +0061 R(3,2) = nz*ny*omcos + nx*sinp; +0062 R(3,3) = nz*nz*omcos + cosp; +0063 +0064 else % for multiple axes and angles +0065 +0066 % n should contain only unit vectors! +0067 % get the individual components +0068 nx = n(:,1); +0069 ny = n(:,2); +0070 nz = n(:,3); +0071 % matrix terms needed +0072 omcos = 1-cos(phi); +0073 cosp = cos(phi); +0074 sinp = sin(phi); +0075 +0076 % assemble rotation matrix R +0077 R(1,1,:) = nx.*nx.*omcos + cosp; +0078 R(1,2,:) = nx.*ny.*omcos - nz.*sinp; +0079 R(1,3,:) = nx.*nz.*omcos + ny.*sinp; +0080 +0081 R(2,1,:) = ny.*nx.*omcos + nz.*sinp; +0082 R(2,2,:) = ny.*ny.*omcos + cosp; +0083 R(2,3,:) = ny.*nz.*omcos - nx.*sinp; +0084 +0085 R(3,1,:) = nz.*nx.*omcos - ny.*sinp; +0086 R(3,2,:) = nz.*ny.*omcos + nx.*sinp; +0087 R(3,3,:) = nz.*nz.*omcos + cosp; +0088 end +0089 +0090 return +0091 +0092 %------------- END OF CODE -------------- +0093 +0094 %% License: +0095 % GNU GPLv3 +0096 % +0097 % BLOCHUS +0098 % Copyright (C) 2019 Thomas Hiller +0099 % +0100 % This program is free software: you can redistribute it and/or modify +0101 % it under the terms of the GNU General Public License as published by +0102 % the Free Software Foundation, either version 3 of the License, or +0103 % (at your option) any later version. +0104 % +0105 % This program is distributed in the hope that it will be useful, +0106 % but WITHOUT ANY WARRANTY; without even the implied warranty of +0107 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +0108 % GNU General Public License for more details. +0109 % +0110 % You should have received a copy of the GNU General Public License +0111 % along with this program. If not, see <https://www.gnu.org/licenses/>.
Generated by m2html © 2005
\ No newline at end of file diff --git a/doc/blochus/functions/blochsim/getRotationMatrixFromVectors.html b/doc/blochus/functions/blochsim/getRotationMatrixFromVectors.html index 87c87e3..ca3cb50 100644 --- a/doc/blochus/functions/blochsim/getRotationMatrixFromVectors.html +++ b/doc/blochus/functions/blochsim/getRotationMatrixFromVectors.html @@ -113,46 +113,48 @@

SOURCE CODE ^% cross product of both vectors 0041 c = cross(A,B); -0042 % check if A and B are parallel / antiparallel -0043 if sum(c)==0 -0044 % check if A == B (parallel) -0045 if all(A==B) -0046 % in that case the rotation matrix is obviously identity -0047 R = eye(3); -0048 else % A == -B (antiparallel) -0049 R = -eye(3); -0050 end -0051 else -0052 % skew-symmetric cross-product -0053 ssc = [ 0 -c(3) c(2); -0054 c(3) 0 -c(1); -0055 -c(2) c(1) 0 ]; -0056 % rotation matrix R -0057 R = eye(3) + ssc + ssc^2*(1-dot(A,B))/(norm(c))^2; -0058 end -0059 -0060 return +0042 % angle between A and B +0043 alpha = acos(dot(A,B)); +0044 % check if A and B are parallel / antiparallel +0045 if abs(sum(c))<1e-128 && (alpha==0 || alpha==pi) +0046 % check if A == B (alpha=0 -> parallel) +0047 if alpha==0 +0048 % in that case the rotation matrix is obviously identity +0049 R = eye(3); +0050 else % A == -B (antiparallel) +0051 R = -eye(3); +0052 end +0053 else +0054 % skew-symmetric cross-product +0055 ssc = [ 0 -c(3) c(2); +0056 c(3) 0 -c(1); +0057 -c(2) c(1) 0 ]; +0058 % rotation matrix R +0059 R = eye(3) + ssc + (ssc/norm(c))^2*(1-dot(A,B)); +0060 end 0061 -0062 %------------- END OF CODE -------------- +0062 return 0063 -0064 %% License: -0065 % GNU GPLv3 -0066 % -0067 % BLOCHUS -0068 % Copyright (C) 2019 Thomas Hiller -0069 % -0070 % This program is free software: you can redistribute it and/or modify -0071 % it under the terms of the GNU General Public License as published by -0072 % the Free Software Foundation, either version 3 of the License, or -0073 % (at your option) any later version. -0074 % -0075 % This program is distributed in the hope that it will be useful, -0076 % but WITHOUT ANY WARRANTY; without even the implied warranty of -0077 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -0078 % GNU General Public License for more details. -0079 % -0080 % You should have received a copy of the GNU General Public License -0081 % along with this program. If not, see <https://www.gnu.org/licenses/>. +0064 %------------- END OF CODE -------------- +0065 +0066 %% License: +0067 % GNU GPLv3 +0068 % +0069 % BLOCHUS +0070 % Copyright (C) 2019 Thomas Hiller +0071 % +0072 % This program is free software: you can redistribute it and/or modify +0073 % it under the terms of the GNU General Public License as published by +0074 % the Free Software Foundation, either version 3 of the License, or +0075 % (at your option) any later version. +0076 % +0077 % This program is distributed in the hope that it will be useful, +0078 % but WITHOUT ANY WARRANTY; without even the implied warranty of +0079 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +0080 % GNU General Public License for more details. +0081 % +0082 % You should have received a copy of the GNU General Public License +0083 % along with this program. If not, see <https://www.gnu.org/licenses/>.
Generated by m2html © 2005
\ No newline at end of file diff --git a/doc/blochus/functions/interface/getRampParameters.html b/doc/blochus/functions/interface/getRampParameters.html index 45e3f83..10984c0 100644 --- a/doc/blochus/functions/interface/getRampParameters.html +++ b/doc/blochus/functions/interface/getRampParameters.html @@ -132,13 +132,13 @@

SOURCE CODE ^% now adjust the direction of the pre-polarization field 0063 % 1.) rotation by angle theta [deg] around y-axis for the z-unit vector -0064 % this means an angle of 90° around the y-axis will turn the z-unit into +0064 % this means an angle of 90° around the y-axis will turn the z-unit into 0065 % the x-unit vector 0066 RM = getRotationMatrixFromAngleandAxis(deg2rad(data.prepol.Theta),[0 1 0]); 0067 orient = RM * zunit; 0068 % now this new orientation vector gets rotated by angle phi [deg] around 0069 % the z-axis for the x-unit vector -0070 % this means an angle of 90° around the z-axis will turn the x-unit into +0070 % this means an angle of 90° around the z-axis will turn the x-unit into 0071 % the y-unit vector 0072 RM = getRotationMatrixFromAngleandAxis(deg2rad(data.prepol.Phi),[0 0 1]); 0073 orient = RM * orient; diff --git a/functions/blochsim/fcn_BLOCHUS_ode.m b/functions/blochsim/fcn_BLOCHUS_ode.m index 04e557e..b39f6f4 100644 --- a/functions/blochsim/fcn_BLOCHUS_ode.m +++ b/functions/blochsim/fcn_BLOCHUS_ode.m @@ -57,7 +57,7 @@ switch param.type case 'std' - B = B0*zunit; + B = B0*zunit; % dM/dt dM = gamma*cross(m,B) - ( (m(1)*xunit+m(2)*yunit) / T2 ) - ( ((m(3)-M0(3))*zunit) / T1 ); diff --git a/functions/blochsim/getPulseTimeSeries.m b/functions/blochsim/getPulseTimeSeries.m index 1739567..e63441d 100644 --- a/functions/blochsim/getPulseTimeSeries.m +++ b/functions/blochsim/getPulseTimeSeries.m @@ -124,8 +124,11 @@ if Imod.useQ && Imod.Q > 0 % get line (band) width -> f_L/Q (simple bandwidth for bandpass) Lwidth = abs(param.omega0/2/pi) / Imod.Q; - % apply "Breit-Wigner" formula (here already normalized to 1 by - % multiplying with pi*Lwidth) + % apply Cauchy-Lorentz type formula (here already normalized to + % 1 by multiplying with pi*Lwidth) + % this is basically a Cauchy distribution PDF of the form: + % PDF = 1/pi * ( bw / ((f-f0)^2 + bw^2) ) + % tweaked with some algebra L = 1 ./ ( ((df+Imod.Qdf).^2 ./ Lwidth.^2) + 1 ); I = I.*L; end diff --git a/functions/blochsim/getRampAmplitude.m b/functions/blochsim/getRampAmplitude.m index 7ba0ecf..29acc85 100644 --- a/functions/blochsim/getRampAmplitude.m +++ b/functions/blochsim/getRampAmplitude.m @@ -58,6 +58,9 @@ Bp = Bmax .* (0.5+(cos(pi*t./Tramp)./2)); case 'lin' % linear Bp = Bmax.*(1-t./Tramp); + case 'custom' + idata = param.interp; + Bp = Bmax.*interp1(idata.t,idata.I,t); end return diff --git a/functions/blochsim/getRotationMatrixFromAngleandAxis.m b/functions/blochsim/getRotationMatrixFromAngleandAxis.m index 793a597..97700d2 100644 --- a/functions/blochsim/getRotationMatrixFromAngleandAxis.m +++ b/functions/blochsim/getRotationMatrixFromAngleandAxis.m @@ -6,11 +6,11 @@ % getRotationMatrixFromAngleandAxis(phi,n) % % Inputs: -% phi - rotation angle [rad] -% n - rotation axis vector [x y z] +% phi - rotation angle [rad]; size Nx1 +% n - rotation axis vector [x y z]; size Nx3 % % Outputs: -% R - 3x3 rotation matrix +% R - 3x3xN rotation matrix % % Example: % R = getRotationMatrixFromAngleandAxis(pi,[0 0 1]') @@ -35,21 +35,57 @@ %------------- BEGIN CODE -------------- -% make "n" a unit vector -n = n./norm(n); -% get the individual components -nx = n(1); -ny = n(2); -nz = n(3); -% matrix terms needed -omcos = 1-cos(phi); -cosp = cos(phi); -sinp = sin(phi); +% for only one axis-angle pair +if numel(phi) == 1 + % make "n" a unit vector + n = n./norm(n); + % get the individual components + nx = n(1); + ny = n(2); + nz = n(3); + % matrix terms needed + omcos = 1-cos(phi); + cosp = cos(phi); + sinp = sin(phi); + + % assemble rotation matrix R + R(1,1) = nx*nx*omcos + cosp; + R(1,2) = nx*ny*omcos - nz*sinp; + R(1,3) = nx*nz*omcos + ny*sinp; + + R(2,1) = ny*nx*omcos + nz*sinp; + R(2,2) = ny*ny*omcos + cosp; + R(2,3) = ny*nz*omcos - nx*sinp; + + R(3,1) = nz*nx*omcos - ny*sinp; + R(3,2) = nz*ny*omcos + nx*sinp; + R(3,3) = nz*nz*omcos + cosp; -% assemble rotation matrix R -R = [nx*nx*omcos + cosp nx*ny*omcos - nz*sinp nx*nz*omcos + ny*sinp; ... - ny*nx*omcos + nz*sinp ny*ny*omcos + cosp ny*nz*omcos - nx*sinp; ... - nz*nx*omcos - ny*sinp nz*ny*omcos + nx*sinp nz*nz*omcos + cosp]; +else % for multiple axes and angles + + % n should contain only unit vectors! + % get the individual components + nx = n(:,1); + ny = n(:,2); + nz = n(:,3); + % matrix terms needed + omcos = 1-cos(phi); + cosp = cos(phi); + sinp = sin(phi); + + % assemble rotation matrix R + R(1,1,:) = nx.*nx.*omcos + cosp; + R(1,2,:) = nx.*ny.*omcos - nz.*sinp; + R(1,3,:) = nx.*nz.*omcos + ny.*sinp; + + R(2,1,:) = ny.*nx.*omcos + nz.*sinp; + R(2,2,:) = ny.*ny.*omcos + cosp; + R(2,3,:) = ny.*nz.*omcos - nx.*sinp; + + R(3,1,:) = nz.*nx.*omcos - ny.*sinp; + R(3,2,:) = nz.*ny.*omcos + nx.*sinp; + R(3,3,:) = nz.*nz.*omcos + cosp; +end return diff --git a/functions/blochsim/getRotationMatrixFromVectors.m b/functions/blochsim/getRotationMatrixFromVectors.m index d55f399..abb7d2d 100644 --- a/functions/blochsim/getRotationMatrixFromVectors.m +++ b/functions/blochsim/getRotationMatrixFromVectors.m @@ -39,10 +39,12 @@ B = B(:)./norm(B); % cross product of both vectors c = cross(A,B); +% angle between A and B +alpha = acos(dot(A,B)); % check if A and B are parallel / antiparallel -if sum(c)==0 - % check if A == B (parallel) - if all(A==B) +if abs(sum(c))<1e-128 && (alpha==0 || alpha==pi) + % check if A == B (alpha=0 -> parallel) + if alpha==0 % in that case the rotation matrix is obviously identity R = eye(3); else % A == -B (antiparallel) @@ -54,7 +56,7 @@ c(3) 0 -c(1); -c(2) c(1) 0 ]; % rotation matrix R - R = eye(3) + ssc + ssc^2*(1-dot(A,B))/(norm(c))^2; + R = eye(3) + ssc + (ssc/norm(c))^2*(1-dot(A,B)); end return diff --git a/pyBLOCHUS/README.md b/pyBLOCHUS/README.md index 3646d1f..b0d82cb 100644 --- a/pyBLOCHUS/README.md +++ b/pyBLOCHUS/README.md @@ -36,11 +36,11 @@ I tested it successfully under Windows 7 (64bit) and 10 (64bit) with [Anaconda]( 1. Makes sure that all required packages are installed and work properly 2. Download or `clone` the **BLOCHUS** repository and put the folder somewhere on your hard drive. -3. As I recommend uying an IDE like [Spyder](https://www.spyder-ide.org/), make sure that the **pyBLOCHUS** folder is on your PYTHONPATH (e.g. in [Spyder](https://www.spyder-ide.org/): Tools -> PYTHONPATH manager -> Add path) +3. As I recommend using an IDE like [Spyder](https://www.spyder-ide.org/), make sure that the **pyBLOCHUS** folder is on your PYTHONPATH (e.g. in [Spyder](https://www.spyder-ide.org/): Tools -> PYTHONPATH manager -> Add path) 4. Sometimes a [Spyder](https://www.spyder-ide.org/) restart is necessary tu update/sync the PYTHONPATH - - - ### Usage -Have a look in the [example](pyBLOCHUS/examples) scripts for the usage of the different modules. +Have a look in the [example](examples) scripts for the usage of the different modules. diff --git a/pyBLOCHUS/pyBLOCHUS/blochus_pulse.py b/pyBLOCHUS/pyBLOCHUS/blochus_pulse.py index e7b3a31..56b0443 100644 --- a/pyBLOCHUS/pyBLOCHUS/blochus_pulse.py +++ b/pyBLOCHUS/pyBLOCHUS/blochus_pulse.py @@ -484,8 +484,11 @@ def get_pulse_amplitude(self, time): # get line (band) width # -> f_L/Q (simple bandwidth for bandpass) lwidth = abs(self.larmor_f)/self.imod.qual_fac - # apply "Breit-Wigner" formula (here already normalized to 1 by - # multiplying with pi*lwidth) + # apply Cauchy-Lorentz type formula (here already normalized to + # 1 by multiplying with pi*lwidth) + # this is basically a Cauchy distribution PDF of the form: + # PDF = 1/pi * ( bw / ((f-f0)^2 + bw^2) ) + # tweaked with some algebra ltmp = 1.0 / (((freq_offset+self.imod.qual_fac_df)**2.0 / lwidth**2)+1.0) current = current*ltmp diff --git a/scripts/example3.m b/scripts/example3.m index 2ce93a4..8a5e4cd 100644 --- a/scripts/example3.m +++ b/scripts/example3.m @@ -245,7 +245,7 @@ for b = 1:numel(Bstar) for r = 1:numel(theta) - % initial orientation is towards y-axis + % initial orientation is towards x-axis orient = getRotationMatrixFromAngleandAxis(deg2rad(theta(r)),[0 1 0])*zunit; % initial magnetization in the direction of B0+Bp Morient = orient*rampparam.Bmax + odeparam.B0*zunit; diff --git a/scripts/example4.m b/scripts/example4.m index 3bf7306..5754d42 100644 --- a/scripts/example4.m +++ b/scripts/example4.m @@ -140,7 +140,7 @@ % switch-over amplitude for the "linexp" ramp (factor*B0) [T] rampparam.Bstar = SwitchFactor*odeparam.B0; - % initial orientation is towards y-axis + % initial orientation is towards x-axis orient = getRotationMatrixFromAngleandAxis(deg2rad(theta(tt)),[0 1 0])*zunit; % initial magnetization in the direction of B0+Bp Morient = orient*rampparam.Bmax + odeparam.B0*zunit; diff --git a/scripts/example6.m b/scripts/example6.m new file mode 100644 index 0000000..8eb30a5 --- /dev/null +++ b/scripts/example6.m @@ -0,0 +1,198 @@ +% script: example6.m +% This script demonstrates how to calculate a lookup-table for an arbitrarily +% shaped pre-polarization switch-off ramp with BLOCHUS. +% +% The switch-off ramp is imported as 'time' and 'current' data and used +% inside BLOCHUS via interpolation +% +% See also: BLOCHUS +% Author: Thomas Hiller +% email: thomas.hiller[at]leibniz-liag.de +% License: GNU GPLv3 (at end) + +%------------- BEGIN CODE -------------- + +% clear data +clear variables; clc; + +% because this script takes approx. 30min, ask the user before continuing +answer = questdlg({'This script needs about 30 min. to run.',... + 'Do you want to continue?'},'Continue?','Yes','No','Yes'); +% handle response +switch answer + case 'Yes' + % nothing to do + case {'No',''} + return +end + +% output switches +save_data = true; +plot_data = true; + +% --- LOOKUP TABLE SETTINGS ----------------------------------------------- +% global lookup table settings (change as you like): +% three custom ramps +lookup_name = {'ramp_2ms','ramp_3ms','ramp_3ms_raw'}; + +% IMPORTANT NOTE: +% This very coarse (and basically useless) Np=10 x Nt=10 grid takes for the +% 3 ramps roughly 30min! +% local lookup table settings (careful with the discretization) +Np = 10; +Nt = 10; +% pre-polarization factor [B0] +PPfac = logspace(-1,4,Np); +% angle between Bp and B0 [deg] +theta = linspace(0,180,Nt); +% ------------------------------------------------------------------------- + +% basic settings +zunit = [0 0 1]'; % z unit vector +% ODE solver error tolerance +tol = 1e-9; +% ODE solver options +options = odeset('RelTol',tol,'AbsTol',[tol tol tol]); + +% hydrogen proton +nucleus = '1H'; + +% parameter needed for the ODE solver +% simulation type [string] +odeparam.type = 'prepol'; +% equilibrium magnetization [A/m] +odeparam.M0 = zunit; +% primary (Earth's) magnetic field B0 [T] +odeparam.B0 = 4.8e-5; +% relaxation time T1 [s] +odeparam.T1 = 2; +% relaxation time T2 [s] +odeparam.T2 = 1; +% gyromagnetic ratio [rad/s/T] +odeparam.gamma = getGyroRatio(nucleus); + +% pre-polarization switch-off ramp parameter +% relaxation during switch-off [0/1] +odeparam.RDS = 0; + +% because Minit is not normalized to 1, M0 needs to be adjusted +odeparam.M0 = odeparam.B0*zunit; + +data = cell(1,numel(lookup_name)); +for i1 = 1:numel(lookup_name) + % load ramp file + dataR = load([lookup_name{i1},'.mat']); + dataR = dataR.data; + Tramp = dataR.t(end); + + % create the output data matrix + % adiabatic quality p + p = zeros(numel(theta),numel(PPfac)); + % final magnetization vector + Mfinal = zeros(numel(theta),numel(PPfac),3); + % loop over the local parameters + for tt = 1:numel(theta) + for pp = 1:numel(PPfac) + + % switch-off ramp type [string] + rampparam.ramp = 'custom'; + % gyromagnetic ratio [rad/s/T] + rampparam.gamma = odeparam.gamma; + % primary (Earth's) magnetic field B0 [T] + rampparam.B0 = odeparam.B0; + % amplitude of pre-polarization field [T] + rampparam.Bmax = PPfac(pp)*odeparam.B0; + % switch-off ramp time [s] + rampparam.Tramp = Tramp; + % ramp interpolation data + rampparam.interp.t = dataR.t; + rampparam.interp.I = dataR.I; + + % switch-off slope time is NOT NEEDED [s] + rampparam.Tslope = Tramp; + % switch-over amplitude is NOT NEEDED [T] + rampparam.Bstar = odeparam.B0; + + % initial orientation is towards x-axis + orient = getRotationMatrixFromAngleandAxis(deg2rad(theta(tt)),[0 1 0])*zunit; + % initial magnetization in the direction of B0+Bp + Morient = orient*rampparam.Bmax + odeparam.B0*zunit; + % Minit does not get normalized in this case + Minit = Morient; + % orientation of the pre-polarization pulse axis (x-axis) + odeparam.orient = orient; + + % add the ramp parameter to the ode parameter struct + odeparam.rampparam = rampparam; + + % ODE solver call + % OUTPUT: time T and magnetization M + [T,M] = ode45(@(t,m) fcn_BLOCHUS_ode(t,m,odeparam),rampparam.interp.t,Minit,options); + + % normalized magnetization vector at end of switch-off ramp + MMn = M(end,:)./norm(M(end,:)); + % "adiabatic quality" p of the switch-off ramp + % -> orientation of M with respect to z_unit + p(tt,pp) = dot(MMn(end,:),zunit)./norm(zunit); + % final orientation of M + Mfinal(tt,pp,1:3) = M(end,:); + + end + infostr1 = lookup_name{i1}; + infostr2 = ['theta=',num2str(theta(tt)),'°']; + disp([infostr1,' | ',infostr2]); + end + % save the local lookup table for this particular combination of + % ramp shape and ramp time + data{i1}.ramp = lookup_name{i1}; + data{i1}.PPfac = PPfac; + data{i1}.theta = theta; + data{i1}.odeparam = odeparam; + data{i1}.M = Mfinal; + data{i1}.p = p; +end + +% save lookup table +if save_data + save('example6_lookup_table.mat','data'); +end + +% plot data +if plot_data + [xx,yy] = meshgrid(linspace(0,180,Nt+1),logspace(-1,4,Np+1)); + figure; + for i1 = 1:numel(lookup_name) + p = data{i1}.p; + subplot(1,numel(lookup_name),i1); + surf(xx,yy,ones(size(xx)),p'); view(2); + shading flat; + set(gca,'XLim',[0 180],'XTick',linspace(0,180,5),'YScale','log',... + 'YLim',[0.1 1e4],'YTick',logspace(-1,4,6),'FontSize',12); + xlabel('\theta [deg]'); + ylabel('Bp [B0]'); + tstr = lookup_name{i1}; + tstr = strrep(tstr,'_',' '); + title(tstr); + end +end + +%------------- END OF CODE -------------- + +%% License: +% GNU GPLv3 +% +% BLOCHUS +% Copyright (C) 2019 Thomas Hiller +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% This program is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. If not, see . diff --git a/scripts/example6_ref.png b/scripts/example6_ref.png new file mode 100644 index 0000000..6daa3ff Binary files /dev/null and b/scripts/example6_ref.png differ diff --git a/scripts/ramp_2ms.mat b/scripts/ramp_2ms.mat new file mode 100644 index 0000000..946ceec Binary files /dev/null and b/scripts/ramp_2ms.mat differ diff --git a/scripts/ramp_3ms.mat b/scripts/ramp_3ms.mat new file mode 100644 index 0000000..5efb23f Binary files /dev/null and b/scripts/ramp_3ms.mat differ diff --git a/scripts/ramp_3ms_raw.mat b/scripts/ramp_3ms_raw.mat new file mode 100644 index 0000000..9ceaf2b Binary files /dev/null and b/scripts/ramp_3ms_raw.mat differ