tp@0: function outval = EDB1betaoverml(zvec,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec)
tp@0: % EDB1betaoverml - Integrand function which is called for num. int. of ED IR.
tp@0: % EDB1betaoverml is copied directly
tp@0: % into the numerical integration function EDB1quadstep for faster operation
tp@0: % but EDB1betaoverml is also called directly by EDbwedge1st_int.
tp@0: %
tp@0: % Uses no special functions.
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 20040323 (svensson@iet.ntnu.no)
tp@0: %
tp@0: % outval = EDB1betaoverml(zvec,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec);
tp@0: 
tp@0: ml = sqrt( (zvec-zs).^2 + rs.^2 ).*sqrt( (zvec-zr).^2 + rr.^2 );
tp@0: 	
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 = real((real(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 + (zvec-zs).*(zvec-zr))./(rs*rr);
tp@0: zvec = real((real(sqrt(zvec.^2-1) + zvec)).^ny);
tp@0: zvec = 0.5*( zvec + 1./zvec);
tp@0: 
tp@0: outval =          sinnyfivec(1)./(zvec - cosnyfivec(1));
tp@0: outval = outval + sinnyfivec(2)./(zvec - cosnyfivec(2));
tp@0: outval = outval + sinnyfivec(3)./(zvec - cosnyfivec(3));
tp@0: outval = outval + sinnyfivec(4)./(zvec - cosnyfivec(4));         
tp@0: 
tp@0: outval = outval./ml;