#! /usr/bin/env python # def r8_e1 ( x ): #*****************************************************************************80 # ## R8_E1 evaluates the exponential integral E1 for an R8 argument. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 25 April 2016 # # Author: # # Original FORTRAN77 version by Wayne Fullerton. # Python 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 exponential integral E1 evaluated at X. # import numpy as np from r8_csevl import r8_csevl from r8_inits import r8_inits from machine import r8_mach ae10cs = np.array ( [ \ +0.3284394579616699087873844201881E-01, \ -0.1669920452031362851476184343387E-01, \ +0.2845284724361346807424899853252E-03, \ -0.7563944358516206489487866938533E-05, \ +0.2798971289450859157504843180879E-06, \ -0.1357901828534531069525563926255E-07, \ +0.8343596202040469255856102904906E-09, \ -0.6370971727640248438275242988532E-10, \ +0.6007247608811861235760831561584E-11, \ -0.7022876174679773590750626150088E-12, \ +0.1018302673703687693096652346883E-12, \ -0.1761812903430880040406309966422E-13, \ +0.3250828614235360694244030353877E-14, \ -0.5071770025505818678824872259044E-15, \ +0.1665177387043294298172486084156E-16, \ +0.3166753890797514400677003536555E-16, \ -0.1588403763664141515133118343538E-16, \ +0.4175513256138018833003034618484E-17, \ -0.2892347749707141906710714478852E-18, \ -0.2800625903396608103506340589669E-18, \ +0.1322938639539270903707580023781E-18, \ -0.1804447444177301627283887833557E-19, \ -0.7905384086522616076291644817604E-20, \ +0.4435711366369570103946235838027E-20, \ -0.4264103994978120868865309206555E-21, \ -0.3920101766937117541553713162048E-21, \ +0.1527378051343994266343752326971E-21, \ +0.1024849527049372339310308783117E-22, \ -0.2134907874771433576262711405882E-22, \ +0.3239139475160028267061694700366E-23, \ +0.2142183762299889954762643168296E-23, \ -0.8234609419601018414700348082312E-24, \ -0.1524652829645809479613694401140E-24, \ +0.1378208282460639134668480364325E-24, \ +0.2131311202833947879523224999253E-26, \ -0.2012649651526484121817466763127E-25, \ +0.1995535662263358016106311782673E-26, \ +0.2798995808984003464948686520319E-26, \ -0.5534511845389626637640819277823E-27, \ -0.3884995396159968861682544026146E-27, \ +0.1121304434507359382850680354679E-27, \ +0.5566568152423740948256563833514E-28, \ -0.2045482929810499700448533938176E-28, \ -0.8453813992712336233411457493674E-29, \ +0.3565758433431291562816111116287E-29, \ +0.1383653872125634705539949098871E-29, \ -0.6062167864451372436584533764778E-30, \ -0.2447198043989313267437655119189E-30, \ +0.1006850640933998348011548180480E-30, \ +0.4623685555014869015664341461674E-31 ] ) ae11cs = np.array ( [ \ +0.20263150647078889499401236517381, \ -0.73655140991203130439536898728034E-01, \ +0.63909349118361915862753283840020E-02, \ -0.60797252705247911780653153363999E-03, \ -0.73706498620176629330681411493484E-04, \ +0.48732857449450183453464992488076E-04, \ -0.23837064840448290766588489460235E-05, \ -0.30518612628561521027027332246121E-05, \ +0.17050331572564559009688032992907E-06, \ +0.23834204527487747258601598136403E-06, \ +0.10781772556163166562596872364020E-07, \ -0.17955692847399102653642691446599E-07, \ -0.41284072341950457727912394640436E-08, \ +0.68622148588631968618346844526664E-09, \ +0.53130183120506356147602009675961E-09, \ +0.78796880261490694831305022893515E-10, \ -0.26261762329356522290341675271232E-10, \ -0.15483687636308261963125756294100E-10, \ -0.25818962377261390492802405122591E-11, \ +0.59542879191591072658903529959352E-12, \ +0.46451400387681525833784919321405E-12, \ +0.11557855023255861496288006203731E-12, \ -0.10475236870835799012317547189670E-14, \ -0.11896653502709004368104489260929E-13, \ -0.47749077490261778752643019349950E-14, \ -0.81077649615772777976249734754135E-15, \ +0.13435569250031554199376987998178E-15, \ +0.14134530022913106260248873881287E-15, \ +0.49451592573953173115520663232883E-16, \ +0.79884048480080665648858587399367E-17, \ -0.14008632188089809829248711935393E-17, \ -0.14814246958417372107722804001680E-17, \ -0.55826173646025601904010693937113E-18, \ -0.11442074542191647264783072544598E-18, \ +0.25371823879566853500524018479923E-20, \ +0.13205328154805359813278863389097E-19, \ +0.62930261081586809166287426789485E-20, \ +0.17688270424882713734999261332548E-20, \ +0.23266187985146045209674296887432E-21, \ -0.67803060811125233043773831844113E-22, \ -0.59440876959676373802874150531891E-22, \ -0.23618214531184415968532592503466E-22, \ -0.60214499724601478214168478744576E-23, \ -0.65517906474348299071370444144639E-24, \ +0.29388755297497724587042038699349E-24, \ +0.22601606200642115173215728758510E-24, \ +0.89534369245958628745091206873087E-25, \ +0.24015923471098457555772067457706E-25, \ +0.34118376888907172955666423043413E-26, \ -0.71617071694630342052355013345279E-27, \ -0.75620390659281725157928651980799E-27, \ -0.33774612157467324637952920780800E-27, \ -0.10479325703300941711526430332245E-27, \ -0.21654550252170342240854880201386E-28, \ -0.75297125745288269994689298432000E-30, \ +0.19103179392798935768638084000426E-29, \ +0.11492104966530338547790728833706E-29, \ +0.43896970582661751514410359193600E-30, \ +0.12320883239205686471647157725866E-30, \ +0.22220174457553175317538581162666E-31 ] ) ae12cs = np.array ( [ \ +0.63629589796747038767129887806803, \ -0.13081168675067634385812671121135, \ -0.84367410213053930014487662129752E-02, \ +0.26568491531006685413029428068906E-02, \ +0.32822721781658133778792170142517E-03, \ -0.23783447771430248269579807851050E-04, \ -0.11439804308100055514447076797047E-04, \ -0.14405943433238338455239717699323E-05, \ +0.52415956651148829963772818061664E-08, \ +0.38407306407844323480979203059716E-07, \ +0.85880244860267195879660515759344E-08, \ +0.10219226625855003286339969553911E-08, \ +0.21749132323289724542821339805992E-10, \ -0.22090238142623144809523503811741E-10, \ -0.63457533544928753294383622208801E-11, \ -0.10837746566857661115340539732919E-11, \ -0.11909822872222586730262200440277E-12, \ -0.28438682389265590299508766008661E-14, \ +0.25080327026686769668587195487546E-14, \ +0.78729641528559842431597726421265E-15, \ +0.15475066347785217148484334637329E-15, \ +0.22575322831665075055272608197290E-16, \ +0.22233352867266608760281380836693E-17, \ +0.16967819563544153513464194662399E-19, \ -0.57608316255947682105310087304533E-19, \ -0.17591235774646878055625369408853E-19, \ -0.36286056375103174394755328682666E-20, \ -0.59235569797328991652558143488000E-21, \ -0.76030380926310191114429136895999E-22, \ -0.62547843521711763842641428479999E-23, \ +0.25483360759307648606037606400000E-24, \ +0.25598615731739857020168874666666E-24, \ +0.71376239357899318800207052800000E-25, \ +0.14703759939567568181578956800000E-25, \ +0.25105524765386733555198634666666E-26, \ +0.35886666387790890886583637333333E-27, \ +0.39886035156771301763317759999999E-28, \ +0.21763676947356220478805333333333E-29, \ -0.46146998487618942367607466666666E-30, \ -0.20713517877481987707153066666666E-30, \ -0.51890378563534371596970666666666E-31 ] ) ae13cs = np.array ( [ \ -0.60577324664060345999319382737747, \ -0.11253524348366090030649768852718, \ +0.13432266247902779492487859329414E-01, \ -0.19268451873811457249246838991303E-02, \ +0.30911833772060318335586737475368E-03, \ -0.53564132129618418776393559795147E-04, \ +0.98278128802474923952491882717237E-05, \ -0.18853689849165182826902891938910E-05, \ +0.37494319356894735406964042190531E-06, \ -0.76823455870552639273733465680556E-07, \ +0.16143270567198777552956300060868E-07, \ -0.34668022114907354566309060226027E-08, \ +0.75875420919036277572889747054114E-09, \ -0.16886433329881412573514526636703E-09, \ +0.38145706749552265682804250927272E-10, \ -0.87330266324446292706851718272334E-11, \ +0.20236728645867960961794311064330E-11, \ -0.47413283039555834655210340820160E-12, \ +0.11221172048389864324731799928920E-12, \ -0.26804225434840309912826809093395E-13, \ +0.64578514417716530343580369067212E-14, \ -0.15682760501666478830305702849194E-14, \ +0.38367865399315404861821516441408E-15, \ -0.94517173027579130478871048932556E-16, \ +0.23434812288949573293896666439133E-16, \ -0.58458661580214714576123194419882E-17, \ +0.14666229867947778605873617419195E-17, \ -0.36993923476444472706592538274474E-18, \ +0.93790159936721242136014291817813E-19, \ -0.23893673221937873136308224087381E-19, \ +0.61150624629497608051934223837866E-20, \ -0.15718585327554025507719853288106E-20, \ +0.40572387285585397769519294491306E-21, \ -0.10514026554738034990566367122773E-21, \ +0.27349664930638667785806003131733E-22, \ -0.71401604080205796099355574271999E-23, \ +0.18705552432235079986756924211199E-23, \ -0.49167468166870480520478020949333E-24, \ +0.12964988119684031730916087125333E-24, \ -0.34292515688362864461623940437333E-25, \ +0.90972241643887034329104820906666E-26, \ -0.24202112314316856489934847999999E-26, \ +0.64563612934639510757670475093333E-27, \ -0.17269132735340541122315987626666E-27, \ +0.46308611659151500715194231466666E-28, \ -0.12448703637214131241755170133333E-28, \ +0.33544574090520678532907007999999E-29, \ -0.90598868521070774437543935999999E-30, \ +0.24524147051474238587273216000000E-30, \ -0.66528178733552062817107967999999E-31 ] ) ae14cs = np.array ( [ \ -0.1892918000753016825495679942820, \ -0.8648117855259871489968817056824E-01, \ +0.7224101543746594747021514839184E-02, \ -0.8097559457557386197159655610181E-03, \ +0.1099913443266138867179251157002E-03, \ -0.1717332998937767371495358814487E-04, \ +0.2985627514479283322825342495003E-05, \ -0.5659649145771930056560167267155E-06, \ +0.1152680839714140019226583501663E-06, \ -0.2495030440269338228842128765065E-07, \ +0.5692324201833754367039370368140E-08, \ -0.1359957664805600338490030939176E-08, \ +0.3384662888760884590184512925859E-09, \ -0.8737853904474681952350849316580E-10, \ +0.2331588663222659718612613400470E-10, \ -0.6411481049213785969753165196326E-11, \ +0.1812246980204816433384359484682E-11, \ -0.5253831761558460688819403840466E-12, \ +0.1559218272591925698855028609825E-12, \ -0.4729168297080398718476429369466E-13, \ +0.1463761864393243502076199493808E-13, \ -0.4617388988712924102232173623604E-14, \ +0.1482710348289369323789239660371E-14, \ -0.4841672496239229146973165734417E-15, \ +0.1606215575700290408116571966188E-15, \ -0.5408917538957170947895023784252E-16, \ +0.1847470159346897881370231402310E-16, \ -0.6395830792759094470500610425050E-17, \ +0.2242780721699759457250233276170E-17, \ -0.7961369173983947552744555308646E-18, \ +0.2859308111540197459808619929272E-18, \ -0.1038450244701137145900697137446E-18, \ +0.3812040607097975780866841008319E-19, \ -0.1413795417717200768717562723696E-19, \ +0.5295367865182740958305442594815E-20, \ -0.2002264245026825902137211131439E-20, \ +0.7640262751275196014736848610918E-21, \ -0.2941119006868787883311263523362E-21, \ +0.1141823539078927193037691483586E-21, \ -0.4469308475955298425247020718489E-22, \ +0.1763262410571750770630491408520E-22, \ -0.7009968187925902356351518262340E-23, \ +0.2807573556558378922287757507515E-23, \ -0.1132560944981086432141888891562E-23, \ +0.4600574684375017946156764233727E-24, \ -0.1881448598976133459864609148108E-24, \ +0.7744916111507730845444328478037E-25, \ -0.3208512760585368926702703826261E-25, \ +0.1337445542910839760619930421384E-25, \ -0.5608671881802217048894771735210E-26, \ +0.2365839716528537483710069473279E-26, \ -0.1003656195025305334065834526856E-26, \ +0.4281490878094161131286642556927E-27, \ -0.1836345261815318199691326958250E-27, \ +0.7917798231349540000097468678144E-28, \ -0.3431542358742220361025015775231E-28, \ +0.1494705493897103237475066008917E-28, \ -0.6542620279865705439739042420053E-29, \ +0.2877581395199171114340487353685E-29, \ -0.1271557211796024711027981200042E-29, \ +0.5644615555648722522388044622506E-30, \ -0.2516994994284095106080616830293E-30, \ +0.1127259818927510206370368804181E-30, \ -0.5069814875800460855562584719360E-31 ] ) e11cs = np.array ( [ \ -0.16113461655571494025720663927566180E+02, \ +0.77940727787426802769272245891741497E+01, \ -0.19554058188631419507127283812814491E+01, \ +0.37337293866277945611517190865690209, \ -0.56925031910929019385263892220051166E-01, \ +0.72110777696600918537847724812635813E-02, \ -0.78104901449841593997715184089064148E-03, \ +0.73880933562621681878974881366177858E-04, \ -0.62028618758082045134358133607909712E-05, \ +0.46816002303176735524405823868362657E-06, \ -0.32092888533298649524072553027228719E-07, \ +0.20151997487404533394826262213019548E-08, \ -0.11673686816697793105356271695015419E-09, \ +0.62762706672039943397788748379615573E-11, \ -0.31481541672275441045246781802393600E-12, \ +0.14799041744493474210894472251733333E-13, \ -0.65457091583979673774263401588053333E-15, \ +0.27336872223137291142508012748799999E-16, \ -0.10813524349754406876721727624533333E-17, \ +0.40628328040434303295300348586666666E-19, \ -0.14535539358960455858914372266666666E-20, \ +0.49632746181648636830198442666666666E-22, \ -0.16208612696636044604866560000000000E-23, \ +0.50721448038607422226431999999999999E-25, \ -0.15235811133372207813973333333333333E-26, \ +0.44001511256103618696533333333333333E-28, \ -0.12236141945416231594666666666666666E-29, \ +0.32809216661066001066666666666666666E-31, \ -0.84933452268306432000000000000000000E-33 ] ) e12cs = np.array ( [ \ -0.3739021479220279511668698204827E-01, \ +0.4272398606220957726049179176528E-01, \ -0.130318207984970054415392055219726, \ +0.144191240246988907341095893982137E-01, \ -0.134617078051068022116121527983553E-02, \ +0.107310292530637799976115850970073E-03, \ -0.742999951611943649610283062223163E-05, \ +0.453773256907537139386383211511827E-06, \ -0.247641721139060131846547423802912E-07, \ +0.122076581374590953700228167846102E-08, \ -0.548514148064092393821357398028261E-10, \ +0.226362142130078799293688162377002E-11, \ -0.863589727169800979404172916282240E-13, \ +0.306291553669332997581032894881279E-14, \ -0.101485718855944147557128906734933E-15, \ +0.315482174034069877546855328426666E-17, \ -0.923604240769240954484015923200000E-19, \ +0.255504267970814002440435029333333E-20, \ -0.669912805684566847217882453333333E-22, \ +0.166925405435387319431987199999999E-23, \ -0.396254925184379641856000000000000E-25, \ +0.898135896598511332010666666666666E-27, \ -0.194763366993016433322666666666666E-28, \ +0.404836019024630033066666666666666E-30, \ -0.807981567699845120000000000000000E-32 ] ) eta = 0.1 * r8_mach ( 3 ) ntae10 = r8_inits ( ae10cs, 50, eta ) ntae11 = r8_inits ( ae11cs, 60, eta ) ntae12 = r8_inits ( ae12cs, 41, eta ) nte11 = r8_inits ( e11cs, 29, eta ) nte12 = r8_inits ( e12cs, 25, eta ) ntae13 = r8_inits ( ae13cs, 50, eta ) ntae14 = r8_inits ( ae14cs, 64, eta ) xmax = - np.log ( r8_mach ( 1 ) ) xmax = xmax - np.log ( xmax ) if ( x <= - 32.0 ): value = np.exp ( - x ) / x * ( 1.0 \ + r8_csevl ( 64.0 / x + 1.0, ae10cs, ntae10 ) ) elif ( x <= - 8.0 ): value = np.exp ( - x ) / x * ( 1.0 \ + r8_csevl ( ( 64.0 / x + 5.0 ) / 3.0, ae11cs, ntae11 ) ) elif ( x <= - 4.0 ): value = np.exp ( - x ) / x * (1.0 \ + r8_csevl ( 16.0 / x + 3.0, ae12cs, ntae12 ) ) elif ( x <= - 1.0 ): value = - np.log ( - x ) \ + r8_csevl ( ( 2.0 * x + 5.0 ) / 3.0, e11cs, nte11 ) elif ( x == 0.0 ): print ( '' ) print ( 'R8_E1 - Fatal error!' ) print ( ' X is zero.' ) exit ( 'R8_E1 - Fatal error!' ) elif ( x <= 1.0 ): value = ( - np.log ( abs ( x ) ) - 0.6875 + x ) \ + r8_csevl ( x, e12cs, nte12 ) elif ( x <= 4.0 ): value = np.exp ( - x ) / x * ( 1.0 \ + r8_csevl ( ( 8.0 / x - 5.0 ) / 3.0, ae13cs, ntae13 ) ) elif ( x <= xmax ): value = np.exp ( - x ) / x * ( 1.0 \ + r8_csevl ( 8.0 / x - 1.0, ae14cs, ntae14 ) ) else: value = 0.0 return value def r8_e1_test ( ): #*****************************************************************************80 # #% R8_E1_TEST tests R8_E1. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 23 April 2016 # # Author: # # John Burkardt # import platform from e1_values import e1_values print ( '' ) print ( 'R8_E1_TEST:' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' R8_E1 evaluates the E1(x) exponential function.' ) print ( '' ) print ( ' X E1(X) R8_E1(X) Diff' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx1 = e1_values ( n_data ) if ( n_data == 0 ): break fx2 = r8_e1 ( x ) print ( ' %14.4f %14.6g %14.6g %14.6g' % ( x, fx1, fx2, abs ( fx1 - fx2 ) ) ) # # Terminate. # print ( '' ) print ( 'R8_E1_TEST:' ) print ( ' Normal end of execution.' ) return if ( __name__ == '__main__' ): from timestamp import timestamp timestamp ( ) r8_e1_test ( ) timestamp ( )