function [ f, g ] = r8_sifg ( x ) %*****************************************************************************80 % %% R8_SIFG is a utility routine. % % 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 F, G. % persistent f1cs persistent f2cs persistent g1cs persistent g2cs persistent g3cs persistent nf1 persistent nf2 persistent ng1 persistent ng2 persistent ng3 persistent xbig persistent xbnd persistent xbndg persistent xmaxf persistent xmaxg if ( isempty ( nf1 ) ) f1cs = [ ... -0.1191081969051363610348201965828918, ... -0.0247823144996236247590074150823133, ... 0.0011910281453357821268120363054457, ... -0.0000927027714388561748308600360706, ... 0.0000093373141568270996868204582766, ... -0.0000011058287820557143938979426306, ... 0.0000001464772071460162169336550799, ... -0.0000000210694496287689532601227548, ... 0.0000000032293492366848236382857374, ... -0.0000000005206529617529375828014986, ... 0.0000000000874878884570278750268316, ... -0.0000000000152176187056123668294574, ... 0.0000000000027257192405419573900583, ... -0.0000000000005007053075968556290255, ... 0.0000000000000940240902726068511779, ... -0.0000000000000180014444791803678336, ... 0.0000000000000035062621432741785826, ... -0.0000000000000006935282926769149709, ... 0.0000000000000001390925136454216568, ... -0.0000000000000000282486885074170585, ... 0.0000000000000000058031305693579081, ... -0.0000000000000000012046901573375820, ... 0.0000000000000000002525052443655940, ... -0.0000000000000000000533980268805594, ... 0.0000000000000000000113855786274122, ... -0.0000000000000000000024462861505259, ... 0.0000000000000000000005293659320439, ... -0.0000000000000000000001153184940277, ... 0.0000000000000000000000252786568318, ... -0.0000000000000000000000055738645378, ... 0.0000000000000000000000012358245621, ... -0.0000000000000000000000002754350842, ... 0.0000000000000000000000000616906808, ... -0.0000000000000000000000000138817443, ... 0.0000000000000000000000000031375329, ... -0.0000000000000000000000000007121249, ... 0.0000000000000000000000000001622778, ... -0.0000000000000000000000000000371206, ... 0.0000000000000000000000000000085221, ... -0.0000000000000000000000000000019633, ... 0.0000000000000000000000000000004538, ... -0.0000000000000000000000000000001052, ... 0.0000000000000000000000000000000245 ]'; f2cs = [ ... -0.03484092538970132330836049733745577, ... -0.01668422056779596873246786312278676, ... 0.00067529012412377385045207859239727, ... -0.00005350666225447013628785577557429, ... 0.00000626934217790075267050759431626, ... -0.00000095266388019916680677790414293, ... 0.00000017456292242509880425504427666, ... -0.00000003687954030653093307097646628, ... 0.00000000872026777051395264075816938, ... -0.00000000226019703919738748530423167, ... 0.00000000063246249765250612520444877, ... -0.00000000018889118884717869240911480, ... 0.00000000005967746729997813372620472, ... -0.00000000001980443117372239011196007, ... 0.00000000000686413954772103383713264, ... -0.00000000000247310193070199106074890, ... 0.00000000000092263594549941404196042, ... -0.00000000000035523634999261784497297, ... 0.00000000000014076049625351591461820, ... -0.00000000000005726228499747652794311, ... 0.00000000000002386537545413171810106, ... -0.00000000000001017141890764597142232, ... 0.00000000000000442594531078364424968, ... -0.00000000000000196344933049189761979, ... 0.00000000000000088688748314810461024, ... -0.00000000000000040743345027311546948, ... 0.00000000000000019016837215675339859, ... -0.00000000000000009009707297478042442, ... 0.00000000000000004329211274095668667, ... -0.00000000000000002108144465322479526, ... 0.00000000000000001039637907026452274, ... -0.00000000000000000518891007948931936, ... 0.00000000000000000261955324869899371, ... -0.00000000000000000133690399951301570, ... 0.00000000000000000068941057702931664, ... -0.00000000000000000035905362610437250, ... 0.00000000000000000018878077255791706, ... -0.00000000000000000010016125265594380, ... 0.00000000000000000005360725691578228, ... -0.00000000000000000002893198974944827, ... 0.00000000000000000001574065100202625, ... -0.00000000000000000000863027106431206, ... 0.00000000000000000000476715602862288, ... -0.00000000000000000000265222739998504, ... 0.00000000000000000000148582865063866, ... -0.00000000000000000000083797235923135, ... 0.00000000000000000000047565916422711, ... -0.00000000000000000000027169073353112, ... 0.00000000000000000000015612738881686, ... -0.00000000000000000000009024555078347, ... 0.00000000000000000000005246097049119, ... -0.00000000000000000000003066450818697, ... 0.00000000000000000000001801996250957, ... -0.00000000000000000000001064443050752, ... 0.00000000000000000000000631942158881, ... -0.00000000000000000000000377013812246, ... 0.00000000000000000000000225997542918, ... -0.00000000000000000000000136100844814, ... 0.00000000000000000000000082333232003, ... -0.00000000000000000000000050025986091, ... 0.00000000000000000000000030526245684, ... -0.00000000000000000000000018705164021, ... 0.00000000000000000000000011508404393, ... -0.00000000000000000000000007108714611, ... 0.00000000000000000000000004408065533, ... -0.00000000000000000000000002743760867, ... 0.00000000000000000000000001714144851, ... -0.00000000000000000000000001074768860, ... 0.00000000000000000000000000676259777, ... -0.00000000000000000000000000426981348, ... 0.00000000000000000000000000270500637, ... -0.00000000000000000000000000171933331, ... 0.00000000000000000000000000109636138, ... -0.00000000000000000000000000070132573, ... 0.00000000000000000000000000045001784, ... -0.00000000000000000000000000028963835, ... 0.00000000000000000000000000018697009, ... -0.00000000000000000000000000012104646, ... 0.00000000000000000000000000007859065, ... -0.00000000000000000000000000005116867, ... 0.00000000000000000000000000003340627, ... -0.00000000000000000000000000002186851, ... 0.00000000000000000000000000001435340, ... -0.00000000000000000000000000000944523, ... 0.00000000000000000000000000000623117, ... -0.00000000000000000000000000000412101, ... 0.00000000000000000000000000000273208, ... -0.00000000000000000000000000000181558, ... 0.00000000000000000000000000000120934, ... -0.00000000000000000000000000000080737, ... 0.00000000000000000000000000000054022, ... -0.00000000000000000000000000000036227, ... 0.00000000000000000000000000000024348, ... -0.00000000000000000000000000000016401, ... 0.00000000000000000000000000000011074, ... -0.00000000000000000000000000000007497, ... 0.00000000000000000000000000000005091, ... -0.00000000000000000000000000000003470, ... 0.00000000000000000000000000000002377 ]'; g1cs = [ ... -0.3040578798253495954499726682091083, ... -0.0566890984597120587731339156118269, ... 0.0039046158173275643919984071554082, ... -0.0003746075959202260618619339867489, ... 0.0000435431556559843679552220840065, ... -0.0000057417294453025046561970723475, ... 0.0000008282552104502629741937616492, ... -0.0000001278245892594642727883913223, ... 0.0000000207978352948687884439257529, ... -0.0000000035313205921990798042032682, ... 0.0000000006210824236308951068631449, ... -0.0000000001125215474446292649336987, ... 0.0000000000209088917684421605267019, ... -0.0000000000039715831737681727689158, ... 0.0000000000007690431314272089939005, ... -0.0000000000001514696742731613519826, ... 0.0000000000000302892146552359684119, ... -0.0000000000000061399703834708825400, ... 0.0000000000000012600605829510933553, ... -0.0000000000000002615029250939483683, ... 0.0000000000000000548278844891796821, ... -0.0000000000000000116038182129526571, ... 0.0000000000000000024771654107129795, ... -0.0000000000000000005330672753223389, ... 0.0000000000000000001155666075598465, ... -0.0000000000000000000252280547744957, ... 0.0000000000000000000055429038550786, ... -0.0000000000000000000012252208421297, ... 0.0000000000000000000002723664318684, ... -0.0000000000000000000000608707831422, ... 0.0000000000000000000000136724874476, ... -0.0000000000000000000000030856626806, ... 0.0000000000000000000000006995212319, ... -0.0000000000000000000000001592587569, ... 0.0000000000000000000000000364051056, ... -0.0000000000000000000000000083539465, ... 0.0000000000000000000000000019240303, ... -0.0000000000000000000000000004446816, ... 0.0000000000000000000000000001031182, ... -0.0000000000000000000000000000239887, ... 0.0000000000000000000000000000055976, ... -0.0000000000000000000000000000013100, ... 0.0000000000000000000000000000003074, ... -0.0000000000000000000000000000000723 ]'; g2cs = [ ... -0.1211802894731646263541834046858267, ... -0.0316761386394950286701407923505610, ... 0.0013383199778862680163819429492182, ... -0.0000895511011392252425531905069518, ... 0.0000079155562961718213115249467924, ... -0.0000008438793322241520181418982080, ... 0.0000001029980425677530146647227274, ... -0.0000000139295750605183835795834444, ... 0.0000000020422703959875980400677594, ... -0.0000000003196534694206427035434752, ... 0.0000000000528147832657267698615312, ... -0.0000000000091339554672671033735289, ... 0.0000000000016426251238967760444819, ... -0.0000000000003055897039322660002410, ... 0.0000000000000585655825785779717892, ... -0.0000000000000115229197730940120563, ... 0.0000000000000023209469119988537310, ... -0.0000000000000004774355834177535025, ... 0.0000000000000001000996765800180573, ... -0.0000000000000000213533778082256704, ... 0.0000000000000000046277190777367671, ... -0.0000000000000000010175807410227657, ... 0.0000000000000000002267657399884672, ... -0.0000000000000000000511630776076426, ... 0.0000000000000000000116767014913108, ... -0.0000000000000000000026935427672470, ... 0.0000000000000000000006275665841146, ... -0.0000000000000000000001475880557531, ... 0.0000000000000000000000350145314739, ... -0.0000000000000000000000083757732152, ... 0.0000000000000000000000020191815152, ... -0.0000000000000000000000004903567705, ... 0.0000000000000000000000001199123348, ... -0.0000000000000000000000000295170610, ... 0.0000000000000000000000000073113112, ... -0.0000000000000000000000000018217843, ... 0.0000000000000000000000000004565148, ... -0.0000000000000000000000000001150151, ... 0.0000000000000000000000000000291267, ... -0.0000000000000000000000000000074125, ... 0.0000000000000000000000000000018953, ... -0.0000000000000000000000000000004868, ... 0.0000000000000000000000000000001256, ... -0.0000000000000000000000000000000325 ]'; g3cs = [ ... -0.0280574367809472928402815264335299, ... -0.0137271597162236975409100508089556, ... 0.0002894032638760296027448941273751, ... -0.0000114129239391197145908743622517, ... 0.0000006813965590726242997720207302, ... -0.0000000547952289604652363669058052, ... 0.0000000055207429918212529109406521, ... -0.0000000006641464199322920022491428, ... 0.0000000000922373663487041108564960, ... -0.0000000000144299088886682862611718, ... 0.0000000000024963904892030710248705, ... -0.0000000000004708240675875244722971, ... 0.0000000000000957217659216759988140, ... -0.0000000000000207889966095809030537, ... 0.0000000000000047875099970877431627, ... -0.0000000000000011619070583377173759, ... 0.0000000000000002956508969267836974, ... -0.0000000000000000785294988256492025, ... 0.0000000000000000216922264368256612, ... -0.0000000000000000062113515831676342, ... 0.0000000000000000018384568838450977, ... -0.0000000000000000005610887482137276, ... 0.0000000000000000001761862805280062, ... -0.0000000000000000000568111050541451, ... 0.0000000000000000000187786279582313, ... -0.0000000000000000000063531694151124, ... 0.0000000000000000000021968802368238, ... -0.0000000000000000000007754666550395, ... 0.0000000000000000000002791018356581, ... -0.0000000000000000000001023178525247, ... 0.0000000000000000000000381693403919, ... -0.0000000000000000000000144767895606, ... 0.0000000000000000000000055779512634, ... -0.0000000000000000000000021817239071, ... 0.0000000000000000000000008656646309, ... -0.0000000000000000000000003482157895, ... 0.0000000000000000000000001419188130, ... -0.0000000000000000000000000585714314, ... 0.0000000000000000000000000244660482, ... -0.0000000000000000000000000103387099, ... 0.0000000000000000000000000044177299, ... -0.0000000000000000000000000019080079, ... 0.0000000000000000000000000008326038, ... -0.0000000000000000000000000003669553, ... 0.0000000000000000000000000001632875, ... -0.0000000000000000000000000000733357, ... 0.0000000000000000000000000000332327, ... -0.0000000000000000000000000000151906, ... 0.0000000000000000000000000000070020, ... -0.0000000000000000000000000000032539, ... 0.0000000000000000000000000000015240, ... -0.0000000000000000000000000000007193, ... 0.0000000000000000000000000000003420, ... -0.0000000000000000000000000000001638, ... 0.0000000000000000000000000000000790, ... -0.0000000000000000000000000000000383 ]'; tol = 0.1 * r8_mach ( 3 ); nf1 = r8_inits ( f1cs, 43, tol ); nf2 = r8_inits ( f2cs, 99, tol ); ng1 = r8_inits ( g1cs, 44, tol ); ng2 = r8_inits ( g2cs, 44, tol ); ng3 = r8_inits ( g3cs, 56, tol ); xbig = sqrt ( 1.0 / r8_mach ( 3 ) ); xmaxf = exp ( min ( - log ( r8_mach ( 1 ) ), ... log ( r8_mach ( 2 ) ) ) - 0.01 ); xmaxg = 1.0 / sqrt ( r8_mach ( 1 ) ); xbnd = sqrt ( 50.0 ); xbndg = sqrt ( 200.0 ); end if ( x < 4.0 ) fprintf ( 1, '\n' ); fprintf ( 1, 'R8_SIFG - Fatal error!\n' ); fprintf ( 1, ' Approximation invalid for X < 4.\n' ); error ( 'R8_SIFG - Fatal error!' ) elseif ( x <= xbnd ) f = ( 1.0 + r8_csevl ( ( 1.0 / x / x - 0.04125 ) ... / 0.02125, f1cs, nf1 ) ) / x; g = ( 1.0 + r8_csevl ( ( 1.0 / x / x - 0.04125 ) ... / 0.02125, g1cs, ng1 ) ) / x / x; elseif ( x <= xbig ) f = ( 1.0 + r8_csevl ( 100. / x / x - 1.0, f2cs, nf2 ) ) / x; if ( x <= xbndg ) g = ( 1.0 + r8_csevl ( ( 10000.0 / x / x - 125.0 ) ... / 75.0, g2cs, ng2 ) ) / x / x; else g = ( 1.0 + r8_csevl ( 400.0 / x / x - 1.0, g3cs, ng3 ) ) / x / x; end else if ( x < xmaxf ) f = 1.0 / x; else f = 0.0; end if ( x < xmaxg ) g = 1.0 / x / x; else g = 0.0; end end return end