1function grad=makeTrapezoid(channel, varargin)
2%makeTrapezoid Create a trapezoid gradient event.
3%
4% PURPOSE
5% Build a trapezoidal gradient event struct (rise / flat-top / fall)
6% for a given logical channel. The returned struct is consumed by
7% mr.Sequence/addBlock to add the gradient to a sequence.
8%
9% SIGNATURES
10% g = mr.makeTrapezoid(channel, ...) % uses mr.opts() defaults
11% g = mr.makeTrapezoid(channel, system, ...) % system as 2nd positional arg
12% g = mr.makeTrapezoid(channel, ..., 'system', system) % system as name/value
13% g = mr.makeTrapezoid(channel, ..., 'Duration', d, 'Area', a)
14% g = mr.makeTrapezoid(channel, ..., 'FlatTime', ft, 'FlatArea', fa)
15% g = mr.makeTrapezoid(channel, ..., 'FlatTime', ft, 'Amplitude', amp)
16% g = mr.makeTrapezoid(channel, ..., 'Area', a) % shortest possible timing
17%
18% Exactly one of 'Area', 'FlatArea', or 'Amplitude' must be supplied.
19% Timing is determined by 'FlatTime' if given, else 'Duration' if given,
20% else the shortest realizable timing for the requested 'Area'.
21% Parameter names are case-insensitive.
22%
23% INPUTS
24% channel [required] char, 'x'|'y'|'z'
25% system [optional] struct from mr.opts; defaults to mr.opts() if omitted
26% 'Duration' [name/value] double, total duration including ramps, seconds, >0
27% 'Area' [name/value] double, total gradient area including ramps, 1/m
28% 'FlatTime' [name/value] double, flat-top duration, seconds
29% 'FlatArea' [name/value] double, flat-top-only area, 1/m
30% 'Amplitude' [name/value] double, flat-top amplitude, Hz/m
31% 'maxGrad' [name/value] double, override system.maxGrad, Hz/m
32% 'maxSlew' [name/value] double, override system.maxSlew, Hz/m/s
33% 'riseTime' [name/value] double, force rise time, seconds
34% 'fallTime' [name/value] double, force fall time, seconds (requires riseTime)
35% 'delay' [name/value] double, pre-event delay, seconds, default 0
36%
37% OUTPUT
38% grad struct with fields:
39% .type char, always 'trap' (includes degenerate triangle, flatTime=0)
40% .channel char, 'x'|'y'|'z'
41% .amplitude double, flat-top amplitude, Hz/m
42% .riseTime double, ramp-up duration, seconds
43% .flatTime double, flat-top duration, seconds (0 for triangular)
44% .fallTime double, ramp-down duration, seconds
45% .area double, total area including ramps, 1/m
46% (= amplitude * (flatTime + riseTime/2 + fallTime/2))
47% .flatArea double, flat-top-only area, 1/m (= amplitude * flatTime)
48% .delay double, pre-event delay, seconds
49% .first double, gradient value at t=0, Hz/m (always 0 for trap)
50% .last double, gradient value at end, Hz/m (always 0 for trap)
51%
52% ERRORS
53% makeTrapezoid:invalidArguments
54% - 'fallTime' specified without 'riseTime'.
55% - Not exactly one of 'Area' / 'FlatArea' / 'Amplitude' supplied.
56% - 'FlatTime' supplied without 'FlatArea' or 'Amplitude'.
57% - Neither 'Area' nor 'Duration' supplied.
58% makeTrapezoid:invalidDuration
59% - Requested area cannot be realized within the requested duration
60% under maxGrad/maxSlew. Error message reports the minimum
61% achievable duration in microseconds.
62% makeTrapezoid:invalidAmplitude
63% - Computed amplitude exceeds maxGrad.
64% Assertion failure
65% - With explicit riseTime+duration: duration < riseTime+fallTime,
66% or computed amplitude exceeds maxGrad.
67%
68% NOTES
69% - All ramp/flat times are rounded up to system.gradRasterTime.
70% - Internal storage uses Hz/m and Hz/m/s regardless of the units
71% passed to mr.opts (mr.opts converts on input). Use mr.convert
72% if you need physical units (mT/m, T/m/s, etc.).
73% - Caches an inputParser in a persistent variable for performance;
74% no other global state.
75%
76% EXAMPLE
77% sys = mr.opts('MaxGrad', 30, 'GradUnit', 'mT/m', ...
78% 'MaxSlew', 170, 'SlewUnit', 'T/m/s');
79% Nx = 256; fov = 256e-3; deltak = 1/fov;
80% % Readout gradient with fixed flat-top area
81% gx = mr.makeTrapezoid('x', sys, 'FlatArea', Nx*deltak, 'FlatTime', 6.4e-3);
82% % Matching prephaser: half the (negative) area of the readout
83% gxPre = mr.makeTrapezoid('x', sys, 'Area', -gx.area/2, 'Duration', 1e-3);
84%
85% SEE ALSO
86% mr.opts, mr.makeExtendedTrapezoid, mr.makeArbitraryGrad,
87% mr.calcDuration, mr.Sequence/addBlock
89persistent parser
91if isempty(parser)
92 validChannels = {'x','y','z'};
93 parser = mr.aux.InputParserCompat;
94 parser.FunctionName = 'makeTrapezoid';
95 parser.addRequired('channel',...
96 @(x) any(validatestring(x,validChannels)));
97 parser.addOptional('system',[],@isstruct);
98 parser.addParamValue('duration',0,@(x)(isnumeric(x) && x>0));
99 parser.addParamValue('area',[],@isnumeric);
100 parser.addParamValue('flatTime',[],@isnumeric);
101 parser.addParamValue('flatArea',[],@isnumeric);
102 parser.addParamValue('amplitude',[],@isnumeric);
103 parser.addParamValue('maxGrad',0,@isnumeric);
104 parser.addParamValue('maxSlew',0,@isnumeric);
105 parser.addParamValue('riseTime',0,@isnumeric);
106 parser.addParamValue('fallTime',0,@isnumeric);
107 parser.addParamValue('delay',0,@isnumeric);
108end
109parse(parser,channel,varargin{:});
110opt = parser.Results;
112if isempty(opt.system)
113 system=mr.opts();
114else
115 system=opt.system;
116end
118maxSlew=system.maxSlew;
119%riseTime=system.riseTime;
120maxGrad=system.maxGrad;
121fallTime = [];
122riseTime = [];
124if opt.maxGrad>0
125 maxGrad=opt.maxGrad;
126end
127if opt.maxSlew>0
128 maxSlew=opt.maxSlew;
129end
130if opt.riseTime>0
131 riseTime=opt.riseTime;
132end
133if opt.fallTime>0
134 if isempty(riseTime)
135 error('makeTrapezoid:invalidArguments','Must always supply ''riseTime'' if ''fallTime'' is specified explicitly.');
136 end
137 fallTime=opt.fallTime;
138end
141if (isempty(opt.area)+isempty(opt.flatArea)+isempty(opt.amplitude))~=2
142 error('makeTrapezoid:invalidArguments','Must supply either ''area'', ''flatArea'' or ''amplitude'', and only one of the three may be specified');
143end
144if ~isempty(opt.flatTime) % MZ was: opt.flatTime>0
145 if ~isempty(opt.amplitude)
146 amplitude = opt.amplitude;
147 else
148 if isempty(opt.flatArea)
149 error('makeTrapezoid:invalidArguments','When ''flatTime'' is provided either ''flatArea'' or ''amplitude'' must be provided as well; you may consider providing ''duration'', ''area'' and optionally ramp times instead.');
150 end
151 amplitude = opt.flatArea/opt.flatTime;
152 end
153 if isempty(riseTime)
154 riseTime = abs(amplitude)/maxSlew;
155 riseTime = ceil(riseTime/system.gradRasterTime)*system.gradRasterTime;
156 if riseTime==0
157 riseTime=system.gradRasterTime;
158 end
159 end
160 if isempty(fallTime)
161 fallTime = riseTime;
162 end
163 flatTime = opt.flatTime;
164elseif opt.duration>0
165 if ~isempty(opt.amplitude)
166 amplitude = opt.amplitude;
167 else
168 if isempty(riseTime)
169 dC = 1/abs(2*maxSlew) + 1/abs(2*maxSlew);
170 possible = opt.duration^2 > 4*abs(opt.area)*dC;
171 if ~possible
172 [~, t1, t2, t3]=calcShortestParamsForArea(opt.area,maxSlew,maxGrad,system.gradRasterTime);
173 error('makeTrapezoid:invalidDuration',['Requested area is too large for this gradient. Minimum required duration for this area (accounting for the gradient raster time) is ' num2str((t1+t2+t3)*1e6) 'us']);
174 end
175 amplitude = ( opt.duration - sqrt(opt.duration^2 - 4*abs(opt.area)*dC) )/(2*dC);
176 else
177 if isempty(fallTime)
178 fallTime = riseTime;
179 end
180 amplitude = opt.area/(opt.duration-0.5*riseTime-0.5*fallTime);
181 possible = opt.duration>=(riseTime+fallTime) & abs(amplitude)<maxGrad;
182 assert(possible,['Requested area is too large for this gradient duration. Probably amplitude is violated (' num2str(round(abs(amplitude)/maxGrad*100)) '%)']);
183 end
184 end
185 if isempty(riseTime)
186 riseTime = ceil(abs(amplitude)/maxSlew/system.gradRasterTime)*system.gradRasterTime;
187 if(riseTime==0)
188 riseTime=system.gradRasterTime;
189 end
190 end
191 if isempty(fallTime)
192 fallTime = riseTime;
193 end
194 flatTime = opt.duration-riseTime-fallTime;
195 if isempty(opt.amplitude)
196 % Adjust amplitude (after rounding) to achieve given area
197 amplitude = opt.area/(riseTime/2 + fallTime/2 + flatTime);
198 end
199else
200 if isempty(opt.area)
201 error('makeTrapezoid:invalidArguments','Must supply area or duration');
202 else
203 % call the local function to calculate the shortest timing
204 [amplitude, riseTime, flatTime, fallTime]=calcShortestParamsForArea(opt.area,maxSlew,maxGrad,system.gradRasterTime);
205 end
206end
207if abs(amplitude)>maxGrad
208 if isempty(opt.area)
209 error('makeTrapezoid:invalidAmplitude',['Amplitude violation (' num2str(round(abs(amplitude)/maxGrad*100)) '%%)']);
210 else
211 % this error can only be produced by the failed trapezoid
212 % calculation with the specified area (leading to exceedingly high
213 % amplitude), the triangular blip error should have occured around
214 % line 103
215 [~, t1, t2, t3]=calcShortestParamsForArea(opt.area,maxSlew,maxGrad,system.gradRasterTime);
216 error('makeTrapezoid:invalidDuration',['Requested duration is too short for the area to be realized within system limits. Minimum duration for this trapezoid (accounting for the gradient raster time) is ' num2str((t1+t2+t3)*1e6) ' us']);
217 end
218end
220grad.type = 'trap';
221grad.channel = opt.channel;
222grad.amplitude = amplitude;
223grad.riseTime = riseTime;
224grad.flatTime = flatTime;
225grad.fallTime = fallTime;
226grad.area = amplitude*(flatTime + riseTime/2 + fallTime/2);
227grad.flatArea = amplitude*flatTime;
228grad.delay = opt.delay;
229grad.first = 0;
230grad.last = 0;
232end
234function [amplitude, riseTime, flatTime, fallTime] = calcShortestParamsForArea(area,maxSlew,maxGrad,gradRasterTime)
235 % find the shortest possible duration
236 % first check if the area can be realized as a triangle
237 % if not we calculate a trapezoid
238 riseTime=ceil(sqrt(abs(area)/maxSlew)/gradRasterTime)*gradRasterTime;
239 if riseTime < gradRasterTime % the "area" was probably 0 or almost 0 ...
240 riseTime=gradRasterTime;
241 end
242 amplitude=area/riseTime;
243 tEff=riseTime;
244 if abs(amplitude)>maxGrad
245 tEff=ceil(abs(area)/maxGrad/gradRasterTime)*gradRasterTime;
246 amplitude=area/tEff;
247 riseTime=ceil(abs(amplitude)/maxSlew/gradRasterTime)*gradRasterTime;
248 if(riseTime==0)
249 riseTime=gradRasterTime;
250 end
251 end
252 flatTime=tEff-riseTime;
253 fallTime=riseTime;
254end