function [ ampl, theta ] = r8_b0mp ( x ) %*****************************************************************************80 % %% R8_B0MP evaluates the modulus and phase for the Bessel J0 and Y0 functions. % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 02 October 2011 % % Author: % % Original FORTRAN77 version by Wayne Fullerton. % MATLAB version by John Burkardt. % % Reference: % % Wayne Fullerton, % Portable Special Function Routines, % in Portability of Numerical Software, % edited by Wayne Cowell, % Lecture Notes in Computer Science, Volume 57, % Springer 1977, % ISBN: 978-3-540-08446-4, % LC: QA297.W65. % % Parameters: % % Input, real X, the argument. % % Output, real AMPL, THETA, the modulus and phase. % persistent bm0cs persistent bm02cs persistent bt02cs persistent bth0cs persistent nbm0 persistent nbm02 persistent nbt02 persistent nbth0 persistent xmax if ( isempty ( nbm0 ) ) bm0cs = [ ... +0.9211656246827742712573767730182E-01, ... -0.1050590997271905102480716371755E-02, ... +0.1470159840768759754056392850952E-04, ... -0.5058557606038554223347929327702E-06, ... +0.2787254538632444176630356137881E-07, ... -0.2062363611780914802618841018973E-08, ... +0.1870214313138879675138172596261E-09, ... -0.1969330971135636200241730777825E-10, ... +0.2325973793999275444012508818052E-11, ... -0.3009520344938250272851224734482E-12, ... +0.4194521333850669181471206768646E-13, ... -0.6219449312188445825973267429564E-14, ... +0.9718260411336068469601765885269E-15, ... -0.1588478585701075207366635966937E-15, ... +0.2700072193671308890086217324458E-16, ... -0.4750092365234008992477504786773E-17, ... +0.8615128162604370873191703746560E-18, ... -0.1605608686956144815745602703359E-18, ... +0.3066513987314482975188539801599E-19, ... -0.5987764223193956430696505617066E-20, ... +0.1192971253748248306489069841066E-20, ... -0.2420969142044805489484682581333E-21, ... +0.4996751760510616453371002879999E-22, ... -0.1047493639351158510095040511999E-22, ... +0.2227786843797468101048183466666E-23, ... -0.4801813239398162862370542933333E-24, ... +0.1047962723470959956476996266666E-24, ... -0.2313858165678615325101260800000E-25, ... +0.5164823088462674211635199999999E-26, ... -0.1164691191850065389525401599999E-26, ... +0.2651788486043319282958336000000E-27, ... -0.6092559503825728497691306666666E-28, ... +0.1411804686144259308038826666666E-28, ... -0.3298094961231737245750613333333E-29, ... +0.7763931143074065031714133333333E-30, ... -0.1841031343661458478421333333333E-30, ... +0.4395880138594310737100799999999E-31 ]'; bm02cs = [ ... +0.9500415145228381369330861335560E-01, ... -0.3801864682365670991748081566851E-03, ... +0.2258339301031481192951829927224E-05, ... -0.3895725802372228764730621412605E-07, ... +0.1246886416512081697930990529725E-08, ... -0.6065949022102503779803835058387E-10, ... +0.4008461651421746991015275971045E-11, ... -0.3350998183398094218467298794574E-12, ... +0.3377119716517417367063264341996E-13, ... -0.3964585901635012700569356295823E-14, ... +0.5286111503883857217387939744735E-15, ... -0.7852519083450852313654640243493E-16, ... +0.1280300573386682201011634073449E-16, ... -0.2263996296391429776287099244884E-17, ... +0.4300496929656790388646410290477E-18, ... -0.8705749805132587079747535451455E-19, ... +0.1865862713962095141181442772050E-19, ... -0.4210482486093065457345086972301E-20, ... +0.9956676964228400991581627417842E-21, ... -0.2457357442805313359605921478547E-21, ... +0.6307692160762031568087353707059E-22, ... -0.1678773691440740142693331172388E-22, ... +0.4620259064673904433770878136087E-23, ... -0.1311782266860308732237693402496E-23, ... +0.3834087564116302827747922440276E-24, ... -0.1151459324077741271072613293576E-24, ... +0.3547210007523338523076971345213E-25, ... -0.1119218385815004646264355942176E-25, ... +0.3611879427629837831698404994257E-26, ... -0.1190687765913333150092641762463E-26, ... +0.4005094059403968131802476449536E-27, ... -0.1373169422452212390595193916017E-27, ... +0.4794199088742531585996491526437E-28, ... -0.1702965627624109584006994476452E-28, ... +0.6149512428936330071503575161324E-29, ... -0.2255766896581828349944300237242E-29, ... +0.8399707509294299486061658353200E-30, ... -0.3172997595562602355567423936152E-30, ... +0.1215205298881298554583333026514E-30, ... -0.4715852749754438693013210568045E-31 ]'; bt02cs = [ ... -0.24548295213424597462050467249324, ... +0.12544121039084615780785331778299E-02, ... -0.31253950414871522854973446709571E-04, ... +0.14709778249940831164453426969314E-05, ... -0.99543488937950033643468850351158E-07, ... +0.85493166733203041247578711397751E-08, ... -0.86989759526554334557985512179192E-09, ... +0.10052099533559791084540101082153E-09, ... -0.12828230601708892903483623685544E-10, ... +0.17731700781805131705655750451023E-11, ... -0.26174574569485577488636284180925E-12, ... +0.40828351389972059621966481221103E-13, ... -0.66751668239742720054606749554261E-14, ... +0.11365761393071629448392469549951E-14, ... -0.20051189620647160250559266412117E-15, ... +0.36497978794766269635720591464106E-16, ... -0.68309637564582303169355843788800E-17, ... +0.13107583145670756620057104267946E-17, ... -0.25723363101850607778757130649599E-18, ... +0.51521657441863959925267780949333E-19, ... -0.10513017563758802637940741461333E-19, ... +0.21820381991194813847301084501333E-20, ... -0.46004701210362160577225905493333E-21, ... +0.98407006925466818520953651199999E-22, ... -0.21334038035728375844735986346666E-22, ... +0.46831036423973365296066286933333E-23, ... -0.10400213691985747236513382399999E-23, ... +0.23349105677301510051777740800000E-24, ... -0.52956825323318615788049749333333E-25, ... +0.12126341952959756829196287999999E-25, ... -0.28018897082289428760275626666666E-26, ... +0.65292678987012873342593706666666E-27, ... -0.15337980061873346427835733333333E-27, ... +0.36305884306364536682359466666666E-28, ... -0.86560755713629122479172266666666E-29, ... +0.20779909972536284571238399999999E-29, ... -0.50211170221417221674325333333333E-30, ... +0.12208360279441714184191999999999E-30, ... -0.29860056267039913454250666666666E-31 ]'; bth0cs = [ ... -0.24901780862128936717709793789967, ... +0.48550299609623749241048615535485E-03, ... -0.54511837345017204950656273563505E-05, ... +0.13558673059405964054377445929903E-06, ... -0.55691398902227626227583218414920E-08, ... +0.32609031824994335304004205719468E-09, ... -0.24918807862461341125237903877993E-10, ... +0.23449377420882520554352413564891E-11, ... -0.26096534444310387762177574766136E-12, ... +0.33353140420097395105869955014923E-13, ... -0.47890000440572684646750770557409E-14, ... +0.75956178436192215972642568545248E-15, ... -0.13131556016891440382773397487633E-15, ... +0.24483618345240857495426820738355E-16, ... -0.48805729810618777683256761918331E-17, ... +0.10327285029786316149223756361204E-17, ... -0.23057633815057217157004744527025E-18, ... +0.54044443001892693993017108483765E-19, ... -0.13240695194366572724155032882385E-19, ... +0.33780795621371970203424792124722E-20, ... -0.89457629157111779003026926292299E-21, ... +0.24519906889219317090899908651405E-21, ... -0.69388422876866318680139933157657E-22, ... +0.20228278714890138392946303337791E-22, ... -0.60628500002335483105794195371764E-23, ... +0.18649748964037635381823788396270E-23, ... -0.58783732384849894560245036530867E-24, ... +0.18958591447999563485531179503513E-24, ... -0.62481979372258858959291620728565E-25, ... +0.21017901684551024686638633529074E-25, ... -0.72084300935209253690813933992446E-26, ... +0.25181363892474240867156405976746E-26, ... -0.89518042258785778806143945953643E-27, ... +0.32357237479762298533256235868587E-27, ... -0.11883010519855353657047144113796E-27, ... +0.44306286907358104820579231941731E-28, ... -0.16761009648834829495792010135681E-28, ... +0.64292946921207466972532393966088E-29, ... -0.24992261166978652421207213682763E-29, ... +0.98399794299521955672828260355318E-30, ... -0.39220375242408016397989131626158E-30, ... +0.15818107030056522138590618845692E-30, ... -0.64525506144890715944344098365426E-31, ... +0.26611111369199356137177018346367E-31 ]'; eta = 0.1 * r8_mach ( 3 ); nbm0 = r8_inits ( bm0cs, 37, eta ); nbt02 = r8_inits ( bt02cs, 39, eta ); nbm02 = r8_inits ( bm02cs, 40, eta ); nbth0 = r8_inits ( bth0cs, 44, eta ); xmax = 1.0 / r8_mach ( 4 ); end if ( x < 4.0 ) fprintf ( 1, '\n' ); fprintf ( 1, 'R8_B0MP - Fatal error!\n' ); fprintf ( 1, ' X < 4.\n' ); error ( 'R8_B0MP - Fatal error!' ) elseif ( x <= 8.0 ) z = ( 128.0 / x / x - 5.0 ) / 3.0; ampl = ( 0.75 + r8_csevl ( z, bm0cs, nbm0 ) ) / sqrt ( x ); theta = x - 0.25 * pi + r8_csevl ( z, bt02cs, nbt02 ) / x; else z = 128.0 / x / x - 1.0; ampl = ( 0.75 + r8_csevl ( z, bm02cs, nbm02) ) / sqrt ( x ); theta = x - 0.25 * pi + r8_csevl ( z, bth0cs, nbth0 ) / x; end return end