tp@0: function Q = EDB1quadstep(a,b,fa,fc,fb,tol,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec)
tp@0: % EDB1quadstep - Recursive numerical integration routine for EDB1wedge1st_int.
tp@0: % Modified from Matlabs QUADSTEP function. The EDB1quadstep version
tp@0: % solves all time steps at the same time.
tp@0: %
tp@0: % Uses no special function 
tp@0: %
tp@0: % ----------------------------------------------------------------------------------------------
tp@0: %   This file is part of the Edge Diffraction Toolbox by Peter Svensson.                       
tp@0: %                                                                                              
tp@0: %   The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify       
tp@0: %   it under the terms of the GNU General Public License as published by the Free Software     
tp@0: %   Foundation, either version 3 of the License, or (at your option) any later version.        
tp@0: %                                                                                              
tp@0: %   The Edge Diffraction Toolbox is distributed in the hope that it will be useful,       
tp@0: %   but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS  
tp@0: %   FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.             
tp@0: %                                                                                              
tp@0: %   You should have received a copy of the GNU General Public License along with the           
tp@0: %   Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.                 
tp@0: % ----------------------------------------------------------------------------------------------
tp@0: % Peter Svensson (svensson@iet.ntnu.no) 20040321
tp@0: %
tp@0: % EDB1quadstep(a,b,fa,fc,fb,tol,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec)
tp@0: 
tp@0: % Evaluate integrand twice in interior of subinterval [a,b].
tp@0: c = (a + b)/2;
tp@0: 
tp@0: nslots = length(a);
tp@0: 
tp@0: % We use the matrix x for both the x-interval values and the corresponding
tp@0: % y-values
tp@0: 
tp@0: x = [(a + c)/2, (c + b)/2];
tp@0: x = reshape(x,2*nslots,1);
tp@0: 
tp@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tp@0: %% Instead of using a function call, we plug in the entire function here!
tp@0: %% x = betaoverml_all(x,rs,rr,zr,ny,sinnyfivec,cosnyfivec);
tp@0: 
tp@0: ml = sqrt( (x-zs).^2 + rs^2 ).*sqrt( (x-zr).^2 + rr^2 );
tp@0: %------------------------------------------------
tp@0: % We would have liked to use the following code:
tp@0: %       y = (ml + (zvec-zs).*(zvec-zr))/(rs*rr);
tp@0: %       expnyeta = (sqrt(y.^2-1) + y).^ny;
tp@0: %       coshnyeta = 0.5*( expnyeta + 1./expnyeta);
tp@0: % but to save three vectors, zvec will be re-used for
tp@0: % y and expnyeta and coshnyeta
tp@0: 
tp@0: zvec = (ml + (x-zs).*(x-zr))/(rs*rr);
tp@0: zvec = (real(sqrt(zvec.^2-1)) + zvec).^ny;
tp@0: zvec = 0.5*( zvec + 1./zvec);
tp@0: x =     sinnyfivec(1)./(zvec - cosnyfivec(1));
tp@0: x = x + sinnyfivec(2)./(zvec - cosnyfivec(2));
tp@0: x = x + sinnyfivec(3)./(zvec - cosnyfivec(3));
tp@0: x = x + sinnyfivec(4)./(zvec - cosnyfivec(4));         
tp@0: 
tp@0: x = x./ml;
tp@0: 
tp@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tp@0: 
tp@0: x = reshape(x,nslots,2);
tp@0: 
tp@0: fd = x(:,1);
tp@0: fe = x(:,2);
tp@0: ftemp = fa + fb + 2*fc;
tp@0: 
tp@0: %-----------------------------------------
tp@0: % We use the output vector Q for another purpose
tp@0: 
tp@0: tempvec = (b-a)/6;
tp@0: 
tp@0: % Three point Simpson's rule.
tp@0: Q1 = tempvec.*(ftemp + 2*fc);
tp@0: 
tp@0: % Five point double Simpson's rule.
tp@0: %Q2 = (h/12).*(fa + 4*fd + 2*fc + 4*fe + fb);
tp@0: 
tp@0: % We would have liked to use the following code:
tp@0: % % % Q2 = tempvec/2.*(ftemp + 4*fd + 4*fe);
tp@0: % % % % One step of Romberg extrapolation.
tp@0: % % % Q = Q2 + (Q2 - Q1)/15;
tp@0: % ...... but we save one vector by re-using ftemp
tp@0: % for Q2
tp@0: 
tp@0: ftemp = tempvec/2.*(ftemp + 4*fd + 4*fe);
tp@0: % ftemp = tempvec/2.*(ftemp + 4*x(:,1) + 4*x(:,2));
tp@0: % One step of Romberg extrapolation.
tp@0: Q = ftemp + (ftemp - Q1)/15;
tp@0: 
tp@0: %-----------------------------------------
tp@0: % Check if all or some of the sample intervals 
tp@0: % have reached the tolerance.
tp@0: 
tp@0: tempvec = find(abs(ftemp-Q)>tol);
tp@0: 
tp@0: % Check accuracy of integral over this subinterval.
tp@0: 
tp@0: if isempty(tempvec)
tp@0:    return
tp@0: 
tp@0: % Subdivide into two subintervals.
tp@0: else
tp@0:    Qac = EDB1quadstep(a(tempvec),c(tempvec),fa(tempvec),fd(tempvec),fc(tempvec),tol,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec);
tp@0:    Qcb = EDB1quadstep(c(tempvec),b(tempvec),fc(tempvec),fe(tempvec),fb(tempvec),tol,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec);
tp@0:    Q(tempvec) = Qac + Qcb;
tp@0: end