function value = r8_besi1e ( x ) %*****************************************************************************80 % %% R8_BESI1E evaluates the exponentially scaled Bessel function I1(X). % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 28 September 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 VALUE, the exponentially scaled Bessel function I1(X). % persistent ai12cs persistent ai1cs persistent bi1cs persistent ntai1 persistent ntai12 persistent nti1 persistent xmin persistent xsml if ( isempty ( nti1 ) ) ai12cs = [ ... +0.2857623501828012047449845948469E-01, ... -0.9761097491361468407765164457302E-02, ... -0.1105889387626237162912569212775E-03, ... -0.3882564808877690393456544776274E-05, ... -0.2512236237870208925294520022121E-06, ... -0.2631468846889519506837052365232E-07, ... -0.3835380385964237022045006787968E-08, ... -0.5589743462196583806868112522229E-09, ... -0.1897495812350541234498925033238E-10, ... +0.3252603583015488238555080679949E-10, ... +0.1412580743661378133163366332846E-10, ... +0.2035628544147089507224526136840E-11, ... -0.7198551776245908512092589890446E-12, ... -0.4083551111092197318228499639691E-12, ... -0.2101541842772664313019845727462E-13, ... +0.4272440016711951354297788336997E-13, ... +0.1042027698412880276417414499948E-13, ... -0.3814403072437007804767072535396E-14, ... -0.1880354775510782448512734533963E-14, ... +0.3308202310920928282731903352405E-15, ... +0.2962628997645950139068546542052E-15, ... -0.3209525921993423958778373532887E-16, ... -0.4650305368489358325571282818979E-16, ... +0.4414348323071707949946113759641E-17, ... +0.7517296310842104805425458080295E-17, ... -0.9314178867326883375684847845157E-18, ... -0.1242193275194890956116784488697E-17, ... +0.2414276719454848469005153902176E-18, ... +0.2026944384053285178971922860692E-18, ... -0.6394267188269097787043919886811E-19, ... -0.3049812452373095896084884503571E-19, ... +0.1612841851651480225134622307691E-19, ... +0.3560913964309925054510270904620E-20, ... -0.3752017947936439079666828003246E-20, ... -0.5787037427074799345951982310741E-22, ... +0.7759997511648161961982369632092E-21, ... -0.1452790897202233394064459874085E-21, ... -0.1318225286739036702121922753374E-21, ... +0.6116654862903070701879991331717E-22, ... +0.1376279762427126427730243383634E-22, ... -0.1690837689959347884919839382306E-22, ... +0.1430596088595433153987201085385E-23, ... +0.3409557828090594020405367729902E-23, ... -0.1309457666270760227845738726424E-23, ... -0.3940706411240257436093521417557E-24, ... +0.4277137426980876580806166797352E-24, ... -0.4424634830982606881900283123029E-25, ... -0.8734113196230714972115309788747E-25, ... +0.4045401335683533392143404142428E-25, ... +0.7067100658094689465651607717806E-26, ... -0.1249463344565105223002864518605E-25, ... +0.2867392244403437032979483391426E-26, ... +0.2044292892504292670281779574210E-26, ... -0.1518636633820462568371346802911E-26, ... +0.8110181098187575886132279107037E-28, ... +0.3580379354773586091127173703270E-27, ... -0.1692929018927902509593057175448E-27, ... -0.2222902499702427639067758527774E-28, ... +0.5424535127145969655048600401128E-28, ... -0.1787068401578018688764912993304E-28, ... -0.6565479068722814938823929437880E-29, ... +0.7807013165061145280922067706839E-29, ... -0.1816595260668979717379333152221E-29, ... -0.1287704952660084820376875598959E-29, ... +0.1114548172988164547413709273694E-29, ... -0.1808343145039336939159368876687E-30, ... -0.2231677718203771952232448228939E-30, ... +0.1619029596080341510617909803614E-30, ... -0.1834079908804941413901308439210E-31 ]'; ai1cs = [ ... -0.2846744181881478674100372468307E-01, ... -0.1922953231443220651044448774979E-01, ... -0.6115185857943788982256249917785E-03, ... -0.2069971253350227708882823777979E-04, ... +0.8585619145810725565536944673138E-05, ... +0.1049498246711590862517453997860E-05, ... -0.2918338918447902202093432326697E-06, ... -0.1559378146631739000160680969077E-07, ... +0.1318012367144944705525302873909E-07, ... -0.1448423418183078317639134467815E-08, ... -0.2908512243993142094825040993010E-09, ... +0.1266388917875382387311159690403E-09, ... -0.1664947772919220670624178398580E-10, ... -0.1666653644609432976095937154999E-11, ... +0.1242602414290768265232168472017E-11, ... -0.2731549379672432397251461428633E-12, ... +0.2023947881645803780700262688981E-13, ... +0.7307950018116883636198698126123E-14, ... -0.3332905634404674943813778617133E-14, ... +0.7175346558512953743542254665670E-15, ... -0.6982530324796256355850629223656E-16, ... -0.1299944201562760760060446080587E-16, ... +0.8120942864242798892054678342860E-17, ... -0.2194016207410736898156266643783E-17, ... +0.3630516170029654848279860932334E-18, ... -0.1695139772439104166306866790399E-19, ... -0.1288184829897907807116882538222E-19, ... +0.5694428604967052780109991073109E-20, ... -0.1459597009090480056545509900287E-20, ... +0.2514546010675717314084691334485E-21, ... -0.1844758883139124818160400029013E-22, ... -0.6339760596227948641928609791999E-23, ... +0.3461441102031011111108146626560E-23, ... -0.1017062335371393547596541023573E-23, ... +0.2149877147090431445962500778666E-24, ... -0.3045252425238676401746206173866E-25, ... +0.5238082144721285982177634986666E-27, ... +0.1443583107089382446416789503999E-26, ... -0.6121302074890042733200670719999E-27, ... +0.1700011117467818418349189802666E-27, ... -0.3596589107984244158535215786666E-28, ... +0.5448178578948418576650513066666E-29, ... -0.2731831789689084989162564266666E-30, ... -0.1858905021708600715771903999999E-30, ... +0.9212682974513933441127765333333E-31, ... -0.2813835155653561106370833066666E-31 ]'; bi1cs = [ ... -0.19717132610998597316138503218149E-02, ... +0.40734887667546480608155393652014, ... +0.34838994299959455866245037783787E-01, ... +0.15453945563001236038598401058489E-02, ... +0.41888521098377784129458832004120E-04, ... +0.76490267648362114741959703966069E-06, ... +0.10042493924741178689179808037238E-07, ... +0.99322077919238106481371298054863E-10, ... +0.76638017918447637275200171681349E-12, ... +0.47414189238167394980388091948160E-14, ... +0.24041144040745181799863172032000E-16, ... +0.10171505007093713649121100799999E-18, ... +0.36450935657866949458491733333333E-21, ... +0.11205749502562039344810666666666E-23, ... +0.29875441934468088832000000000000E-26, ... +0.69732310939194709333333333333333E-29, ... +0.14367948220620800000000000000000E-31 ]'; eta = 0.1 * r8_mach ( 3 ); nti1 = r8_inits ( bi1cs, 17, eta ); ntai1 = r8_inits ( ai1cs, 46, eta ); ntai12 = r8_inits ( ai12cs, 69, eta ); xmin = 2.0 * r8_mach ( 1 ); xsml = sqrt ( 8.0 * r8_mach ( 3 ) ); end y = abs ( x ); if ( y <= xmin ) value = 0.0; elseif ( y <= xsml ) value = 0.5 * x * exp ( - y ); elseif ( y <= 3.0 ) value = x * ( 0.875 + r8_csevl ( y * y / 4.5 - 1.0, bi1cs, nti1 ) ) ... * exp ( - y ); elseif ( y <= 8.0 ) value = ( 0.375 + r8_csevl ( ( 48.0 / y - 11.0) / 5.0, ... ai1cs, ntai1 ) ) / sqrt ( y ); if ( x < 0.0 ) value = - value; end else value = ( 0.375 + r8_csevl ( 16.0 / y - 1.0, ai12cs, ntai12 ) ) ... / sqrt ( y ); if ( x < 0.0 ) value = - value; end end return end