C ALGORITHM 794, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 25,NO. 2, June, 1999, P. 240--250. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Doc/ # Doc/PackingList # Doc/calgo.ps # Fortran77/ # Fortran77/Drivers/ # Fortran77/Drivers/Dp/ # Fortran77/Drivers/Dp/data # Fortran77/Drivers/Dp/driver.f # Fortran77/Drivers/Dp/res # Fortran77/Drivers/Dp/res.err # Fortran77/Src/ # Fortran77/Src/Dp/ # Fortran77/Src/Dp/specfun.f # Fortran77/Src/Dp/src.f # This archive created: Fri Oct 22 09:59:44 1999 export PATH; PATH=/bin:$PATH if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' if test -f 'PackingList' then echo shar: will not over-write existing file "'PackingList'" else cat << SHAR_EOF > 'PackingList' Algorithm:794 Language:Fortran77 Precision:Dp Src:specfun.f src.f SrcLibs:port Driver_0:driver.f @Src Data_0:stdin=data Res_0:stdout=res SHAR_EOF fi # end of overwriting check if test -f 'calgo.ps' then echo shar: will not over-write existing file "'calgo.ps'" else cat << SHAR_EOF > 'calgo.ps' %!PS-Adobe-2.0 %%Creator: dvips 5.54 Copyright 1986, 1994 Radical Eye Software %%Title: calgo.dvi %%CreationDate: Wed Feb 10 16:30:08 1999 %%Pages: 4 %%PageOrder: Ascend %%BoundingBox: 0 0 612 792 %%EndComments %DVIPSCommandLine: C:\HANKEL\TOMS\DVIPS.EXE calgo %DVIPSParameters: dpi=300, compressed, comments removed %DVIPSSource: TeX output 1999.02.10:1623 %%BeginProcSet: texc.pro /TeXDict 250 dict def TeXDict begin /N{def}def /B{bind def}N /S{exch}N /X{S N}B /TR{translate}N /isls false N /vsize 11 72 mul N /hsize 8.5 72 mul N /landplus90{false}def /@rigin{isls{[0 landplus90{1 -1}{-1 1} ifelse 0 0 0]concat}if 72 Resolution div 72 VResolution div neg scale isls{landplus90{VResolution 72 div vsize mul 0 exch}{Resolution -72 div hsize mul 0}ifelse TR}if Resolution VResolution vsize -72 div 1 add mul TR matrix currentmatrix dup dup 4 get round 4 exch put dup dup 5 get round 5 exch put setmatrix}N /@landscape{/isls true N}B /@manualfeed{ statusdict /manualfeed true put}B /@copies{/#copies X}B /FMat[1 0 0 -1 0 0]N /FBB[0 0 0 0]N /nn 0 N /IE 0 N /ctr 0 N /df-tail{/nn 8 dict N nn begin /FontType 3 N /FontMatrix fntrx N /FontBBox FBB N string /base X array /BitMaps X /BuildChar{CharBuilder}N /Encoding IE N end dup{/foo setfont}2 array copy cvx N load 0 nn put /ctr 0 N[}B /df{/sf 1 N /fntrx FMat N df-tail}B /dfs{div /sf X /fntrx[sf 0 0 sf neg 0 0]N df-tail}B /E{ pop nn dup definefont setfont}B /ch-width{ch-data dup length 5 sub get} B /ch-height{ch-data dup length 4 sub get}B /ch-xoff{128 ch-data dup length 3 sub get sub}B /ch-yoff{ch-data dup length 2 sub get 127 sub}B /ch-dx{ch-data dup length 1 sub get}B /ch-image{ch-data dup type /stringtype ne{ctr get /ctr ctr 1 add N}if}B /id 0 N /rw 0 N /rc 0 N /gp 0 N /cp 0 N /G 0 N /sf 0 N /CharBuilder{save 3 1 roll S dup /base get 2 index get S /BitMaps get S get /ch-data X pop /ctr 0 N ch-dx 0 ch-xoff ch-yoff ch-height sub ch-xoff ch-width add ch-yoff setcachedevice ch-width ch-height true[1 0 0 -1 -.1 ch-xoff sub ch-yoff .1 add]/id ch-image N /rw ch-width 7 add 8 idiv string N /rc 0 N /gp 0 N /cp 0 N{ rc 0 ne{rc 1 sub /rc X rw}{G}ifelse}imagemask restore}B /G{{id gp get /gp gp 1 add N dup 18 mod S 18 idiv pl S get exec}loop}B /adv{cp add /cp X}B /chg{rw cp id gp 4 index getinterval putinterval dup gp add /gp X adv}B /nd{/cp 0 N rw exit}B /lsh{rw cp 2 copy get dup 0 eq{pop 1}{dup 255 eq{pop 254}{dup dup add 255 and S 1 and or}ifelse}ifelse put 1 adv} B /rsh{rw cp 2 copy get dup 0 eq{pop 128}{dup 255 eq{pop 127}{dup 2 idiv S 128 and or}ifelse}ifelse put 1 adv}B /clr{rw cp 2 index string putinterval adv}B /set{rw cp fillstr 0 4 index getinterval putinterval adv}B /fillstr 18 string 0 1 17{2 copy 255 put pop}for N /pl[{adv 1 chg} {adv 1 chg nd}{1 add chg}{1 add chg nd}{adv lsh}{adv lsh nd}{adv rsh}{ adv rsh nd}{1 add adv}{/rc X nd}{1 add set}{1 add clr}{adv 2 chg}{adv 2 chg nd}{pop nd}]dup{bind pop}forall N /D{/cc X dup type /stringtype ne{] }if nn /base get cc ctr put nn /BitMaps get S ctr S sf 1 ne{dup dup length 1 sub dup 2 index S get sf div put}if put /ctr ctr 1 add N}B /I{ cc 1 add D}B /bop{userdict /bop-hook known{bop-hook}if /SI save N @rigin 0 0 moveto /V matrix currentmatrix dup 1 get dup mul exch 0 get dup mul add .99 lt{/QV}{/RV}ifelse load def pop pop}N /eop{SI restore showpage userdict /eop-hook known{eop-hook}if}N /@start{userdict /start-hook known{start-hook}if pop /VResolution X /Resolution X 1000 div /DVImag X /IE 256 array N 0 1 255{IE S 1 string dup 0 3 index put cvn put}for 65781.76 div /vsize X 65781.76 div /hsize X}N /p{show}N /RMat[1 0 0 -1 0 0]N /BDot 260 string N /rulex 0 N /ruley 0 N /v{/ruley X /rulex X V}B /V {}B /RV statusdict begin /product where{pop product dup length 7 ge{0 7 getinterval dup(Display)eq exch 0 4 getinterval(NeXT)eq or}{pop false} ifelse}{false}ifelse end{{gsave TR -.1 -.1 TR 1 1 scale rulex ruley false RMat{BDot}imagemask grestore}}{{gsave TR -.1 -.1 TR rulex ruley scale 1 1 false RMat{BDot}imagemask grestore}}ifelse B /QV{gsave newpath transform round exch round exch itransform moveto rulex 0 rlineto 0 ruley neg rlineto rulex neg 0 rlineto fill grestore}B /a{moveto}B /delta 0 N /tail{dup /delta X 0 rmoveto}B /M{S p delta add tail}B /b{S p tail} B /c{-4 M}B /d{-3 M}B /e{-2 M}B /f{-1 M}B /g{0 M}B /h{1 M}B /i{2 M}B /j{ 3 M}B /k{4 M}B /w{0 rmoveto}B /l{p -4 w}B /m{p -3 w}B /n{p -2 w}B /o{p -1 w}B /q{p 1 w}B /r{p 2 w}B /s{p 3 w}B /t{p 4 w}B /x{0 S rmoveto}B /y{ 3 2 roll p a}B /bos{/SS save N}B /eos{SS restore}B end %%EndProcSet TeXDict begin 40258431 52099146 1000 300 300 (C:\HANKEL\TOMS/calgo.dvi) @start /Fa 2 50 df0 D<000F131E393BC06180396060804038 403100D8801A1320130EA3130B3940118040903820C0C03930C07B80390F001E001B0D7E 8C21>49 D E /Fb 1 109 df<12701230A31260A412C0A312D012E0A2040E7E8D0A>108 D E /Fc 1 91 df90 D E /Fd 39 91 df<126012F0AA127012 00A4126012F0A212600414799312>33 D40 D<128012C0126012301218121C 120C120EA21207A7120EA2120C121C12181230126012C0128008197C9612>I<1207A3EA E738EAFFF8EA7FF0EA1FC0A2EA7FF0EAFFF8EAE738EA0700A30D0E7E9012>I<126012F0 12F8127812181230A212E012C00509798312>44 DI<126012F0 A212600404798312>I<13181338A21370A213E0A2EA01C0A3EA0380A2EA0700A2120EA2 5AA35AA25AA25AA25A0D1A7E9612>II<1206A2120E121E12FE12EE 120EACEAFFE0A20B147D9312>II<126012F0A212 601200A6126012F0A21260040E798D12>58 D<1318137813F8EA01E0EA07C0EA0F80EA1E 005A12F85A7E123C7EEA0F80EA07C0EA01E0EA00F8137813180D137E9312>60 DI<12C012F07E123C121FEA0F80EA 03C0EA01E0EA00F8137813F8EA01E0EA03C0EA0F80EA1F00123C12F85A12C00D137E9312 >I65 DIIIIIIIII76 DI< EAFEFEA2EA3E38123A123BA6EA39B8A6123813F812FEA20F147F9312>III82 DII<38FE3F80A238380E00AE6C5A6C5AEA07F06C5A1114809312>IIIIII E /Fe 1 2 df<127812FCA4127806067B9111> 1 D E /Ff 29 122 df<12E0A312601240A312C003087C820B>44 D<12E0A303037C820B>46 D<5A1207B4FCA21207B2EA7FF0A20C187D9713>49 DII<1378A213B8A21201EA0338A212071206120E12 1C121812381230127012E0B5FCA2EA0038A610187F9713>I57 D<12E0A31200AB12E0A303117C900B>I<133813 7CA2135C13CEA2EA01CF138F1387000313801307380703C0A21206380E01E0A2EA0FFF48 13F0EA1C004813F81478A248137C143C126000E0131E171A7F991A>65 D69 DI<00F013F0ABB5FCA2EAF000AD14 1A7D991B>72 D<00F0137814F0EB01E0EB03C0EB0780EB0F00131E5B5B5BEAF1E0EAF3F0 12F7EAFF78EAFE7CEAFC3C487E12F07F14801307EB03C014E01301EB00F014F8151A7D99 1B>75 D<12F0B3A6EAFFFEA20F1A7D9915>I<00FC1370A27EA212EFA2EAE780A2EAE3C0 A2EAE1E0A2EAE0F0A21378A2133CA2131EA2130FA2EB07F0A21303A2141A7D991B>78 D84 D<3AF003F00380A300F813703A7807780700 A21306393C0E3806EC3C0EA2EB0C1CD81E1C130CEC1E1C1318000EEB0E18120F9038380F 3813300007EB0730A213200003EB032013A001E013E013C000016D5A211A7F9924>87 D97 D<12E0A9EAE7C0EAFFF0EAF878EAF038EAE01C130C130EA513 1CA2EAF0381370EAFFE0EAE7C00F1A7D9914>I<130EA9EA07CEEA0FFEEA1C1EEA380E12 70A212E0A51270A2EA381EEA3C3EEA1FEEEA07CE0F1A7F9914>100 DI<12F0A41200A61270B1041B7E9A0A>105 D110 D I114 DI117 DI121 D E /Fg 39 90 df<13E01201EA0380EA0700120E5AA25AA25AA3 5AA91270A37EA27EA27E7EEA0380EA01E012000B217A9C16>40 D<12C07E12707E7E7EA2 7EA2EA0380A3EA01C0A9EA0380A3EA0700A2120EA25A5A5A5A5A0A217B9C16>II<1238127C127EA2123E120E121E121C127812F01260070B798416>44 D<127012F8A312700505788416>46 DII<12035AA25A5AB4FCA212E71207AEEAFFF8A30D197B 9816>III<137C13FC13DC1201EA039CA2EA07 1C120F120E121E123C1238127812F0B512E0A338001C00A53801FFC0A313197F9816>I< EA3FFE127FA20070C7FCA7EA77F0EA7FFC7FEA780FEA300738000380A2126012F0A238E0 0700EA781EEA3FFC6C5AEA07E011197E9816>I<127012F8A312701200A8127012F8A312 700512789116>58 D<1238127CA312381200A812381278127CA2123C121CA21238127012 F012400618799116>III<13E0487EA213B0A2EA03B8A31318EA071CA5EA 0E0EA2EA0FFEA2487EEA1C07A3387E0FC038FF1FE0387E0FC013197F9816>65 DI<3801F180EA07FBEA0FFFEA1F0FEA3C 07EA38031270A200F0C7FC5AA77E38700380A21238383C0700EA1F0FEA0FFE6C5AEA01F0 11197E9816>II<387FFFC0B5FC7EEA1C01A490C7 FCA2131CA2EA1FFCA3EA1C1CA290C7FC14E0A5EA7FFFB5FC7E13197F9816>I<387FFFE0 B5FC7EEA1C00A41400A2131CA2EA1FFCA3EA1C1CA290C7FCA6EA7F80487E6C5A13197F98 16>I<3801F180EA07FBEA0FFFEA1F0FEA3C07EA38031270A200F0C7FC5AA4EB1FC014E0 14C038F00380127013071238123CEA1E0FEA0FFFEA07FBEA01F313197F9816>I<387F07 F038FF8FF8387F07F0381C01C0A7EA1FFFA3EA1C01A9387F07F038FF8FF8387F07F01519 809816>II<387F0F E038FF8FF0387F0FE0381C0780EB0F00130E5B133C5B5B5BEA1DF0121F7F1338EA1E1C12 1C7FA27FA2EB0380387F07E038FF8FF0387F07E01419809816>75 DI<38FC07E0EAFE0FA2383A0B 80EA3B1BA513BBEA39B3A413F3EA38E3A21303A538FE0FE0A313197F9816>I<387E07F0 38FF0FF8387F07F0381D81C0A313C1121C13E1A213611371A313311339A21319131D130D A3EA7F07EAFF87EA7F031519809816>III< EA7FF0EAFFFC6C7EEA1C0FEB07801303A41307EB0F00EA1FFE5B7FEA1C0E7FA414101438 A2387F03F0EAFF83387F01E01519809816>82 DI<387FFFE0B5FCA2EAE0E0A400001300AFEA07FC487E6C5A13 197F9816>I<387F07F038FF8FF8387F07F0381C01C0B0380E0380A23807070013FF6C5A EA00F81519809816>I<38FE0FE0A338380380EA3C07001C1300A3EA1E0FEA0E0EA46C5A A4EA031813B8A3EA01B013F0A26C5A13197F9816>I<387E03F038FF07F8387E03F03838 00E0A4381C01C0A3137113F9A213D9A2000C1380A3EA0DDD138DA338078F00A213071519 809816>I<387F1F80EB3FC0EB1F80380E1E00131C12075BEA03B813F012015B12001201 7F120313B81207131C120FEA0E0EA2487E387E0FC038FF1FE0387E0FC013197F9816>I< 38FE0FE0EAFF1FEAFE0F381C0700A2EA0E0EA26C5AA3EA03B8A2EA01F0A26C5AA8EA03F8 487E6C5A13197F9816>I E /Fh 10 121 df23 D<124012E0124003037D820A>58 D<13201360A213C0A3EA0180A3EA0300A31206A25AA35AA35AA35AA35AA30B1D7E9511> 61 D 97 D<133C130C1318A41330EA07B0EA0C701210EA30601260A3EAC0C013C8A21241EA62 D0EA3C700E147E9311>100 DI<123C120C1218A41230A41260A412C012C8A312 D0126006147F930A>108 D110 D<1207EA1880EA19C0EA3180EA3800121E7EEA0380124112 E1EAC1001282127C0A0D7E8C10>115 D120 D E /Fi 4 54 df<120FEA30C0EA6060A2EA4020EAC030A9EA4020EA6060A2EA30C0EA0F000C137E9211> 48 D<120C121C12EC120CAFEAFFC00A137D9211>I<121FEA60C01360EAF07013301260EA 0070A2136013C012011380EA02005AEA08101210EA2020EA7FE012FF0C137E9211>I53 D E /Fj 10 121 df<007E1360000E13E0A3381C01C0A2EB 03801400485A130E130C5B485A5BEA71C00073C7FC12FC12F013127E9115>23 D<1310A3EB1F8013F03801CF00EA038048C7FC120EA6EA06FCEA0384EA06FC0008C7FC12 185A12201260A212E0A31270127CEA3F80EA1FE0EA03F8C67E131E130E130CEA0108EA00 F01125809C12>I<380FFFF85A5A386084001241EA81041201EA030CA212021206A2120E 120CEA1C0EA21238EA180615127E9118>I<126012F0A2126004047C830C>58 D<126012F0A212701210A41220A212401280040C7C830C>I<12C012F0123C120FEA03C0 EA00F01338130E6D7EEB01E0EB0078141EEC0780A2EC1E001478EB01E0EB0780010EC7FC 133813F0EA03C0000FC8FC123C12F012C0191A7D9620>62 D74 D100 D102 D<380787803808C8403810F0 C03820F1E0EBE3C03840E1803800E000A2485AA43863808012F3EB810012E5EA84C6EA78 7813127E9118>120 D E /Fk 12 118 df<1230127812F0126005047C830D>46 D97 D<13F8EA0704120CEA1802EA38041230EA7008EA7FF0EAE000A5EA6004 1308EA30101360EA0F800F127C9113>101 D104 DI107 DI110 D<13F8EA030CEA0E06487E121812 3000701380A238E00700A3130EA25BEA60185BEA30E0EA0F8011127C9115>I114 D<12035AA3120EA4EAFFE0EA1C00A35AA45AA4EAE080A2EAE100A2126612380B1A7C990E >116 D<381C0180EA2E03124EA2388E0700A2121CA2EA380EA438301C80A3EA383C3818 4D00EA0F8611127C9116>I E /Fl 64 123 df11 D<137E3801C180EA0301380703C0120EEB018090C7FCA5B512C0EA0E01B0387F87F8151D 809C17>II< 126012F012F812681208A31210A2122012401280050C7C9C0C>39 D<1380EA0100120212065AA25AA25AA35AA412E0AC1260A47EA37EA27EA27E12027EEA00 80092A7C9E10>I<7E12407E12307EA27EA27EA37EA41380AC1300A41206A35AA25AA25A 12205A5A092A7E9E10>I<1306ADB612E0A2D80006C7FCAD1B1C7E9720>43 D<126012F0A212701210A41220A212401280040C7C830C>II<12 6012F0A2126004047C830C>I48 D<5A1207123F12C71207B3A5EAFFF80D1C7C9B15>III<130CA2131C133CA2135C13DC139CEA011C12 0312021204120C1208121012301220124012C0B512C038001C00A73801FFC0121C7F9B15 >I61 D<1306A3130FA3EB1780A2EB37C01323 A2EB43E01341A2EB80F0A338010078A2EBFFF83802003CA3487FA2000C131F80001E5BB4 EBFFF01C1D7F9C1F>65 DI<90381F8080EBE0613801 801938070007000E13035A14015A00781300A2127000F01400A8007014801278A212386C EB0100A26C13026C5B380180083800E030EB1FC0191E7E9C1E>IIII<39FFF0FFF0390F 000F00AC90B5FCEB000FAD39FFF0FFF01C1C7F9B1F>72 DI<39FFF01FE0390F000780EC060014045C5C5C5C5C49C7FC13021306130FEB 17801327EB43C0EB81E013016D7E1478A280143E141E80158015C039FFF03FF01C1C7F9B 20>75 DIIIII82 D<3807E080EA1C19EA30051303EA600112E01300A36C13007E127CEA7FC0EA3FF8EA1FFE EA07FFC61380130FEB07C0130313011280A300C01380A238E00300EAD002EACC0CEA83F8 121E7E9C17>I<007FB512C038700F010060130000401440A200C014201280A300001400 B1497E3803FFFC1B1C7F9B1E>I<39FFF01FF0390F000380EC0100B3A26C130213800003 5BEA01C03800E018EB7060EB0F801C1D7F9B1F>I<39FFE00FF0391F0003C0EC01806C14 00A238078002A213C000035BA2EBE00C00011308A26C6C5AA213F8EB7820A26D5AA36D5A A2131F6DC7FCA21306A31C1D7F9B1F>I<3AFFE1FFC0FF3A1F003E003C001E013C13186C 6D1310A32607801F1320A33A03C0278040A33A01E043C080A33A00F081E100A39038F900 F3017913F2A2017E137E013E137CA2013C133C011C1338A20118131801081310281D7F9B 2B>I<39FFF07FC0390FC01E003807800CEBC00800035B6C6C5A13F000005BEB7880137C 013DC7FC133E7F7F80A2EB13C0EB23E01321EB40F0497E14783801007C00027F141E0006 131F001F148039FF807FF01C1C7F9B1F>I<39FFF003FC390F8001E00007EB00C06D1380 0003EB01006D5A000113026C6C5A13F8EB7808EB7C18EB3C10EB3E20131F6D5A14C06D5A ABEB7FF81E1C809B1F>I97 D<12FC121CAA137CEA1D87 381E0180381C00C014E014601470A6146014E014C0381E018038190700EA10FC141D7F9C 17>IIII< 13F8EA018CEA071E1206EA0E0C1300A6EAFFE0EA0E00B0EA7FE00F1D809C0D>II<12FC121CAA137C1387EA1D03001E1380121CAD38FF9FF0141D7F9C17>I<1218 123CA21218C7FCA712FC121CB0EAFF80091D7F9C0C>I<13C0EA01E0A2EA00C01300A7EA 07E01200B3A21260EAF0C012F1EA6180EA3E000B25839C0D>I<12FC121CAAEB0FE0EB07 80EB06005B13105B5B13E0121DEA1E70EA1C781338133C131C7F130F148038FF9FE0131D 7F9C16>I<12FC121CB3A9EAFF80091D7F9C0C>I<39FC7E07E0391C838838391D01901800 1EEBE01C001C13C0AD3AFF8FF8FF8021127F9124>IIII<3803E080EA0E19EA1805EA3807EA7003A212E0A61270A2EA38071218EA0E 1BEA03E3EA0003A7EB1FF0141A7F9116>III<1204A4120CA2121C123C EAFFE0EA1C00A91310A5120CEA0E20EA03C00C1A7F9910>I<38FC1F80EA1C03AD130712 0CEA0E1B3803E3F014127F9117>I<38FF07E0383C0380381C0100A2EA0E02A2EA0F06EA 0704A2EA0388A213C8EA01D0A2EA00E0A3134013127F9116>I<39FF3FC7E0393C0703C0 001CEB01801500130B000E1382A21311000713C4A213203803A0E8A2EBC06800011370A2 EB8030000013201B127F911E>I<38FF0FE0381E0700EA1C06EA0E046C5AEA039013B0EA 01E012007F12011338EA021C1204EA0C0E487E003C138038FE1FF014127F9116>I<38FF 07E0383C0380381C0100A2EA0E02A2EA0F06EA0704A2EA0388A213C8EA01D0A2EA00E0A3 1340A25BA212F000F1C7FC12F312661238131A7F9116>II E /Fm 12 118 df67 D97 D102 D<1238127CA312381200A412FCA2123CAB12FFA208187F97 0B>105 D<39FC7E0FC039FD8730E0393E07C0F0A2003C1380A939FF1FE3FCA21E0F7E8E 23>109 DIII114 DI<120CA4121CA2EA3FE012FFEA3C00A71330A4EA1E60EA07 C00C157F9410>II E /Fn 45 122 df<126012F0A2126004047D830A>46 D<1206120E12FE120EB1EAFFE00B 157D9412>49 DI<126012F0A21260 1200A6126012F0A21260040E7D8D0A>58 D<13101338A3135CA3138EA3EA0107A2380203 80A33807FFC0EA0401A2380800E0A2001813F0123838FE03FE17177F961A>65 DIIIIII<38FF83FE381C0070AA381FFFF0 381C0070AA38FF83FE17177F961A>II<38FF80 FE381C0078146014401480EB0100130613085B13381378139CEA1D0E121EEA1C07EB0380 EB01C0A2EB00E014701478147C38FF80FF18177F961B>75 D<00FC13FE001E1338001F13 101217EA1380EA11C0A2EA10E013701338A2131C130E130F1307EB0390EB01D0A2EB00F0 14701430123800FE131017177F961A>78 D<13FCEA0303380E01C0381C00E04813700030 1330007013380060131800E0131CA700701338A200301330003813706C13E0380E01C038 030300EA00FC16177E961B>II82 DI<387FFFF8386038180040 1308A200801304A300001300AF3803FF8016177F9619>I<38FF80FE381C00381410B06C 132012066C13403801818038007E0017177F961A>I<3AFF07FC3F803A3C00E00E00001C 1404A2EB0170000E5CA2EB023800075CA2EB041CD803845BA2EB880ED801C85BA2EBD007 D800F05BA3EBE003016090C7FCA221177F9624>87 D<12FCA212C0B3AB12FCA206217D98 0A>91 D<12FCA2120CB3AB12FCA2062180980A>93 D97 D<12F81238A8EA39F0EA3E0CEA380613077F1480A414005B1306EA361CEA21F011177F96 14>II<133E130EA8EA07CEEA1C3EEA300E1270126012E0A412601270EA301EEA182E38 07CF8011177F9614>IIII<12F812 38A813F8EA3B1CEA3C0E1238AA38FE3F8011177F9614>I<12301278A212301200A512F8 1238AC12FE07177F960A>I<1203EA0780A2EA0300C7FCA5EA1F801203AF1243EAE30012 E7127C091D82960B>I<12F81238A8133E13381330134013801239EA3FC0EA39E0123813 F01378133CA2EAFE7F10177F9613>I<12F81238B3A312FE07177F960A>I<38F8F83E383B 1CC7393C0F0380EA380EAA39FE3F8FE01B0E7F8D1E>IIII 114 DI<1208A31218A21238EAFFC0EA3800A71340A4EA1C80EA0F000A14 7F930E>II121 D E /Fo 4 73 df0 D24 DI<017E130648B4130CD8061F131C486C133C003814381230D8600E137800C0 1470EA001E15F015E0131C1401013C13C048B5FC5A38003803017813801370140713F013 E015001201495A1502D803801306159CD8070013F00004EB07C01F1E809B23>72 D E /Fp 35 122 df<12E0A312601240A312C003087C820C>44 DI<12E0A303037C820C>I<130113031306A3130CA31318A31330A31360A213C0A3EA0180 A3EA0300A31206A25AA35AA35AA35AA35AA210297E9E15>I<12E0A31200AC12E0A30312 7C910C>58 D66 D68 D70 DI<00FCEB07E0A300EE130DA300E71319A3EB803900E3 1331EBC071A200E11361A2EBE0E1A200E013C113F1EB7181A3EB3B01A3131EA313001B1D 7C9C24>77 D83 DI<00F01370B3A5 007813E0A2383C01C0381E0380EA0F073807FE00EA01F8141E7C9C1D>I<00F001F81370 A3007801B81360D9019C13E0A3D83C03EB01C0141E140E001E15809038070F0313061407 000F1500010E1387130C0007EB0386A2019C138E019813CE0003EB01CCA339019000C801 D013D801F013F85B00001470241D7F9C27>87 D97 D99 D<1307ABEA07C7EA1FF7EA3FFFEA3C1FEA7807127012E0A61270EA 780FEA3C1FEA3FFFEA1FF7EA07C7101D7F9C15>II<13FC 12011203EA0700120EA7EAFFE0A2EA0E00B00E1D809C0D>I<3807C3C0EA0FFF5A383838 00487EA56C5AEA3FF05BEA77C00070C7FCA2EA3FFC13FF481380EA700738E001C0A3EAF0 03387C0F80383FFF006C5AEA07F8121B7F9115>I<12E0ABEAE3E0EAEFF0EAFFF8EAF83C EAF01C12E0AD0E1D7D9C15>I<12F0A41200A71270B2041D7E9C0A>I<12E0AB133C137813 F0EAE1E0EAE3C0EAE780EAEF00B4FC138012FBEAF9C0EAF1E012E013F013781338133C13 1E0F1D7D9C14>107 D<12E0B3AB031D7D9C0A>I<38E3F03F39EFF8FF80D8FFFD13C039F8 1F81E038F00F00EAE00EAD1B127D9122>III I114 DI<121CA6EAFFE0A2EA1C00AC1320EA1FF0120FEA07C00C187F970F>III<39E03E0380A3D870371300EB7707 A213733838E38EA33818E18C381CC1CC001D13DCA2380D80D8000F13F8A2000713701912 7F911C>I121 D E /Fq 31 122 df<12FCA61200B312FCA6061F7A9E13>58 D<147EA214FFA3903801EF80 A3903803CFC014C7A290380787E01483010F7FA21403496C7EA2131E90383E00FCA2133C 017C137EA2137801F87FA25B0001EC1F80A25B48B612C0A34815E09038C000075B000FEC 03F0A248C713F81501A2003E15FC1500A24815FE167EA248157F163F28327EB12D>65 D68 DII<00FCEC07E0B3A4B7FCA400FCC71207B3A623327AB130>72 D<12FCB3B3AE06327AB113>I<00FCEC07F0ED0FE016C0ED1F80ED3F00157E5D4A5A4A5A 4A5A140F4A5A5D4AC7FC147E5C495A495A495A130F131F497EA2497E13FD38FDF9FC38FF F0FEEBE07EEBC07F8001807F496C7E48130F48801407816E7E1401816E7E157E157F8116 80ED1FC0150F16E0150716F0ED03F825327AB12F>75 D<12FCB3B3A9B612C0A51A327AB1 24>I78 D80 D84 D<13FF000713C0001F13E04813F0EB01 F8383800FC0030137C0020137EC7123EA5EB07FE13FF1203120F381FE03EEA3F00127C12 FC5AA36C137E14FEEA7F0313FF6C13BE381FFE3EEA07F0171F7D9E20>97 D<12F8B3EB1F80EBFFE000FB7FB57EEB81FC38FE007E487F487F1580140FA2EC07C0A814 0F1580A2141F15006C133E6C13FE38FF03FCEBFFF800FB5B00F813C0013FC7FC1A327AB1 23>II101 DI<017F13F83901FFC7FC4813FF5A390FC1FC00EA1F80EB007CA2 003E7FA66C5BA2EB80FC380FC1F8EBFFF0485B001D5BD81C7FC7FC003CC8FCA37E381FFF F814FF6C14804814C04814E0393E000FF048130300FCEB01F8481300A46C1301007EEB03 F06CEB07E0EBE03F000FB512806C1400000113FC38003FE01E2E7E9E22>I<12F8B3EB3F 80EBFFE000F913F000FB13F8EAFF8313004813FC48137CA35AB3A316327AB123>I<12FC A61200AB127CB3AD06307CAF10>I<12F8B3147F14FE495A495A495A495A495A495A91C7 FC137E5B12F912FB12FF7F7F13BFEB1F80486C7E12FC486C7E6D7EA26D7E6D7EA2147E80 A2EC1F80EC0FC01A327BB121>107 D<12F8B3B3AE05327BB110>I<3AF81FC007F090397F F01FFC3AF9FFF87FFE00FB14FF3AFF81FDE07F9039007FC01F48028013804890383F000F A348133EB3A3291F7A9E36>I<38F83F80EBFFE000F913F000FB13F8EAFF8313004813FC 48137CA35AB3A3161F7A9E23>II<38F8 1F80EBFFE000FB7FB57EEB83FC48C67E48133F5AEC1F80140FA215C01407A7140F1580A2 141FEC3F006C137E6C13FE38FF03FCEBFFF800FB5B00F813C0013FC7FC90C8FCAE1A2D7A 9E23>I114 D<48B4FC000F13E04813F85AEA7E01387C0030481300A47E127EEA7FE0EA3FFE381FFF80 6C13E0000313F038003FF81303EB00FC147CA31240126000F813F8EAFE03B512F06C13E0 001F13C03801FE00161F7E9E1A>II<00F8137CB3A514FCA21301EAFC07EA7FFFEBFE 7CEA3FFCEA0FE0161F7A9E23>I<00F8EB01F06C1303007C14E0127E003EEB07C0A27EEC 0F801380120FEC1F00EA07C0A2143EEA03E0143C0001137C13F01478000013F813F8EB78 F0A2EB79E0133DA2EB1DC0131FA26D5AA391C7FC5B131EA2133E133CA2EA2078EA38F8EA 3FF0A25BEA0F801C2D7F9E1F>121 D E end %%EndProlog %%BeginSetup %%Feature: *Resolution 300dpi TeXDict begin %%PaperSize: Letter %%EndSetup %%Page: 1 1 1 0 bop 225 191 a Fq(Numerical)21 b(Hank)n(el)h(transfo)n(rm)g(b)n(y)g (the)h(F)n(o)n(rtran)g(p)n(rogram)225 274 y(HANKEL)225 357 y(P)n(a)n(rt)g(I)r(I:)f(T)-6 b(echnical)22 b(Details)225 472 y Fp(Thomas)13 b(Wieder)225 530 y(Da)o(rmstadt)g(Universit)o(y)h (of)g(T)m(echnology)m(,)f(Germany)225 588 y(FB)g(Materialwissenschaft,) j(F)o(G)d(Strukturfo)o(rschung)225 646 y(http://www.tu-da)o (rmstadt.de/)p Fo(\030)p Fp(wieder)p 225 714 1495 1 v 225 772 a Fn(Categories)d(and)g(Sub)r(ject)g(Descriptors:)j(F.2.1)d([)p Fm(Computation)i(of)i(transforms)p Fn(])225 831 y(General)c(T)m(erms:)k (Numerical)c(analysis)225 889 y(Additional)f(Key)j(W)m(ords)f(and)f (Phrases:)k(Hank)o(el)d(transform)p 225 939 V 267 1084 a Fl(Some)f(input)h(and)g(output)g(is)g(done)h(via)e(\014les.)18 b(If)11 b(sampled)f(data)h(whic)o(h)g(are)g(stored)i(in)d(a)h(\014le) 225 1134 y(are)j(to)g(b)q(e)g(transformed,)f(then)h(these)h(data)f(m)o (ust)e(b)q(e)j(read)f(from)e(this)i(\014le)f(\(e.g.)g Fk(hankel.in)p Fl(\).)225 1184 y(This)i(input)g(\014le)h(will)e(b)q(e)i (op)q(ened)g(b)o(y)f(DHHANKEL)h(\(DHHANKEL)g(will)e(ask)h(for)g(the)h (\014le)225 1233 y(name\).)21 b(Results)16 b(will)d(b)q(e)j(written)g (to)f(an)g(output)h(\014le)f(\(e.g.)f Fk(hankel.out)p Fl(\))i(and)f(a)g(rep)q(ort)i(on)225 1283 y(the)i(calculation)f(will)g (b)q(e)h(written)g(to)g(an)g(error)g(log\014le)f(\(e.g)h Fk(hankel.err)p Fl(\).)32 b(Both)19 b(output)225 1333 y(\014les)13 b(m)o(ust)g(b)q(e)g(op)q(ened)h(b)o(y)f(the)h(user's)g (calling)e(program)f(b)q(efore)j(the)g(call)e(to)h(DHHANKEL.)225 1383 y(A)k(call)f(to)h(subroutine)h(DHINIT)f(m)o(ust)e(preceede)20 b(an)o(y)c(call)h(to)f(HANKEL.)i(By)f(DHINIT,)225 1433 y(the)f(logical)e(n)o(um)o(b)q(ers)h(\(LUIN,)g(LUERR,)f(LUOUT\))i(of)f (the)h(\014les)g(will)e(b)q(e)i(set.)23 b(If)15 b(the)h(user)225 1483 y(w)o(an)o(ts)f(to)g(ha)o(v)o(e)h(all)e(output)h(to)g(the)h (standard)g(output)g(\(i.e.)22 b(the)16 b(terminal\),)d(then)j(the)g (call)225 1532 y(to)e(DHINIT)g(w)o(ould)f(read)h(as)g(CALL)g (DHINIT\(10,1,1\).)267 1582 y(The)20 b(user)g(has)g(either)g(to)g(pro)o (vide)f(a)h(SUBR)o(OUTINE)g(DHFUNC)g(whic)o(h)f(giv)o(es)h Fj(f)t Fl(\()p Fj(x)p Fl(\))225 1632 y(as)91 b(a)g(function)g(of)f Fj(x)h Fl(or)g(to)g(pro)o(vide)g(NEND)h(data)225 1682 y(pairs)18 b(\()p Fj(x)371 1688 y Fi(1)389 1682 y Fj(;)7 b(f)t Fl(\()p Fj(x)472 1688 y Fi(1)491 1682 y Fl(\)\))p Fj(;)g(:)g(:)g(:)e(;)i Fl(\()p Fj(x)656 1688 y Fh(n)678 1682 y Fj(;)g(f)t Fl(\()p Fj(x)761 1688 y Fh(n)783 1682 y Fl(\)\))p Fj(;)g(:)g(:)g(:)f(;)h Fl(\()p Fj(x)949 1688 y Fh(nend)1024 1682 y Fj(;)g(f)t Fl(\()p Fj(x)1107 1688 y Fh(nend)1183 1682 y Fl(\)\))18 b(from)e(whic)o(h)h(HANKEL)i(will)225 1732 y(in)o(terp)q(olate)14 b Fj(f)t Fl(\()p Fj(x)p Fl(\).)20 b(F)m(or)14 b Fj(f)t Fl(\()p Fj(x)p Fl(\))h(giv)o(en)f(in)g(DHFUNC)g (\(this)h(case)g(is)f(designated)h(b)o(y)f(ID)o(A)m(T=1)225 1782 y(within)f(HANKEL\))i(the)f(call)g(to)f(HANKEL)i(is)225 1858 y Fg(CALL)21 b(DHINIT\(10,11,12\))225 1908 y(CALL)g(HANKEL)g (\(XI,NU,IVERBOSE)o(,IERR)o(,HTFU)o(N\))225 1984 y Fl(where)14 b(XI)f(corresp)q(onds)i(to)e Fj(\030)r Fl(,)f(NU)h(means)f Fj(\027)s Fl(,)g(IVERBOSE)i(is)f(set)g(equal)g(to)g(0,)f(1,)g(or)h(2,)g (and)225 2034 y(IERR)j(is)h(an)f(error)i(\015ag.)26 b(F)m(or)17 b(IVERBOSE)g(=)g(0)g(no)f(output)h(will)f(b)q(e)h(written)g(\(ho)o(w)o (ev)o(er,)225 2084 y(b)q(oth)j(the)h(error)g(log\014le)e(and)h(the)h (output)g(\014le)f(m)o(ust)f(ha)o(v)o(e)h(b)q(een)h(op)q(ened)g(an)o (yw)o(a)o(y\),)g(for)225 2134 y(IVERBOSE)e(=)f(1)g(some)f(information)e (on)i(the)i(calculation)e(pro)q(cedure)i(will)e(b)q(e)h(prin)o(ted,)225 2184 y(and)d(IVERBOSE)i(=)f(2)f(is)g(for)g(debugging)g(purp)q(oses)i (only)m(.)22 b(If)15 b(IERR)g(is)h(equal)f(to)g(zero)i(on)225 2233 y(return,)d(then)g(no)f(error)i(has)e(b)q(een)i(detected.)20 b(HTFUN)14 b(is)f(the)h(desired)g(result)h Fo(H)1536 2239 y Fh(\027)1556 2233 y Fl(\()p Fj(\030)r(;)7 b(f)t Fl(\()p Fj(x)p Fl(\)\).)225 2283 y(The)14 b(structure)i(of)e(DHFUNC)g (is)f(displa)o(y)o(ed)h(in)f(table)h(1.)267 2333 y(T)m(abulated)g(data) h(can)g(b)q(e)h(passed)g(to)f(HANKEL)h(either)g(via)f(an)g(input)g (\014le)g(\(this)g(case)h(is)225 2383 y(designated)g(b)o(y)e(ID)o(A)m (T=2\))h(or)g(from)e(the)i(calling)f(program)f(\(i.e.)21 b(the)16 b(user's)g(program\))d(via)225 2433 y(the)g(call)e(to)h (HANKEL)h(\(this)f(case)h(is)f(designated)h(b)o(y)e(ID)o(A)m(T=3\).)17 b(If)12 b(the)h(data)e(are)i(stored)g(in)225 2483 y(an)g(input)h (\014le)f(\(ID)o(A)m(T=2\),)g(then)h(the)g(user)h(has)f(to)f(call)g (SUBR)o(OUTINE)h(DHFUNK2)g(once)225 2532 y(\(and)g(only)f(once\))i(in)e (adv)n(ance)h(of)f(the)i(call)e(to)h(HANKEL.)g(The)g(calling)f (sequence)j(will)c(b)q(e)p eop %%Page: 2 2 2 1 bop 225 125 a Ff(2)71 b Fe(\001)78 b Ff(T.)13 b(Wieder)547 398 y Fn(T)m(able)e(1.)35 b(The)11 b(structure)f(of)h(SUBR)o(OUTINE)i (DHFUNC.)331 464 y Fd(SUBROUTINE)h(DHFUNC\(NU,)o(X,Y)o(,XL)o(IM)o(IT,)o (IER)o(R\))331 506 y(IMPLICIT)h(NONE)331 547 y(INTEGER)g(IERR,LUIN,)o (LUO)o(UT,)o(LUE)o(RR)331 589 y(DOUBLE)g(PRECISION)g(X,Y,NU,XLI)o(MI)o (T)331 630 y(COMMON)g(/INOUTE/)g(LUIN,LUOUT,)o(LU)o(ERR)225 672 y(C)225 713 y(C)88 b(**********)o(***)o(**)o(***)o(***)o(***)o(***) o(**)o(***)o(***)o(***)o(***)o(**)o(***)o(***)o(***)o(***)o(**)o(***)o (***)225 755 y(C)225 796 y(C)g(INPUT)16 b(VARIABLES)o(:)f(NU,)h(X,)h (LUIN,LUOUT,)o(LUE)o(RR)225 838 y(C)88 b(OUPUT)16 b(VARIABLES)o(:)f(Y,) i(XLIMIT,)e(IERR)225 879 y(C)225 921 y(C)88 b(NU)17 b(=)g(ORDER)f(OF)h (TRANSFORM)d(WITH)i(-1/2)h(<)g(NU)225 962 y(C)88 b(X)35 b(=)17 b(PARAMETER)d(OF)j(TRANSFORM)e(WITH)h(0)h(<)h(X)225 1004 y(C)88 b(LUIN)16 b(=)h(LOGICAL)e(NUMBER)h(OF)h(INPUT)f(FILE)225 1045 y(C)211 b(\(WILL)16 b(HAVE)g(BEEN)h(SET)f(BY)h(PREVIOUS)e(CALL)h (TO)h(SUBROUTINE)d(DHINIT\))225 1087 y(C)88 b(LUOUT)16 b(=)h(LOGICAL)e(NUMBER)h(OF)h(OUTPUT)e(FILE)225 1128 y(C)211 b(\(WILL)16 b(HAVE)g(BEEN)h(SET)f(BY)h(PREVIOUS)e(CALL)h(TO)h (SUBROUTINE)d(DHINIT\))225 1170 y(C)88 b(LUERR)16 b(=)h(LOGICAL)e (NUMBER)h(OF)h(ERROR)e(LOG)i(FILE)225 1211 y(C)211 b(\(WILL)16 b(HAVE)g(BEEN)h(SET)f(BY)h(PREVIOUS)e(CALL)h(TO)h(SUBROUTINE)d (DHINIT\))225 1253 y(C)88 b(IERR)16 b(=)h(ERROR)f(FLAG,)g(IERR)g(=)i(0) f(==>)f(NO)h(ERROR)f(HAS)h(BEEN)f(DETECTED)225 1295 y(C)225 1336 y(C)88 b(**********)o(***)o(**)o(***)o(***)o(***)o(***)o(**)o(***) o(***)o(***)o(***)o(**)o(***)o(***)o(***)o(***)o(**)o(***)o(***)225 1378 y(C)225 1419 y(C)g(INITIALIZE)14 b(THE)i(ERROR)g(FLAG.)225 1461 y(C)88 b(***)16 b(IERR=0)g(***)g(MEANS)g(THAT)h(NO)f(ERRORS)g (OCCURRED)f(SO)i(FAR.)225 1502 y(C)331 1544 y(IERR=0)225 1585 y(C)225 1627 y(C)88 b(THE)16 b(INDEFINITE)e(INTEGRAL)h(OVER)h(THE) h(FUNCTION)e(***)h(FUN)h(***)225 1668 y(C)88 b(WILL)16 b(BE)h(SPLIT)f(INTO)g(A)h(DEFINITE)e(INTEGRAL)g(FROM)h(0)h(TO)225 1710 y(C)88 b(***)16 b(XLIMIT)g(***)g(AND)h(AN)g(INDEFINITE)d(INTEGRAL) h(FROM)225 1751 y(C)88 b(***)16 b(XLIMIT)g(***)g(TO)h(INFINITY.)225 1793 y(C)88 b(***)16 b(XLIMIT)g(***)g(IS)h(THE)g(UPPER)f(LIMIT)g(OF)h (THE)f(DEFINITE)f(INTEGRAL.)225 1834 y(C)88 b(YOU)16 b(NEED)h(NOT)f(TO)h(SPECIFY)e(***)i(XLIMIT)e(***.)i(IN)f(THAT)h(CASE)f (YOU)225 1876 y(C)88 b(JUST)16 b(GIVE)g(A)i(NEGATIVE)c(VALUE,)i(E.G.)g (***)h(XLIMIT=-1)o(.0D)o(0)e(***)225 1917 y(C)331 1959 y(XLIMIT=-1.)o(0D0)225 2000 y(C)225 2042 y(C)88 b(HERE)16 b(YOU)h(GIVE)f(YOUR)g(FUNCTION)f(***)h(FUN)h(***)225 2083 y(C)88 b(Y=FUN\(X,NU)o(\))225 2125 y(C)331 2166 y(Y=...)225 2208 y(C)225 2249 y(C)g(REMEMBER:)14 b(IF)j(AN)g(ERROR)f (HAS)h(OCCURRED,)d(THEN)i(PUT)h(***)f(IERR=1)g(***)g(!)225 2291 y(C)331 2332 y(RETURN)331 2374 y(END)p eop %%Page: 3 3 3 2 bop 994 125 a Ff(HANKEL,)12 b(version:)k(F)o(eb)o(rua)o(ry)c(1999) 80 b Fe(\001)70 b Ff(3)356 233 y Fg(CALL)21 b(DHINIT\(10,11,12\))225 283 y(C)356 332 y(DUMMY=0.0D0)356 382 y(IERR4=0)225 432 y(C)356 482 y(CALL)g(DHFUNK2)f(\(DUMMY,DUMMY,IERR4)o(\))225 532 y(C)225 581 y(C)109 b(SET)21 b(INPUT)g(VARIBALES:)e (XI,NU,IVERBOSE,)g(E.G.:)356 631 y(XI=1.0D0)356 681 y(NU=5.0D0)356 731 y(IVERBOSE=1)225 781 y(C)356 831 y(CALL)i(HANKEL)f (\(XI,NU,IVERBOSE,IER)o(R,HTF)o(UN\))225 880 y(C)225 930 y(C)109 b(REPEATED)20 b(CALLS)h(TO)g(HANKEL)g(ARE)g(POSSIBLE)f (WITHOUT)g(NEW)h(CALL)g(TO)h(DHFUNK2)225 999 y Fl(If)9 b(the)h(data)g(are)g(stored)g(in)g(arra)o(ys)f(XD)o(A)m(T,)g(YD)o(A)m (T)g(and)g(will)f(b)q(e)i(transfered)i(via)c(a)i(subroutine)225 1049 y(call)g(\(ID)o(A)m(T=3\),)h(then)h(the)f(user)i(has)e(to)g(call)f (SUBR)o(OUTINE)i(DHFUNK3)f(once)h(\(and)f(only)225 1099 y(once\))k(in)e(adv)n(ance)h(of)f(the)i(call)e(to)h(HANKEL.)g(The)g (calling)f(sequence)j(will)c(b)q(e)356 1168 y Fg(CALL)21 b(DHINIT\(10,11,12\))225 1218 y(C)356 1268 y(DUMMY=0.0D0)356 1318 y(IERR4=0)225 1367 y(C)109 b(AS)21 b(AN)h(EXAMPLE,)e(WE)h(PUT)g(Y) h(=)f(X**0.5)g(FOR)g(X)h(<=)f(1;)g(0)h(ELSEWHRE:)356 1417 y(NEND=100)356 1467 y(DO)f(10)h(N=1,NEND)356 1517 y(XDAT\(N\)=DBLE\(N\)/)o(DBLE\()o(NEND\))356 1567 y(IF)f (\(XDAT\(N\).LE.1.0D0\))d(THEN)356 1616 y(YDAT\(N\)=XDAT\(N\)*)o(*0.5D) o(0)356 1666 y(ELSE)j(IF)g(\(XDAT\(N\).GT.1.0D0\))d(THEN)356 1716 y(Y=0.0D0)356 1766 y(END)j(IF)225 1816 y(C)356 1866 y(CALL)g(DHFUNK3)f(\(DUMMY,DUMMY,XDAT,)o(YDAT,)o(NEND)o(,IERR)o(4\))225 1915 y(C)225 1965 y(C)109 b(SET)21 b(INPUT)g(VARIBALES:)e (XI,NU,IVERBOSE,)g(E.G.:)356 2015 y(XI=1.0D0)356 2065 y(NU=5.0D0)356 2115 y(IVERBOSE=1)225 2164 y(C)356 2214 y(CALL)i(HANKEL)f(\(XI,NU,IVERBOSE,IER)o(R,HTF)o(UN\))225 2264 y(C)225 2314 y(C)109 b(REPEATED)20 b(CALLS)h(TO)g(HANKEL)g(ARE)g (POSSIBLE)f(WITHOUT)g(NEW)h(CALL)g(TO)h(DHFUNK3)267 2383 y Fl(Success)17 b(or)f(failure)f(of)f(HANKEL)j(for)e(a)g(giv)o(en)g Fj(f)t Fl(\()p Fj(x)p Fl(\))i(will)d(dep)q(end)i(on)g(the)g(c)o(hoice)g (of)f Fj(x)1707 2389 y Fh(l)225 2433 y Fl(and)i(the)h(implemen)o (tation)c(of)j(DHFUNC,)g(but)g(if)g(no)g(clue)h(on)f Fj(x)1302 2439 y Fh(l)1332 2433 y Fl(is)g(a)o(v)n(ailable,)f(then)i (the)225 2483 y(user)j(puts)f(XLIMIT=-1.0.)35 b(Then)20 b(HANKEL)g(will)f(set)h Fj(x)1211 2489 y Fh(l)1243 2483 y Fl(equal)g(to)f Fj(x)1439 2489 y Fh(as)1494 2483 y Fl(whic)o(h)h(has)g(a)225 2532 y(preset)d(v)n(alue)e(in)h(HANKEL.)g (The)g(presen)o(t)h(v)n(alue)e(for)g Fj(x)1147 2538 y Fh(as)1198 2532 y Fl(is)h Fj(x)1266 2538 y Fh(as)1316 2532 y Fl(=)f(100.)22 b(No)16 b(theoretical)p eop %%Page: 4 4 4 3 bop 225 125 a Ff(4)71 b Fe(\001)78 b Ff(T.)13 b(Wieder)225 233 y Fl(justi\014cation)h(can)g(b)q(e)h(o\013ered)g(for)f(that)g(v)n (alue,)f(but)i(our)f(p)q(ositiv)o(e)g(exp)q(erience)i(with)e(sev)o (eral)225 283 y(test)j(functions.)22 b(One)16 b(w)o(ould)f(exp)q(ect)i (that)f(shifting)e Fj(x)1118 289 y Fh(l)1146 283 y Fl(to)h(larger)h(v)n (alues)f(should)g(increase)225 332 y(the)f(accuracy)h(of)f(appro)o (ximation)d(\(1\))246 470 y Fo(H)281 476 y Fh(\027)301 470 y Fl(\()p Fj(\030)r(;)c(f)t Fl(\()p Fj(x)p Fl(\)\))13 b Fo(\031)508 413 y Fc(Z)550 424 y Fh(x)569 428 y Fb(l)531 508 y Fi(0)583 470 y Fl(\()p Fj(x\030)r Fl(\))659 453 y Fi(1)p Fh(=)p Fi(2)711 470 y Fj(f)t Fl(\()p Fj(x)p Fl(\))p Fj(J)814 476 y Fh(\027)836 470 y Fl(\()p Fj(x\030)r Fl(\))7 b Fj(dx)r Fl(+)1001 413 y Fc(Z)1042 424 y Fa(1)1023 508 y Fh(x)1042 512 y Fb(l)1077 470 y Fl(\()1100 442 y(2)p 1098 460 26 2 v 1098 498 a Fj(\031)1128 470 y Fl(\))1144 453 y Fi(1)p Fh(=)p Fi(2)1196 470 y Fj(f)t Fl(\()p Fj(x)p Fl(\))g(cos)r(\()p Fj(x\030)t Fo(\000)1441 442 y Fj(\027)s(\031)p 1441 460 49 2 v 1455 498 a Fl(2)1496 470 y Fo(\000)1536 442 y Fj(\031)p 1536 460 26 2 v 1538 498 a Fl(4)1566 470 y(\))g Fj(dx:)19 b Fl(\(1\))267 570 y(F)m(or)f(most)g(of)g(the)i (examined)e(test)i(functions)f(this)g(exp)q(ectation)h(is)e (ful\014lled,)h(but)g(test)225 620 y(function)g(2)f(with)h Fj(f)t Fl(\()p Fj(x)p Fl(\))i(=)f Fj(x)709 605 y Fa(\000)p Fi(0)p Fh(:)p Fi(5)787 620 y Fl(exp\()p Fo(\000)p Fj(x)p Fl(\))f(is)g(a)g(coun)o(terexample.)33 b(While)18 b(INTHP)h(w)o(as)225 670 y(successful)f(for)e Fj(x)506 676 y Fh(as)556 670 y Fl(=)g(100)p Fj(:)p Fl(0,)e(INTHP)j(failed)e(\(under)i(our)f(usage\)) h(completely)e(for)g Fj(x)1636 676 y Fh(as)1687 670 y Fl(=)225 720 y(1000)p Fj(:)p Fl(0)c(\(ga)o(v)o(e)i(zero)h(v)n(alues)f (for)g(all)e(transforms\).)18 b(Th)o(us,)13 b(if)f Fj(f)t Fl(\()p Fj(x)p Fl(\))g(=)g(0)h(for)g Fj(x)e(>)h(x)1520 726 y Fh(d)1539 720 y Fl(,)g(then)i(one)225 769 y(should)f(set)i Fj(x)445 775 y Fh(l)469 769 y Fl(=)c Fj(x)536 775 y Fh(d)555 769 y Fl(.)18 b Fj(x)609 775 y Fh(as)658 769 y Fl(ma)o(y)12 b(b)q(e)i(c)o(hanged)g(within)e(the)i(source)h(co)q(de,)f(but)g Fj(x)1497 775 y Fh(l)1522 769 y Fl(will)e(alw)o(a)o(ys)225 819 y(o)o(v)o(erride)i Fj(x)407 825 y Fh(as)443 819 y Fl(.)p eop %%Trailer end userdict /end-hook known{end-hook}if %%EOF SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Fortran77' then mkdir 'Fortran77' fi cd 'Fortran77' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'data' then echo shar: will not over-write existing file "'data'" else cat << SHAR_EOF > 'data' res.err res 1 1.0 5.0 2 SHAR_EOF fi # end of overwriting check if test -f 'driver.f' then echo shar: will not over-write existing file "'driver.f'" else cat << SHAR_EOF > 'driver.f' C BEGIN OF PROGRAM *** DRIVER *** PROGRAM DRIVER C C ################################################################# C C FIRST VERSION: 04.01.1998 C PREVIOUS VERSION: 08.02.1999 C LATEST VERSION: 11.02.1999 C C PROGRAMMED BY: C THOMAS WIEDER C E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE C HTTP://WWW.UNI-KASSEL.DE/~WIEDER C C ACKNOWLEDGMENTS: C C MANY THANKS TO DR. T.R. HOPKINS, THE UNIVERSITY OF KENT, FOR C DETECTING A BUG AND FOR SEVERAL SUGGESTIONS TO IMPROVE THE FORTRAN C CODE (1998). C MANY THANKS TO DR. RICHARD HANSON, ACM-TOMS, FOR SEVERAL C SUGGESTIONS TO IMPROVE THE PROGRAM DESCRIPTION "NUMERICAL C HANKEL TRANSFORM BY THE FORTRAN PROGRAM HANKEL" (1998). C C PURPOSE: C C DRIVER ROUTINE. C EXAMPLE FOR USAGE OF SUBROUTINE *** HANKEL ***. C C MOST IMPORTANT VARIABLES: C C XI = POSITION AT WHICH THE HANKEL TRANSFORM *** HT *** C OF FUNCTION *** FUN *** HAS TO BE EVALUATED. C NU = ORDER OF HANKEL TRANSFORM C = ORDER OF BESSEL FUNCTION *** DHJNUX ***. C C ( infty C HT(FUN,NU,XI) = | (XI * X)**(1/2) JNUX(XI * X) FUN(X) dX C 0 ) C C HT = HANKEL TRANSFORM OF FUNCTION *** FUN ***. C C THE FUNCTION *** FUN *** IS OF THE FORM: Y = FUN(X) C C IMODE = LOGICAL FLAG C IMODE EQUAL TO 1 --> FUNCTION *** FUN *** IS GIVEN ANALYTICALLY C WITHIN FUNCTION *** DHFUNC ***. C IMODE EQUAL TO 2 --> FUNCTION *** FUN *** IS GIVEN AS TABULATED C DATA: C X_1,Y_1, ..., X_N,Y_N, ..., X_NEND,Y_NEND. C THESE DATA ARE PASSED TO *** HANKEL *** C VIA AN INPUT FILE AND A CALL TO SUBROUTINE C *** FUNK2 *** C PRIOR TO THE CALL OF *** HANKEL ***. C IMODE EQUAL TO 3 --> FUNCTION *** FUN *** IS GIVEN AS TABULATED C DATA C X_1,Y_1, ..., X_N,Y_N, ..., X_NEND,Y_NEND. C THESE DATA ARE PASSED TO *** HANKEL *** C VIA A CALL TO SUBROUTINE *** FUNK3 *** C PRIOR TO THE CALL OF *** HANKEL ***. C C IVERBO = 0 --> SUBROUTINE *** HANKEL *** RUNS IN QUIET MODE. C IVERBO = 1 --> *** SUBROUTINE HANKEL *** SENDS SOME REPORT C TO UNIT * (STANDARD OUTPUT). C IVERBO > 1 --> *** SUBROUTINE HANKEL *** SENDS SOME REPORT C TO UNIT * AND ADDITIONAL REPORT TO C FILE *** FILEER ***. C USUALLY YOU WILL CALL *** HANKEL *** WITH *** IVERBO *** = 0. C *** IVERBO *** 1 IS FOR DEBUGGING PURPOSES ONLY! C C ################################################################# C C IMPLICIT NONE C C ----------------------------------------------------------------- C C .. Scalars in Common .. CHARACTER*80 FILEER,FILEIN,FILEOU C .. C .. Local Scalars .. DOUBLE PRECISION DUMMY,HT,JNUXV,NU,SLTN,XI INTEGER ICON,IERR,IERRJ,IMODE,IVERBO,LUERR,LUIN,LUOUT,N,NEND C .. C .. Local Arrays .. DOUBLE PRECISION XDAT(1000),YDAT(1000) C .. C .. External Subroutines .. EXTERNAL DHINIT,DHJNUX,ERRMSS,FUNK1,FUNK2,FUNK3,HANKEL C .. C .. Intrinsic Functions .. INTRINSIC DBLE,FLOAT C .. C .. Common blocks .. COMMON /NAMEN/FILEIN,FILEOU,FILEER integer i1mach, nout, nin nin =i1mach(1) nout = i1mach(2) C .. 10 WRITE(NOUT,FMT=9000) WRITE(NOUT,FMT=9010) C C ------------------------------------------------------------------- C C SET SOME PARAMETERS: C C *** IVERBO *** CONTROLS THE VERBOSITY OF OUTPUT IVERBO = 1 C C *** LUIN *** IS THE LOGICAL NUMBER OF THE INPUT FILE LUIN = 66 C C *** LUOUT *** IS THE LOGICAL NUBER OF THE OUTPUT FILE LUOUT = 67 C C *** LUERR *** IS THE LOGICAL NUMBER OF THE ERROR LOG FILE LUERR = 68 C C SET *** DUMMY *** TO AN ARBITRARY REAL NUMBER DUMMY = 0.0D0 C C ------------------------------------------------------------------ C C OPEN FILES C WRITE(NOUT,FMT=9020) READ(NIN,FMT=*,ERR=20,END=20) FILEER GO TO 30 20 FILEER = 'hankel.err' 30 WRITE(NOUT,FMT=9030) FILEER WRITE(NOUT,FMT=9040) READ(NIN,FMT=9050,ERR=40,END=40) FILEOU GO TO 50 40 FILEOU = 'hankel.out' 50 WRITE(NOUT,FMT=9060) FILEOU C FILE *** FILEER *** STORES MESSAGES FROM *** PROGRAM HANKEL ***. C SEARCH FOR ERROR MESSAGES IN FILE *** FILEER *** ! OPEN (UNIT=LUERR,FILE=FILEER,STATUS='UNKNOWN',ACCESS='SEQUENTIAL', + ERR=90) C FILE *** FILEOU *** STORES RESULTS FROM *** PROGRAM HANKEL ***. OPEN (UNIT=LUOUT,FILE=FILEOU,STATUS='UNKNOWN',ACCESS='SEQUENTIAL', + ERR=90) C C ------------------------------------------------------------------- C C INITIALIZE: C C INITIALIZE THE LOGICAL NUMBERS FOR INPUT AND OUTPUT FILES C CALL DHINIT(LUIN,LUOUT,LUERR) C C ------------------------------------------------------------------- C C READ SOME INPUT DATA: C C TELL WHERE THE FUNCTION *** FUN *** IS: 60 WRITE(NOUT,FMT=9070) READ(NIN,FMT=*) IMODE IF (IMODE.EQ.0) STOP C C *** NU *** IS THE ORDER OF THE HANKEL TRANSFORM: WRITE(NOUT,FMT=9080) READ(NIN,FMT=*) NU C C SET PARAMETER *** XI ***: WRITE(NOUT,FMT=9090) READ(NIN,FMT=*) XI C C ------------------------------------------------------------------ C C DO THE HANKEL TRANSFORM: C C ------------------------------------------------------------------ C IF (IMODE.EQ.1) THEN C C FUNCTION *** FUN *** IS CODED WITHIN SUBROUTINE *** DHFUNC ***, C THEREFORE, FIRST CALL SUBROUTINE *** FUNK1 *** !!!!!!!!!!!!!!!!! C CALL FUNK1(DUMMY,DUMMY,IERR) C C DO NOT CHANGE *** IMODE *** AFTER CALL OF *** FUNK1 *** ! C THEN CALL SUBROUTINE *** HANKEL ***, C TO DO THE HANKEL TRANSFORM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C CALL HANKEL(XI,NU,IVERBO,IERR,HT) C IF (IERR.NE.0) THEN WRITE(NOUT,FMT=9100) END IF C C ------------------------------------------------------------------ C ELSE IF (IMODE.EQ.2) THEN C C DATA ARE PROVIDED IN AN INPUT FILE. C 70 WRITE(NOUT,FMT=9110) READ(NIN,FMT=9120,ERR=70) FILEIN C C FIRST CALL SUBROUTINE *** FUNK2 *** C TO READ THE DATA FROM FILE!!!!!!!!!!!!!!!!! C CALL FUNK2(DUMMY,DUMMY,IERR) C C DO NOT CHANGE *** IMODE *** AFTER CALL OF *** FUNK2 *** ! C THEN CALL SUBROUTINE *** HANKEL ***, C TO DO THE HANKEL TRANSFORM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C CALL HANKEL(XI,NU,IVERBO,IERR,HT) C IF (IERR.NE.0) THEN WRITE(NOUT,FMT=9100) END IF C C ------------------------------------------------------------------ C ELSE IF (IMODE.EQ.3) THEN C C DATA ARE TRANSFERED VIA CALL. C C GENERATE DATA. C NEND = 200 C C AS AN EXAMPLE: C FUN(X)=X**(NU+0.5D0) FOR X < 1 C AND C FUN(X)=0 FOR X=> 1 DO 80 N = 1,NEND XDAT(N) = DBLE(FLOAT(N)/FLOAT(NEND)) C JUST AN EXAMPLE: YDAT(N) = XDAT(N)** (NU+0.5D0) 80 CONTINUE C C DATA ARE PASSED VIA CALL. C THEREFORE, FIRST CALL SUBROUTINE *** FUNK3 *** !!!!!!!!!!!!!!!!! C CALL FUNK3(DUMMY,DUMMY,XDAT,YDAT,NEND,IERR) C C DO NOT CHANGE *** IMODE *** AFTER CALL OF *** FUNK3 *** ! C THEN CALL SUBROUTINE *** HANKEL ***, C TO DO THE HANKEL TRANSFORM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C CALL HANKEL(XI,NU,IVERBO,IERR,HT) C IF (IERR.NE.0) THEN WRITE(NOUT,FMT=9100) END IF C ELSE C WRONG INPUT FOR *** IMODE ***, READ AGAIN: GO TO 60 END IF C C COMPARE ANALYTICAL AND NUMERICAL SOLUTION: C IF (IERR.EQ.0 .AND. IVERBO.NE.0) THEN C CALCULATE THE ANALYTICAL SOLUTION FOR COMPARISION: IERRJ = 0 C CALL DHJNUX(XI,NU+1.0D0,IERRJ,JNUXV) IF (IERRJ.NE.0) THEN WRITE(NOUT,FMT=9130) END IF C IF (IVERBO.GT.1) THEN SLTN = XI** (-0.5D0)*JNUXV WRITE(NOUT,FMT=9140) NU,XI,SLTN,HT END IF C END IF C C ------------------------------------------------------------------- C WRITE(NOUT,FMT=9150) READ(NIN,FMT=*) ICON IF (ICON.EQ.1) GO TO 10 C WRITE(NOUT,FMT=9160) WRITE(NOUT,FMT=9170) C C ------------------------------------------------------------------- C C ERROR HANDLING: C GO TO 100 90 WRITE(NOUT,FMT=9180) FILEER IERR = -10 CALL ERRMSS(IERR) C 100 STOP 9000 FORMAT (/,' WELCOME TO H A N K E L !',/, + ' YOU ARE COMMUNICATING WITH D R I V E R,',/, + ' THE DRIVER ROUTINE TO H A N K E L.') 9010 FORMAT (/,' AS AN EXAMPLE, WE WILL NUMERICALLY HANKEL TRANSFORM', + /,' FUN(X) = X**(NU+1/2) FOR X < 1 AND = 0 ELSEWHERE',/, + ' AT POSITION XI FOR TRANSFORMATION ORDER NU.',/, + ' THE ANALYTICAL SOLUTION IS',/, + ' HT(FUN(X),NU,XI) = XI**(-1/2)*J(NU+1,XI)',/, + ' (J(NU,X) IS THE BESSEL FUNCTION OF THE FIRST KIND.)') 9020 FORMAT (/,' MESSAGE DRIVER: GIVE NAME FOR PROTOCOL FILE',/, + ' (E.G. hankel.err):') 9030 FORMAT (' MESSAGE DRIVER: ERROR LOG FILE =',A80) 9040 FORMAT (/,' MESSAGE DRIVER: GIVE NAME FOR OUTPUT FILE',/, + ' (E.G. hankel.out):') 9050 FORMAT (A80) 9060 FORMAT (' MESSAGE DRIVER: RESULT FILE =',A80) 9070 FORMAT (/,' FUNCTION *** FUN(X) *** IS EITHER CODED WITHIN',/, + ' FUNCTION *** DHFUNC *** = 1',/,' OR',/, + ' TABULATED IN A FILE = 2 (VERY SLOW!)',/,' OR',/, + ' PASSED AS TABLE VIA CALL TO HANKEL = 3 (SLOW!)',/,' OR ', + /,' STOP = 0',/,' GIVE CHOICE (1 OR 2 OR 3 OR 0):',/) 9080 FORMAT (/,' GIVE THE ORDER *** NU ***:',/,' (E.G. 1.0):',/) 9090 FORMAT (/,' GIVE THE ARGUMENT *** XI ***:',/,' (E.G. 5.0):',/) 9100 FORMAT (1X,' SERIOUS ERROR! HANKEL TRANSFORMATION FAILED!') 9110 FORMAT (/,' MESSAGE DRIVER: GIVE NAME OF INPUT FILE:',/) 9120 FORMAT (A80) 9130 FORMAT (/, +' AN ERROR DURING FUNCTION EVALUATION OCCURED!' + ,/,' THE HANKEL TRANSFORM HAS FAILED!') 9140 FORMAT (/, + ' RESULTS FOR TEST FUNCTION FUN(X)=X**(NU+0.5D0) FOR X = 1.0:' + ,/, + ' (IT IS ASSUMED THAT *** IFUNC=0 *** IN *** DHFUNC *** )', + /, + ' NU XI ANALYTICAL SOLUTION NUMERICAL SOLUTION:' + ,/,D17.9,3X,D17.9,3X,D17.9,3X,D17.9) 9150 FORMAT (/,' DO YOU WANT ANOTHER CALCULATION = 1',/, + ' STOP = 2',/,' GIVE 1 OR 2:',/) 9160 FORMAT (/,' A PROTOCOL HAS BEEN WRITTEN INTO A FILE,',/, + ' RESULTS HAVE BEEN WRITTEN INTO A FILE.') 9170 FORMAT (/,' GOODBYE!') 9180 FORMAT (/,' MESSAGE DHINIT: ERROR ON OPENING FILE:',/,A80) END C END OF PROGRAM *** DRIVER *** C C BEGIN OF SUBROUTINE *** DHFUNC *** SUBROUTINE DHFUNC(NU,X,Y,XLIMIT,IERR) C C FIRST VERSION: 04.01.1998 C PREVIOUS VERSION: 05.20.1999 C LATEST VERSION: 06.02.1999 C C ################################################################# C C PROGRAMMED BY: C THOMAS WIEDER C E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE C HTTP://WWW.UNI-KASSEL.DE/~WIEDER C C PURPOSE: C SUBROUTINE *** DHFUNC *** PROVIDES THE FUNCTION *** FUN *** C WHICH SHOULD BE HANKEL-INVERTED BY SUBROUTINE *** HANKEL ***. C C THE FORM OF SUBROUTINE *** DHFUNC *** IS: C C SUBROUTINE DHFUNC(NU,X,Y,XLIMIT,LUERR,IERR) C IMPLICIT NONE C INTEGER IERR,LUIN,LUOUT,LUERR C DOUBLE PRECISION X,Y,NU,XASYMP,XLIMIT C COMMON /INOUTE/ LUIN,LUOUT,LUERR C C ############################################################## C C INPUT VARIABLES: NU, X,LUERR C OUPUT VARIABLES: Y,XLIMIT C C ############################################################## C C INITIALIZE THE ERROR FLAG. C *** IERR=0 *** MEANS THAT NO ERRORS OCCURRED SO FAR. C (OF COURSE YOU OMIT THE 'C'): C C IERR=0 C C THE INFINITE INTEGRAL OVER THE FUNCTION *** FUN *** C WILL BE SPLIT INTO A FINITE INTEGRAL FROM 0 TO C *** XLIMIT *** AND AN INFINITE INTEGRAL FROM C *** XLIMIT *** TO INFINITY. C *** XLIMIT *** IS THE UPPER LIMIT OF THE FINITE INTEGRAL. C YOU NEED NOT TO SPECIFY *** XLIMIT ***. IN THAT CASE YOU C JUST GIVE A NEGATIVE VALUE, E.G. *** XLIMIT=-1.0D0 *** C (OF COURSE YOU OMIT THE 'C'): C C XLIMIT=-1.0D0 C C HERE YOU GIVE YOUR FUNCTION *** FUN *** C Y=FUN(X,NU) C (OF COURSE YOU OMIT THE 'C'): C C Y=... C C IF AN ERROR HAS OCCURRED, THEN PUT *** IERR=1 *** ! C C RETURN C END C C CALLING SEQUENCE: C CALLED BY FUNCTION *** FUNINT ***. C C MOST IMPORTANT VARIABLES: C C X = INPUT VALUE = INDEPENDENT VARIABLE OF *** FUN ***. C Y = OUTPUT VALUE = DEPENDENT VARIABLE OF *** FUN ***. C Y = FUN(X) C XLIMIT = UPPER LIMIT OF THE FINITE INTEGRAL. C LUERR = LOGICAL NUMBER OF ERROR LOG FILE C IERR = ERROR FLAG C C IERR = 0 --> NO ERROR OCCURED DURING EVALUATION OF C SUBROUTINE *** DHFUNC ***. C C ################################################################# C C IMPLICIT NONE C C ------------------------------------------------------------------- C C THE FOLLOWING VARIBALES ARE NECCESSARY C FOR DEBUGGING PURPOSES ONLY. C THE COMMON BLOCK *** SOLUTI *** IS NECEESSARY C FOR DEBUGGING PURPOSES ONLY. C C ------------------------------------------------------------------ C C SET SOME PARAMETERS: C C .. Scalar Arguments .. DOUBLE PRECISION NU,X,XLIMIT,Y INTEGER IERR C .. C .. Scalars in Common .. DOUBLE PRECISION SLTN,XISLTN INTEGER IFUNC,LUERR,LUIN,LUOUT C .. C .. Local Scalars .. DOUBLE PRECISION JNUXV,PI,XASYMP INTEGER IERRJ C .. C .. External Subroutines .. EXTERNAL DHJNUX C .. C .. Intrinsic Functions .. INTRINSIC ASIN,DCOS,DSIN,EXP integer i1mach, nout C .. C .. Common blocks .. COMMON /INOUTE/LUIN,LUOUT,LUERR COMMON /SOLUTI/XISLTN,SLTN,IFUNC C .. nout = i1mach(2) IERR = 0 XLIMIT = -1.0D0 C Y = 0.0D0 C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C *** 24.02.1998 *** A GOOD VALUE FOR *** XASYMP *** = 100.0 C USUALLY, YOU DO NOT CHANGE THE FOLLOWING LINE: XASYMP = 100.0D0 C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C !!! IFUNC MUST BE EQUAL TO ZERO IN ORDER TO HAVE YOUR FUNCTION C FUN HANKEL TRANSFORMED !!! C FOR DEBUGGING *** HANKEL ** ONE MAY CHOOSE ONE OF THE TEST C FUNCTIONS GIVEN BELOW BY SETTING *** IFUNC *** EQUAL TO VALUES C FROM 1 TO 4. C USUALLY, YOU DO NOT CHANGE THE FOLLOWING LINE: IFUNC = 0 C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C ----------------------------------------------------------------- C C GIVE THE FUNCTION *** FUN ***: C IF (IFUNC.EQ.0) THEN C FOR *** IFUNC=0 *** YOUR FUNCTION WILL BE CALCULATED! C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C HERE YOU HAVE TO TYPE YOUR FUNCTION: C HERE YOU HAVE TO TYPE YOUR FUNCTION: C HERE YOU HAVE TO TYPE YOUR FUNCTION: C C WE PUT X**(NU+0.5D0) FOR X < 1 AS AN EXAMPLE: C C WE SET *** XLIMIT *** = 1 BECAUSE OF THE DEFINITION RANGE. XLIMIT = 1.0D0 C C THIS IS THE FUNCTION *** FUN(X) ***: IF (X.LE.1.0D0) THEN C BE V E R Y CAUTIOUS ABOUT 0**0! IF (X.NE.0.0D0) Y = X** (NU+0.5D0) ELSE IF (X.GT.1.0D0) THEN Y = 0.0D0 END IF C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C ------------------------------------------------------------------ C C HERE SOME TEST FUNCTIONS: C ELSE IF (IFUNC.EQ.1) THEN C AN ALGEBRAIC TEST FUNCTION. C SEE: ERDELYI, W. MAGNUS, F. OBERHETTINGER, F.G. TRICOMI: 8.5. (6) C THIS TEST FUNCTION IS EQUAL TO ZERO ABOVE *** X = 1.0 ***. C THEREFORE WE SET *** XLIMIT = 1.0D0 ***: XLIMIT = 1.0D0 IF (NU.LE.-1.0D0) THEN WRITE (LUERR,FMT=9000) IFUNC,NU WRITE(NOUT,FMT=9000) IFUNC,NU STOP END IF IF (NU.LT.0.0D0 .AND. X.EQ.0.0D0) THEN Y = 0.0D0 RETURN END IF IF (X.LE.1.0D0) THEN C BE V E R Y CAUTIOUS ABOUT 0**0! IF (X.GT.0.0D0) Y = X** (NU+0.5D0) ELSE IF (X.GT.1.0D0) THEN Y = 0.0D0 END IF C C SOLUTION: C XI**(-0.5) * JNUX(XI,NU+1) IERRJ = 0 C CALL DHJNUX(XISLTN,NU+1.0D0,IERRJ,JNUXV) C IF (IERRJ.NE.0) THEN IERR = 1 WRITE(NOUT,FMT=9010) WRITE (LUERR,FMT=9010) END IF SLTN = XISLTN** (-0.5D0)*JNUXV C C ------------------------------------------------------------------ C ELSE IF (IFUNC.EQ.2) THEN C AN EXPONENTIAL TEST FUNCTION. C SEE: ERDELYI, W. MAGNUS, F. OBERHETTINGER, F.G. TRICOMI: 8.6. (1) XLIMIT = 100.0D0 IF (NU.LE.-1.0D0) THEN WRITE (LUERR,FMT=9000) IFUNC,NU WRITE(NOUT,FMT=9000) IFUNC,NU STOP END IF IF (X.EQ.0.0D0) THEN Y = 0.0D0 RETURN END IF C BE V E R Y CAUTIOUS ABOUT 0**0! IF (X.GT.0.0D0) Y = X** (-0.5D0)*EXP(-X) C C SOLUTION: C XI**(0.5-NU) * (1.0+XI**2)**(-0.5) * ((1.0+XI**2)**(-0.5) - 1)**N SLTN = XISLTN** (0.5D0-NU)* (1.0D0+XISLTN**2)** (-0.5D0)* + ((1.0D0+XISLTN**2)** (0.5D0)-1.0D0)**NU C C ------------------------------------------------------------------ C ELSE IF (IFUNC.EQ.3) THEN C AN TRIGONOMETRIC TEST FUNCTION. C SEE: ERDELYI, W. MAGNUS, F. OBERHETTINGER, F.G. TRICOMI: 8.7. (1) PI = 3.141592653D0 IF (XISLTN.EQ.1.0D0) THEN WRITE(NOUT,FMT=9020) IFUNC,XISLTN WRITE (LUERR,FMT=9020) IFUNC,XISLTN STOP END IF IF (NU.LE.-2.0D0) THEN WRITE (LUERR,FMT=9000) IFUNC,NU WRITE(NOUT,FMT=9000) IFUNC,NU STOP END IF IF (X.EQ.0.0D0) THEN Y = 0.0D0 RETURN END IF IF (X.LE.XASYMP) THEN C BE V E R Y CAUTIOUS ABOUT 0**0! IF (X.GT.0.0D0) Y = X** (-0.5D0)*DSIN(X) ELSE IF (X.GT.XASYMP) THEN Y = 0.0D0 END IF C C SOLUTION: C FOR 0 < XI < 1: C COS(0.5 * PI * NU) * XI**(NU+0.5) * (1 - XI**2)**(-0.5) * C (1 + (1 - XI**2)**0.5)**(-NU) C FOR 1 < XI < INFINITY: C (1/NU) * XI**0.5 * SIN(NU * (ASIN(1/XI))) IF (XISLTN.LT.1.0D0) THEN SLTN = DCOS(0.5D0*PI*NU)*XISLTN** (NU+0.5D0)* + (1.0D0-XISLTN**2)** (-0.5D0)* + (1.0D0+ (1.0D0-XISLTN**2)**0.5D0)** (-NU) ELSE IF (XISLTN.GT.1.0D0) THEN SLTN = XISLTN**0.5D0* (XISLTN**2-1.0D0)** (-0.5D0)* + DSIN(NU*ASIN(1.0D0/XISLTN)) END IF C C ------------------------------------------------------------------ C ELSE IF (IFUNC.EQ.4) THEN C A TEST FUNCTION CONTAINING A BESSEL FUNCTION. C SEE: ERDELYI, W. MAGNUS, F. OBERHETTINGER, F.G. TRICOMI: C 8.11. (1) IF (NU.LE.-1.0D0) THEN WRITE (LUERR,FMT=9000) IFUNC,NU WRITE(NOUT,FMT=9000) IFUNC,NU STOP END IF IF (XISLTN.LE.1.0D0) THEN Y = 0.0D0 SLTN = 0.0D0 RETURN END IF IF (X.EQ.0.0D0) THEN Y = 0.0D0 RETURN END IF C C *** AMENDMENT SUGGESTED BY DR. T.R. HOPKINS, 15.04.1998 *** IERRJ = 0 C IF (X.LE.XASYMP) THEN C CALL DHJNUX(X,NU-1.0D0,IERRJ,JNUXV) C C BE V E R Y CAUTIOUS ABOUT 0**0! IF (X.GT.0.0D0) Y = X** (-0.5D0)*JNUXV C ELSE IF (X.GT.XASYMP) THEN Y = 0.0D0 END IF C IF (IERRJ.NE.0) THEN WRITE(NOUT,FMT=9030) IERRJ WRITE(NOUT,FMT=9040) X END IF C C SOLUTION: C ASSUME A=1 IN FORMULA 8.11.(1)! C FOR 0 < XI < 1: C 0.0 C FOR 1 < XI < INFTY: C Y**(-NU+0.5) IF (XISLTN.LT.1.0D0) THEN SLTN = 0.0D0 ELSE IF (XISLTN.GT.1.0D0) THEN SLTN = XISLTN** (-NU+0.5D0) END IF C C END OF IF-CLAUSE *** IFUNC *** END IF C RETURN 9000 FORMAT (/,' MESSAGE FUNC:',/,' TEST FUNCTION',I6, + ' NOT DEFINED FOR *** NU *** =',D17.9) 9010 FORMAT (/,' MESSAGE FUNC:',/, + ' AN ERROR DURING FUNCTION EVALUATIO OCCURED!',/, + ' THE HANKEL TRANSFORM HAS FAILED!') 9020 FORMAT (/,' MESSAGE FUNC:',/,' HANKEL TRANSFORM OF TEST FUNCTION ' + ,I6,/,' NOT DEFINED FOR *** XI *** =',D17.9) 9030 FORMAT (/,' AN ERROR OCCURED DURING CALL OF *** DHJNUX ***',/, + ' *** IERRJ *** =',I6) 9040 FORMAT (' ARGUMENT *** X *** =',D17.9) END C END OF SUBROUTINE *** DHFUNC *** C C BEGIN OF SUBROUTINE *** FUNK1 *** SUBROUTINE FUNK1(X,Y,IERR4) C C FIRST VERSION: 04.01.1998 C PREVIOUS VERSION: 06.10.1998 C LATEST VERSION: 07.10.1998 C C ################################################################## C C PROGRAMMED BY: C THOMAS WIEDER C E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE C C PURPOSE: C SUBROUTINE *** FUNK *** PROVIDES THE FUNCTION *** FUN *** C WHICH SHOULD BE HANKEL-INVERTED BY SUBROUTINE *** HANKEL ***. C C *** FUN *** HAS TO BE PROVIDED ANALYTICALLY BY C SUBROUTINE *** DHFUNC *** !!!!!!!!!!!!!!!!!!!! C C CALLING SEQUENCE: C CALLED BY SUBROUTINE *** HANKEL ***. C CALLS SUBROUTINE *** DHFUNC ***. C C MOST IMPORTANT VARIABLES: C C X = INDEPENDENT REAL VARIABLE OF *** FUN *** (INPUT). C Y = DEPENDENT REAL VARIABLE OF *** FUN *** (OUTPUT). C Y = FUN(X) C C IDAT (OUTPUT) C IDAT EQUAL TO 1 --> FUNCTION *** FUN *** IS GIVEN ANALYTICALLY C IDAT EQUAL TO 2 --> FUNCTION *** FUN *** IS GIVEN AS TABULATED C DATA IN A FILE C IDAT EQUAL TO 3 --> FUNCTION *** FUN *** IS GIVEN AS TABULATED C DATA PASSED BY CALL TO *** HANKEL *** C C IERR4 = ERROR FLAG (OUTPUT) C IERR4 EQUAL TO 0 ON RETURN --> NO ERROR OCCURED. C IERR4 EQUAL TO 1 ON RETURN --> AN ERROR OCCURED DURING C EVALUATION OF FUNCTION *** FUN *** C WITHIN SUBROUTINE *** DHFUNC ***. C C ################################################################## C C IMPLICIT NONE C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C FUNCTION *** FUN *** IS GIVEN ANALYTICALLY WITHIN C SUBROUTINE *** DHFUNC ***. C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C .. Parameters .. INTEGER ITEST PARAMETER (ITEST=0) C .. C .. Scalar Arguments .. DOUBLE PRECISION X,Y INTEGER IERR4 C .. C .. Scalars in Common .. DOUBLE PRECISION NU,XI,XLIMIT INTEGER IDAT,IREAD,LUERR,LUIN,LUOUT C .. C .. Local Scalars .. DOUBLE PRECISION RESULT INTEGER IERR5,N,NEND C .. C .. External Subroutines .. EXTERNAL DHFUNC C .. C .. Common blocks .. COMMON /ARGUME/XI,NU,XLIMIT COMMON /INOUTE/LUIN,LUOUT,LUERR COMMON /LOGIK/IDAT,IREAD integer i1mach, nout nout = i1mach(2) C .. IF (IREAD.NE.1) THEN C C SAVE INPUT C IREAD = 1 END IF C IERR4 = 0 IERR5 = 0 C CALL DHFUNC(NU,X,RESULT,XLIMIT,IERR5) C IF (IERR5.NE.0) THEN WRITE (LUERR,FMT=9000) IERR5 WRITE(NOUT,FMT=9000) IERR5 WRITE (LUERR,FMT=9010) X IERR4 = 1 RETURN END IF C IF (ITEST.EQ.1) THEN NEND = 100 WRITE(NOUT,FMT=9020) X = 0.1D0 DO 10 N = 1,NEND WRITE(NOUT,FMT=*) N,X,RESULT X = X + 0.1D0 10 CONTINUE END IF C C SUCCESSFUL RETURN FROM *** FUNK1 *** IDAT = 1 Y = RESULT C RETURN 9000 FORMAT (/,' MESSAGE FUNK1:',/, + ' ERROR ON CALCULATING FUNCTION VA LUE!',/, + ' ERROR IN SUBROUTINE *** DHFUNC *** OR SIMILAR SUBROUTINE!' + ,/,' IERR5 = ',I6) 9010 FORMAT (1X,' FUNCTION ARGUMENT X =',D17.9) 9020 FORMAT (/, + ' MESSAGE FUNK1: TEST RUN FOR SUBROUTINE *** DHFUNC ** *:' + ) END C END OF SUBROUTINE *** FUNK1 *** C C BEGIN OF SUBROUTINE *** FUNK2 *** SUBROUTINE FUNK2(X,Y,IERR4) C C FIRST VERSION: 04.01.1998 C PREVIOUS VERSION: 10.10.1998 C LATEST VERSION: 13.10.1998 C C ################################################################## C C PROGRAMMED BY: C THOMAS WIEDER C E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE C C PURPOSE: C SUBROUTINE *** FUNK *** PROVIDES THE FUNCTION *** FUN *** C WHICH SHOULD BE HANKEL-INVERTED BY SUBROUTINE *** HANKEL ***. C C *** FUN *** HAS TO BE PROVIDED AS TABULATED DATA X,Y C IN A FILE. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C CALLING SEQUENCE: C CALLED BY SUBROUTINE *** HANKEL ***. C CALLS SUBROUTINE *** DIVDIF ***. C C MOST IMPORTANT VARIABLES: C C X = INDEPENDENT REAL VARIABLE OF *** FUN *** (INPUT). C Y = DEPENDENT REAL VARIABLE OF *** FUN *** (OUTPUT). C Y = FUN(X) C C NEND = NUMBER OF DATA PAIRS (INPUT). C NMAX = MAXIMUM NUMBER OF ALLOWED FUNCTION VALUES (OUTPUT). C C IDAT (OUTPUT) C IDAT EQUAL TO 1 --> FUNCTION *** FUN *** IS GIVEN ANALYTICALLY C IDAT EQUAL TO 2 --> FUNCTION *** FUN *** IS GIVEN AS TABULATED C DATA IN A FILE C IDAT EQUAL TO 3 --> FUNCTION *** FUN *** IS GIVEN AS TABULATED C DATA PASSED BY CALL TO *** HANKEL *** C C IERR4 = ERROR FLAG (OUTPUT) C IERR4 EQUAL TO 0 ON RETURN --> NO ERROR OCCURED. C IERR4 EQUAL TO 1 ON RETURN --> AN ERROR OCCURED DURING C EVALUATION OF FUNCTION *** FUN *** C C ################################################################## C C IMPLICIT NONE C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C FUNCTION *** FUN *** IS GIVEN AS DATA TABLE IN AN INPUT FILE. C THE FORMAT OF THE DATA IS: C X_1 Y_1 C ... ... C X_N Y_N C ... ... C X_NEND Y_NEND C I.E. EACH LINE CONTAINS ONE DATA POINT (X,Y). C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C .. Parameters .. INTEGER NMAX PARAMETER (NMAX=1000) INTEGER ITEST PARAMETER (ITEST=0) C .. C .. Scalar Arguments .. DOUBLE PRECISION X,Y INTEGER IERR4 C .. C .. Scalars in Common .. DOUBLE PRECISION NU,XDATMA,XDATMI,XI,XLIMIT,YDATMA,YDATMI INTEGER IDAT,IREAD,LUERR,LUIN,LUOUT,NEND C .. C .. Arrays in Common .. DOUBLE PRECISION XDAT(NMAX),YDAT(NMAX) C .. C .. Local Scalars .. DOUBLE PRECISION RESULT INTEGER IDEG,IERR6,IERRRW,N C .. C .. Local Arrays .. DOUBLE PRECISION TABLE(NMAX,NMAX) C .. C .. External Subroutines .. EXTERNAL DHRW,DIVDIF C .. C .. Common blocks .. COMMON /ARGUME/XI,NU,XLIMIT COMMON /DATEN/XDAT,YDAT,XDATMI,XDATMA,YDATMI,YDATMA,NEND COMMON /INOUTE/LUIN,LUOUT,LUERR COMMON /LOGIK/IDAT,IREAD integer i1mach, nout nout = i1mach(2) C .. IERR4 = 0 C C READ DATA FROM FILE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C IF (IREAD.NE.1) THEN C C *** LUIN *** = LOGICAL NUMBER OF INPUT FILE C CALL DHRW(IERRRW) C READ DATA ONLY ONCE: IREAD = 1 C IF (IERRRW.NE.0) THEN WRITE(NOUT,FMT=9000) IERRRW WRITE (LUERR,FMT=9000) IERRRW C AN ERROR OCCURED, DO NOT CONTINUE, GO TO STOP! WRITE (LUERR,FMT=9010) DO 10 N = 1,NEND WRITE (LUERR,FMT=*) N,XDAT(N),YDAT(N) 10 CONTINUE STOP END IF C IF (ITEST.EQ.1) THEN WRITE(NOUT,FMT=9020) DO 20 N = 1,NEND WRITE(NOUT,FMT=*) N,XDAT(N),YDAT(N) 20 CONTINUE END IF C C SUCCESSFUL RETURN FROM *** FUNK2 ***, DATA HAVE BEEN READ! IDAT = 2 RETURN C C END OF IF-CLAUSE *** IREAD *** END IF C C INTERPOLATE FUNCTION VALUE FOR ARGUMENT *** X *** !!!!!!!!!!!!!!!! C C DO NOT PERFORM THE DEFINITE INTEGRAL ABOVE THE TABULATED C DATA RANGE: XLIMIT = XDATMA C IF (X.GE.XDATMI .AND. X.LE.XDATMA) THEN C C *** X*** LIES WITHIN TABULATED DATA RANGE. TRY THE INTERPOLATION! C C CHOOSE *** IDEG *** NOT TOO HIGH!!! IDEG = 5 IERR6 = 0 C CALL DIVDIF(XDAT,YDAT,NEND,IDEG,IDEG,TABLE,X,RESULT,IERR6) C CALL DIVDIF (XDAT,YDAT,NEND,IDEG,IDEG,X,RESULT,IERR6) C IF (IERR6.NE.0) THEN WRITE (LUERR,FMT=9030) IERR6 WRITE(NOUT,FMT=9030) IERR6 WRITE (LUERR,FMT=9040) X WRITE (LUERR,FMT=9050) NEND IERR4 = 2 WRITE (LUERR,FMT=9060) IERR4 RETURN END IF C C SUCCESSFUL RETURN FROM *** DIVDIF *** Y = RESULT C ELSE C C *** X *** LIES OUTSIDE THE TABULATED DATA RANGE! C NO INTERPOLATION POSSIBLE! C THIS CASE WILL LEGALLY HAPPEN, SINCE THE HANKEL TRANSFORM C INCLUDES AN INTEGRATION TO INFINITY. C Y = 0.0D0 C C END IF-CLAUSE *** X.LE.XDATMA ***. END IF C RETURN 9000 FORMAT (/,' MESSAGE FUNK2:',/,' ERROR ON READING DATA! IERRRW =', + I6) 9010 FORMAT (/,' MESSAGE FUNK2: DATA READ FROM INPUT FILE:') 9020 FORMAT (' MESSAGE FUNK2: DATA READ FROM FILE:') 9030 FORMAT (/,' MESSAGE FUNK2:',/, + ' ERROR ON INTERPOLATING FUNCTION VALUE!',/, + ' ERROR IN SUBROUTINE *** DIVDIF ***!',/,' IERR6 = ',I6) 9040 FORMAT (1X,' ARGUMENT *** X *** =',D17.9) 9050 FORMAT (1X,' NUMBER OF DATA PAIRS *** NEND *** =',I6) 9060 FORMAT (/,' MESSAGE FUNK2: IERR4 =',I6) END C END OF SUBROUTINE FUNK2 *** C C BEGIN OF SUBROUTINE *** FUNK3 *** SUBROUTINE FUNK3(X,Y,XDATIN,YDATIN,NENDIN,IERR4) C C FIRST VERSION: 04.01.1998 C PREVIOUS VERSION: 10.10.1998 C LATEST VERSION: 13.10.1998 C C ################################################################## C C PROGRAMMED BY: C THOMAS WIEDER C E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE C C PURPOSE: C SUBROUTINE *** FUNK *** PROVIDES THE FUNCTION *** FUN *** C WHICH SHOULD BE HANKEL-INVERTED BY SUBROUTINE *** HANKEL ***. C C *** FUN *** HAS TO BE PROVIDED AS FIELDS X(N),Y(N) PASSED TO C *** FUNK3 *** BY THE CALLING PROGRAM. !!!!!!!!!!!!!!!!!!!! C C CALLING SEQUENCE: C CALLED BY SUBROUTINE *** HANKEL ***. C CALLS SUBROUTINE *** DIVDIF ***. C C MOST IMPORTANT VARIABLES: C C X = INDEPENDENT REAL VARIABLE OF *** FUN *** (INPUT). C Y = DEPENDENT REAL VARIABLE OF *** FUN *** (OUTPUT). C Y = FUN(X) C C NEND = NUMBER OF DATA PAIRS (INPUT). C NMAX = MAXIMUM NUMBER OF ALLOWED FUNCTION VALUES (OUTPUT). C C IDAT (OUTPUT) C IDAT EQUAL TO 1 --> FUNCTION *** FUN *** IS GIVEN ANALYTICALLY C IDAT EQUAL TO 2 --> FUNCTION *** FUN *** IS GIVEN AS TABULATED C DATA IN A FILE C IDAT EQUAL TO 3 --> FUNCTION *** FUN *** IS GIVEN AS TABULATED C DATA PASSED BY CALL TO *** HANKEL *** C C IERR4 = ERROR FLAG (OUTPUT) C IERR4 EQUAL TO 0 ON RETURN --> NO ERROR OCCURED. C IERR4 EQUAL TO 1 ON RETURN --> AN ERROR OCCURED DURING C EVALUATION OF FUNCTION *** FUN *** C C ################################################################## C C IMPLICIT NONE C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C FUNCTION *** FUN *** IS GIVEN AS DATA TABLE IN AN INPUT FILE. C THE FORMAT OF THE DATA IS: C X_1 Y_1 C ... ... C X_N Y_N C ... ... C X_NEND Y_NEND C I.E. EACH LINE CONTAINS ONE DATA POINT (X,Y). C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C .. Parameters .. INTEGER NMAX PARAMETER (NMAX=1000) INTEGER ITEST PARAMETER (ITEST=0) C .. C .. Scalar Arguments .. DOUBLE PRECISION X,Y INTEGER IERR4,NENDIN C .. C .. Array Arguments .. DOUBLE PRECISION XDATIN(NMAX),YDATIN(NMAX) C .. C .. Scalars in Common .. DOUBLE PRECISION NU,XDATMA,XDATMI,XI,XLIMIT,YDATMA,YDATMI INTEGER IDAT,IREAD,LUERR,LUIN,LUOUT,NEND C .. C .. Arrays in Common .. DOUBLE PRECISION XDAT(NMAX),YDAT(NMAX) C .. C .. Local Scalars .. DOUBLE PRECISION RESULT INTEGER IDEG,IERR6,N C .. C .. Local Arrays .. DOUBLE PRECISION TABLE(NMAX,NMAX) C .. C .. External Functions .. DOUBLE PRECISION D1MACH EXTERNAL D1MACH C .. C .. External Subroutines .. EXTERNAL DIVDIF C .. C .. Intrinsic Functions .. INTRINSIC MAX,MIN C .. C .. Common blocks .. COMMON /ARGUME/XI,NU,XLIMIT COMMON /DATEN/XDAT,YDAT,XDATMI,XDATMA,YDATMI,YDATMA,NEND COMMON /INOUTE/LUIN,LUOUT,LUERR COMMON /LOGIK/IDAT,IREAD integer i1mach, nout nout = i1mach(2) C .. IERR4 = 0 C C ACCEPT DATA FROM CALLING PROGRAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C IF (IREAD.NE.1) THEN C C SAVE INPUT DATA C NEND = NENDIN DO 10 N = 1,NEND XDAT(N) = XDATIN(N) YDAT(N) = YDATIN(N) 10 CONTINUE C WRITE (LUERR,FMT=9000) C READ DATA ONLY ONCE: IREAD = 1 C C CHECK DATA C IF (NEND.GT.NMAX) THEN WRITE(NOUT,FMT=9010) NEND STOP END IF C XDATMI = D1MACH(2) XDATMA = -D1MACH(2) YDATMI = D1MACH(2) YDATMA = -D1MACH(2) DO 20 N = 1,NEND XDATMI = MIN(XDATMI,XDAT(N)) XDATMA = MAX(XDATMA,XDAT(N)) YDATMI = MIN(YDATMI,YDAT(N)) YDATMA = MAX(YDATMA,YDAT(N)) 20 CONTINUE C IF (XDATMI.LT.-D1MACH(2)) THEN WRITE(NOUT,FMT=9020) WRITE (LUERR,FMT=9020) STOP END IF IF (XDATMI.GT.D1MACH(2)) THEN WRITE(NOUT,FMT=9030) WRITE (LUERR,FMT=9030) STOP END IF IF (YDATMI.LT.-D1MACH(2)) THEN WRITE(NOUT,FMT=9040) WRITE (LUERR,FMT=9040) STOP END IF IF (YDATMI.GT.D1MACH(2)) THEN WRITE(NOUT,FMT=9050) WRITE (LUERR,FMT=9050) STOP END IF C C SUCCESSFUL RETURN FROM *** DHFUNC3 ***, DATA HAVE BEEN READ! IDAT = 3 RETURN C C END OF IF-CLAUSE *** IREAD *** END IF C IF (ITEST.EQ.1) THEN WRITE(NOUT,FMT=9060) DO 30 N = 1,NEND WRITE(NOUT,FMT=*) N,XDAT(N),YDAT(N) 30 CONTINUE END IF C C INTERPOLATE FUNCTION VALUE FOR ARGUMENT *** X *** !!!!!!!!!!!!!!!! C C DO NOT PERFORM THE DEFINITE INTEGRAL ABOVE THE TABULATED C DATA RANGE: XLIMIT = XDATMA C IF (X.GE.XDATMI .AND. X.LE.XDATMA) THEN C C *** X*** LIES WITHIN TABULATED DATA RANGE. TRY THE INTERPOLATION! C C CHOOSE *** IDEG *** NOT TOO HIGH!!! IDEG = 5 IERR6 = 0 C CALL DIVDIF(XDAT,YDAT,NEND,IDEG,IDEG,TABLE,X,RESULT,IERR6) C CALL DIVDIF (XDAT,YDAT,NEND,IDEG,IDEG,X,RESULT,IERR6) C IF (IERR6.NE.0) THEN WRITE (LUERR,FMT=9070) IERR6 WRITE(NOUT,FMT=9070) IERR6 WRITE (LUERR,FMT=9080) X WRITE (LUERR,FMT=9090) NEND IERR4 = 1 WRITE (LUERR,FMT=9100) IERR4 RETURN END IF C C SUCCESSFUL RETURN FROM *** DIVDIF *** Y = RESULT C ELSE C C *** X *** LIES OUTSIDE THE TABULATED DATA RANGE! C NO INTERPOLATION POSSIBLE! C THIS CASE WILL LEGALLY HAPPEN, SINCE THE HANKEL TRANSFORM C INCLUDES AN INTEGRATION TO INFINITY. C Y = 0.0D0 C C END IF-CLAUSE *** X.LE.XDATMA ***. END IF C RETURN 9000 FORMAT (/,' MESSAGE FUNK3: ',/, + ' DATA ** X(N), Y(N) *** HAVE BEE N READ FROM CALLING PROGRAM!' + ) 9010 FORMAT (/,' MELDUNG RW:',/,' ERROR! MAXIMAL NUMBER =',I6, + ' OF DATA PAIRS EXCEEDED!') 9020 FORMAT (/,' MESSAGE RW: TOO SMALL VALUE FOR ABCISSA *** X ***!',/, + ' PLEASE CHECK YOUR DATA!') 9030 FORMAT (/,' MESSAGE FUNK3: TOO BIG VALUE FOR ABCISSA *** X ***!', + /,' PLEASE CHECK YOUR DATA!') 9040 FORMAT (/, + ' MESSAGE FUNK3: TOO SMALL VALUE FOR ORDINATE *** Y ** *!' + ,/,' PLEASE CHECK YOUR DATA!') 9050 FORMAT (/, + ' MESSAGE FUNK3: TOO BIG VALUE FOR ORDINATE *** X ***! ', + /,' PLEASE CHECK YOUR DATA!') 9060 FORMAT (' MESSAGE FUNK3: DATA RECEIVED FROM CALLING PROGRAM:') 9070 FORMAT (/,' MESSAGE FUNK3:',/, + ' ERROR ON INTERPOLATING FUNCTION VALUE!',/, + ' ERROR IN SUBROUTINE *** DIVDIF ***!',/,' IERR6 = ',I6) 9080 FORMAT (1X,' ARGUMENT *** X *** =',D17.9) 9090 FORMAT (1X,' NUMBER OF DATA PAIRS *** NEND *** =',I6) 9100 FORMAT (/,' MESSAGE FUNK3: IERR4 =',I6) END SHAR_EOF fi # end of overwriting check if test -f 'res' then echo shar: will not over-write existing file "'res'" else cat << SHAR_EOF > 'res' # XI HANKEL_TRANSFORM(XI) 0.500000000E+01 0.208245531E-01 OUTINE TO H A N K E L. AS AN EXAMPLE, WE WILL NUMERICALLY HANKEL TRANSFORM FUN(X) = X**(NU+1/2) FOR X < 1 AND = 0 ELSEWHERE AT POSITION XI FOR TRANSFORMATION ORDER NU. THE ANALYTICAL SOLUTION IS HT(FUN(X),NU,XI) = XI**(-1/2)*J(NU+1,XI) (J(NU,X) IS THE BESSEL FUNCTION OF THE FIRST KIND.) MESSAGE DRIVER: GIVE NAME FOR PROTOCOL FILE (E.G. hankel.err): MESSAGE DRIVER: ERROR LOG FILE =res.err MESSAGE DRIVER: GIVE NAME FOR OUTPUT FILE (E.G. hankel.out): MESSAGE DRIVER: RESULT FILE =res FUNCTION *** FUN(X) *** IS EITHER CODED WITHIN FUNCTION *** DHFUNC *** = 1 OR TABULATED IN A FILE = 2 (VERY SLOW!) OR PASSED AS TABLE VIA CALL TO HANKEL = 3 (SLOW!) OR STOP = 0 GIVE CHOICE (1 OR 2 OR 3 OR 0): GIVE THE ORDER *** NU ***: (E.G. 1.0): GIVE THE ARGUMENT *** XI ***: (E.G. 5.0): #################################### MESSAGE HANKEL: BEGIN OF CALCULATION #################################### *** PROGRAM HANKEL *** WRITTEN BY THOMAS WIEDER E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE DATE OF FIRST VERSION: A.D. 1998 DATE OF LAST CHANGE: FEBRUARY 1999 MESSAGE HANKEL: WILL TRY TO TRANSFORM *** FUNCTION FUNC ***! MESSAGE HANKEL: WILL TRY NUMERICAL HANKEL TRANSFORM...PLEASE WAIT! MESSAGE HANKEL: ARGUMENT XI ORDER NU NUMERICAL HANKEL TRANSFORM HT: 0.500000000D+01 0.100000000D+01 0.208245531D-01 MESSAGE HANKEL: FOR ERROR MESSAGES LOOK INTO FILE: res.err MESSAGE HANKEL: FOR RESULTS LOOK INTO FILE: res IERR1 = 0 MESSAGE HANKEL: SUCCESSFUL END OF CALCULATION! NO ERRORS HAVE BEEN DETECTED DURING CALCULATION! GOOD BYE! DO YOU WANT ANOTHER CALCULATION = 1 STOP = 2 GIVE 1 OR 2: A PROTOCOL HAS BEEN WRITTEN INTO A FILE, RESULTS HAVE BEEN WRITTEN INTO A FILE. GOODBYE! SHAR_EOF fi # end of overwriting check if test -f 'res.err' then echo shar: will not over-write existing file "'res.err'" else cat << SHAR_EOF > 'res.err' MESSAGE HANKEL: ARGUMENT XI ORDER NU NUMERICAL HANKEL TRANSFORM HT: 0.500000000D+01 0.100000000D+01 0.208245531D-01 SHAR_EOF fi # end of overwriting check cd .. cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'specfun.f' then echo shar: will not over-write existing file "'specfun.f'" else cat << SHAR_EOF > 'specfun.f' C BEGIN OF SUBROUTINE *** RJBESL *** SUBROUTINE RJBESL(XINPUT,ALPHA,NB,B,NCALC) C C ################################################################## C C *** T. WIEDER, 05.03.1998 *** C C SOURCE FOR SUBROUTINE *** RJBESL ***: C C SUBPROGRAM LIBRARY SPECFUN C C AVAILABILITIY: PUBLIC DOMAIN C C DEVELOPER: W.J. CODY AND L. STOLTZ, APPLIED MATHEMATICS DIVISION, C ARGONNE NATIONAL LABORATORY, ARGONNE, IL 60439. C C DISTRIBUTOR: NETLIB C C ################################################################## C C--------------------------------------------------------------------- C This routine calculates Bessel functions J sub(N+ALPHA) (X) C for non-negative argument X, and non-negative order N+ALPHA. C C C Explanation of variables in the calling sequence. C C XIN = X *** T. WIEDER, 15.04.98 *** C X - working precision non-negative real argument for which C J's are to be calculated. C ALPHA - working precision fractional part of order for which C J's or exponentially scaled J'r (J*exp(X)) are C to be calculated. 0 <= ALPHA < 1.0. C NB - integer number of functions to be calculated, NB > 0. C The first function calculated is of order ALPHA, and the C last is of order (NB - 1 + ALPHA). C B - working precision output vector of length NB. If RJBESL C terminates normally (NCALC=NB), the vector B contains the C functions J/ALPHA/(X) through J/NB-1+ALPHA/(X), or the C corresponding exponentially scaled functions. C NCALC - integer output variable indicating possible errors. C Before using the vector B, the user should check that C NCALC=NB, i.e., all orders have been calculated to C the desired accuracy. See Error Returns below. C C C******************************************************************* C******************************************************************* C C Explanation of machine-dependent constants C C it = Number of bits in the mantissa of a working precision C variable C NSIG = Decimal significance desired. Should be set to C INT(LOG10(2)*it+1). Setting NSIG lower will result C in decreased accuracy while setting NSIG higher will C increase CPU time without increasing accuracy. The C truncation error is limited to a relative error of C T=.5*10**(-NSIG). C ENTEN = 10.0 ** K, where K is the largest integer such that C ENTEN is machine-representable in working precision C ENSIG = 10.0 ** NSIG C RTNSIG = 10.0 ** (-K) for the smallest integer K such that C K .GE. NSIG/4 C ENMTEN = Smallest ABS(X) such that X/4 does not underflow C XLARGE = Upper limit on the magnitude of X. If ABS(X)=N, C then at least N iterations of the backward recursion C will be executed. The value of 10.0 ** 4 is used on C every machine. C C C Approximate values for some important machines are: C C C it NSIG ENTEN ENSIG C C CRAY-1 (S.P.) 48 15 1.0E+2465 1.0E+15 C Cyber 180/855 C under NOS (S.P.) 48 15 1.0E+322 1.0E+15 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 24 8 1.0E+38 1.0E+8 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 53 16 1.0D+308 1.0D+16 C IBM 3033 (D.P.) 14 5 1.0D+75 1.0D+5 C VAX (S.P.) 24 8 1.0E+38 1.0E+8 C VAX D-Format (D.P.) 56 17 1.0D+38 1.0D+17 C VAX G-Format (D.P.) 53 16 1.0D+307 1.0D+16 C C C RTNSIG ENMTEN XLARGE C C CRAY-1 (S.P.) 1.0E-4 1.84E-2466 1.0E+4 C Cyber 180/855 C under NOS (S.P.) 1.0E-4 1.25E-293 1.0E+4 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 1.0E-2 4.70E-38 1.0E+4 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 1.0E-4 8.90D-308 1.0D+4 C IBM 3033 (D.P.) 1.0E-2 2.16D-78 1.0D+4 C VAX (S.P.) 1.0E-2 1.17E-38 1.0E+4 C VAX D-Format (D.P.) 1.0E-5 1.17D-38 1.0D+4 C VAX G-Format (D.P.) 1.0E-4 2.22D-308 1.0D+4 C C******************************************************************* C******************************************************************* C C Error returns C C In case of an error, NCALC .NE. NB, and not all J's are C calculated to the desired accuracy. C C NCALC .LT. 0: An argument is out of range. For example, C NBES .LE. 0, ALPHA .LT. 0 or .GT. 1, or X is too large. C In this case, B(1) is set to zero, the remainder of the C B-vector is not calculated, and NCALC is set to C MIN(NB,0)-1 so that NCALC .NE. NB. C C NB .GT. NCALC .GT. 0: Not all requested function values could C be calculated accurately. This usually occurs because NB is C much larger than ABS(X). In this case, B(N) is calculated C to the desired accuracy for N .LE. NCALC, but precision C is lost for NCALC .LT. N .LE. NB. If B(N) does not vanish C for N .GT. NCALC (because it is too small to be represented), C and B(N)/B(NCALC) = 10**(-K), then only the first NSIG-K C significant figures of B(N) can be trusted. C C C Intrinsic and other functions required are: C C ABS, AINT, COS, DBLE, GAMMA (or DGAMMA), INT, MAX, MIN, C C REAL, SIN, SQRT C C C Acknowledgement C C This program is based on a program written by David J. Sookne C (2) that computes values of the Bessel functions J or I of real C argument and integer order. Modifications include the restriction C of the computation to the J Bessel function of non-negative real C argument, the extension of the computation to arbitrary positive C order, and the elimination of most underflow. C C References: "A Note on Backward Recurrence Algorithms," Olver, C F. W. J., and Sookne, D. J., Math. Comp. 26, 1972, C pp 941-947. C C "Bessel Functions of Real Argument and Integer Order," C Sookne, D. J., NBS Jour. of Res. B. 77B, 1973, pp C 125-132. C C Latest modification: March 19, 1990 C C Author: W. J. Cody C Applied Mathematics Division C Argonne National Laboratory C Argonne, IL 60439 C C--------------------------------------------------------------------- CS REAL GAMMA, C .. Scalar Arguments .. DOUBLE PRECISION ALPHA,XINPUT INTEGER NB,NCALC C .. C .. Array Arguments .. DOUBLE PRECISION B(NB) C .. C .. Local Scalars .. DOUBLE PRECISION ALP2EM,ALPEM,CAPP,CAPQ,EIGHTH,EM,EN,ENMTEN,ENSIG, + ENTEN,FOUR,GNU,HALF,HALFX,ONE,ONE30,P,PI2,PLAST, + POLD,PSAVE,PSAVEL,RTNSIG,S,SUM,T,T1,TEMPA,TEMPB, + TEMPC,TEST,THREE,THREE5,TOVER,TWO,TWOFIV,TWOPI1, + TWOPI2,VCOS,VSIN,X,XC,XIN,XK,XLARGE,XM,Z,ZERO INTEGER I,J,K,L,M,MAGX,N,NBMX,NEND,NSTART C .. C .. Local Arrays .. DOUBLE PRECISION FACT(25) C .. C .. External Functions .. DOUBLE PRECISION DGAMMA EXTERNAL DGAMMA C .. C .. Intrinsic Functions .. INTRINSIC ABS,AINT,COS,DBLE,INT,MAX,MIN,SIN,SQRT C .. C .. Statement Functions .. DOUBLE PRECISION CONV,FUNC C .. C .. Data statements .. C--------------------------------------------------------------------- C Mathematical constants C C PI2 - 2 / PI C TWOPI1 - first few significant digits of 2 * PI C TWOPI2 - (2*PI - TWOPI) to working precision, i.e., C TWOPI1 + TWOPI2 = 2 * PI to extra precision. C--------------------------------------------------------------------- C DATA PI2, TWOPI1, TWOPI2 /0.636619772367581343075535E0,6.28125E0, C 1 1.935307179586476925286767E-3/ C DATA ZERO, EIGHTH, HALF, ONE /0.0E0,0.125E0,0.5E0,1.0E0/ C DATA TWO, THREE, FOUR, TWOFIV /2.0E0,3.0E0,4.0E0,25.0E0/ C DATA ONE30, THREE5 /130.0E0,35.0E0/ C--------------------------------------------------------------------- C Machine-dependent parameters C--------------------------------------------------------------------- CS DATA ENTEN, ENSIG, RTNSIG /1.0E38,1.0E8,1.0E-2/ CS DATA ENMTEN, XLARGE /1.2E-37,1.0E4/ C--------------------------------------------------------------------- C Factorial(N) C--------------------------------------------------------------------- CS DATA FACT /1.0E0,1.0E0,2.0E0,6.0E0,24.0E0,1.2E2,7.2E2,5.04E3, CS 1 4.032E4,3.6288E5,3.6288E6,3.99168E7,4.790016E8,6.2270208E9, CS 2 8.71782912E10,1.307674368E12,2.0922789888E13,3.55687428096E14, CS 3 6.402373705728E15,1.21645100408832E17,2.43290200817664E18, CS 4 5.109094217170944E19,1.12400072777760768E21, CS 5 2.585201673888497664E22,6.2044840173323943936E23/ DATA PI2,TWOPI1,TWOPI2/0.636619772367581343075535D0,6.28125D0, + 1.935307179586476925286767D-3/ DATA ZERO,EIGHTH,HALF,ONE/0.0D0,0.125D0,0.5D0,1.0D0/ DATA TWO,THREE,FOUR,TWOFIV/2.0D0,3.0D0,4.0D0,25.0D0/ DATA ONE30,THREE5/130.0D0,35.0D0/ DATA ENTEN,ENSIG,RTNSIG/1.0D38,1.0D17,1.0D-4/ DATA ENMTEN,XLARGE/1.2D-37,1.0D4/ DATA FACT/1.0D0,1.0D0,2.0D0,6.0D0,24.0D0,1.2D2,7.2D2,5.04D3, + 4.032D4,3.6288D5,3.6288D6,3.99168D7,4.790016D8,6.2270208D9, + 8.71782912D10,1.307674368D12,2.0922789888D13, + 3.55687428096D14,6.402373705728D15,1.21645100408832D17, + 2.43290200817664D18,5.109094217170944D19, + 1.12400072777760768D21,2.585201673888497664D22, + 6.2044840173323943936D23/ C .. C .. Statement Function definitions .. C C--------------------------------------------------------------------- C Statement functions for conversion and the gamma function. C--------------------------------------------------------------------- CS CONV(I) = REAL(I) CS FUNC(X) = GAMMA(X) CONV(I) = DBLE(I) FUNC(X) = DGAMMA(X) C .. C C *** T. WIEDER, 15.04.1998 *** C *** AMENDMENT SUGGESTED BY DR. T.R. HOPKINS, 15.04.1998 *** X = XINPUT C C--------------------------------------------------------------------- C Check for out of range arguments. C--------------------------------------------------------------------- MAGX = INT(X) IF ((NB.GT.0) .AND. (X.GE.ZERO) .AND. (X.LE.XLARGE) .AND. + (ALPHA.GE.ZERO) .AND. (ALPHA.LT.ONE)) THEN C--------------------------------------------------------------------- C Initialize result array to zero. C--------------------------------------------------------------------- NCALC = NB DO 10 I = 1,NB B(I) = ZERO 10 CONTINUE C--------------------------------------------------------------------- C Branch to use 2-term ascending series for small X and asymptotic C form for large X when NB is not too large. C--------------------------------------------------------------------- IF (X.LT.RTNSIG) THEN C--------------------------------------------------------------------- C Two-term ascending series for small X. C--------------------------------------------------------------------- TEMPA = ONE ALPEM = ONE + ALPHA HALFX = ZERO IF (X.GT.ENMTEN) HALFX = HALF*X IF (ALPHA.NE.ZERO) TEMPA = HALFX**ALPHA/ + (ALPHA*FUNC(ALPHA)) TEMPB = ZERO IF ((X+ONE).GT.ONE) TEMPB = -HALFX*HALFX B(1) = TEMPA + TEMPA*TEMPB/ALPEM IF ((X.NE.ZERO) .AND. (B(1).EQ.ZERO)) NCALC = 0 IF (NB.NE.1) THEN IF (X.LE.ZERO) THEN DO 20 N = 2,NB B(N) = ZERO 20 CONTINUE ELSE C--------------------------------------------------------------------- C Calculate higher order functions. C--------------------------------------------------------------------- TEMPC = HALFX TOVER = (ENMTEN+ENMTEN)/X IF (TEMPB.NE.ZERO) TOVER = ENMTEN/TEMPB DO 30 N = 2,NB TEMPA = TEMPA/ALPEM ALPEM = ALPEM + ONE TEMPA = TEMPA*TEMPC IF (TEMPA.LE.TOVER*ALPEM) TEMPA = ZERO B(N) = TEMPA + TEMPA*TEMPB/ALPEM IF ((B(N).EQ.ZERO) .AND. + (NCALC.GT.N)) NCALC = N - 1 30 CONTINUE END IF END IF ELSE IF ((X.GT.TWOFIV) .AND. (NB.LE.MAGX+1)) THEN C--------------------------------------------------------------------- C Asymptotic series for X .GT. 21.0. C--------------------------------------------------------------------- XC = SQRT(PI2/X) XIN = (EIGHTH/X)**2 M = 11 IF (X.GE.THREE5) M = 8 IF (X.GE.ONE30) M = 4 XM = FOUR*CONV(M) C--------------------------------------------------------------------- C Argument reduction for SIN and COS routines. C--------------------------------------------------------------------- T = AINT(X/ (TWOPI1+TWOPI2)+HALF) Z = ((X-T*TWOPI1)-T*TWOPI2) - (ALPHA+HALF)/PI2 VSIN = SIN(Z) VCOS = COS(Z) GNU = ALPHA + ALPHA DO 50 I = 1,2 S = ((XM-ONE)-GNU)* ((XM-ONE)+GNU)*XIN*HALF T = (GNU- (XM-THREE))* (GNU+ (XM-THREE)) CAPP = S*T/FACT(2*M+1) T1 = (GNU- (XM+ONE))* (GNU+ (XM+ONE)) CAPQ = S*T1/FACT(2*M+2) XK = XM K = M + M T1 = T DO 40 J = 2,M XK = XK - FOUR S = ((XK-ONE)-GNU)* ((XK-ONE)+GNU) T = (GNU- (XK-THREE))* (GNU+ (XK-THREE)) CAPP = (CAPP+ONE/FACT(K-1))*S*T*XIN CAPQ = (CAPQ+ONE/FACT(K))*S*T1*XIN K = K - 2 T1 = T 40 CONTINUE CAPP = CAPP + ONE CAPQ = (CAPQ+ONE)* (GNU*GNU-ONE)* (EIGHTH/X) B(I) = XC* (CAPP*VCOS-CAPQ*VSIN) IF (NB.EQ.1) GO TO 180 T = VSIN VSIN = -VCOS VCOS = T GNU = GNU + TWO 50 CONTINUE C--------------------------------------------------------------------- C If NB .GT. 2, compute J(X,ORDER+I) I = 2, NB-1 C--------------------------------------------------------------------- IF (NB.GT.2) THEN GNU = ALPHA + ALPHA + TWO DO 60 J = 3,NB B(J) = GNU*B(J-1)/X - B(J-2) GNU = GNU + TWO 60 CONTINUE END IF C--------------------------------------------------------------------- C Use recurrence to generate results. First initialize the C calculation of P*S. C--------------------------------------------------------------------- ELSE NBMX = NB - MAGX N = MAGX + 1 EN = CONV(N+N) + (ALPHA+ALPHA) PLAST = ONE P = EN/X C--------------------------------------------------------------------- C Calculate general significance test. C--------------------------------------------------------------------- TEST = ENSIG + ENSIG IF (NBMX.GE.3) THEN C--------------------------------------------------------------------- C Calculate P*S until N = NB-1. Check for possible overflow. C--------------------------------------------------------------------- TOVER = ENTEN/ENSIG NSTART = MAGX + 2 NEND = NB - 1 EN = CONV(NSTART+NSTART) - TWO + (ALPHA+ALPHA) DO 90 K = NSTART,NEND N = K EN = EN + TWO POLD = PLAST PLAST = P P = EN*PLAST/X - POLD IF (P.GT.TOVER) THEN C--------------------------------------------------------------------- C To avoid overflow, divide P*S by TOVER. Calculate P*S until C ABS(P) .GT. 1. C--------------------------------------------------------------------- TOVER = ENTEN P = P/TOVER PLAST = PLAST/TOVER PSAVE = P PSAVEL = PLAST NSTART = N + 1 70 N = N + 1 EN = EN + TWO POLD = PLAST PLAST = P P = EN*PLAST/X - POLD IF (P.LE.ONE) GO TO 70 TEMPB = EN/X C--------------------------------------------------------------------- C Calculate backward test and find NCALC, the highest N such that C the test is passed. C--------------------------------------------------------------------- TEST = POLD*PLAST* (HALF-HALF/ (TEMPB*TEMPB)) TEST = TEST/ENSIG P = PLAST*TOVER N = N - 1 EN = EN - TWO NEND = MIN(NB,N) DO 80 L = NSTART,NEND POLD = PSAVEL PSAVEL = PSAVE PSAVE = EN*PSAVEL/X - POLD IF (PSAVE*PSAVEL.GT.TEST) THEN NCALC = L - 1 GO TO 110 END IF 80 CONTINUE NCALC = NEND GO TO 110 END IF 90 CONTINUE N = NEND EN = CONV(N+N) + (ALPHA+ALPHA) C--------------------------------------------------------------------- C Calculate special significance test for NBMX .GT. 2. C--------------------------------------------------------------------- TEST = MAX(TEST,SQRT(PLAST*ENSIG)*SQRT(P+P)) END IF C--------------------------------------------------------------------- C Calculate P*S until significance test passes. C--------------------------------------------------------------------- 100 N = N + 1 EN = EN + TWO POLD = PLAST PLAST = P P = EN*PLAST/X - POLD IF (P.LT.TEST) GO TO 100 C--------------------------------------------------------------------- C Initialize the backward recursion and the normalization sum. C--------------------------------------------------------------------- 110 N = N + 1 EN = EN + TWO TEMPB = ZERO TEMPA = ONE/P M = 2*N - 4* (N/2) SUM = ZERO EM = CONV(N/2) ALPEM = (EM-ONE) + ALPHA ALP2EM = (EM+EM) + ALPHA IF (M.NE.0) SUM = TEMPA*ALPEM*ALP2EM/EM NEND = N - NB IF (NEND.GT.0) THEN C--------------------------------------------------------------------- C Recur backward via difference equation, calculating (but not C storing) B(N), until N = NB. C--------------------------------------------------------------------- DO 120 L = 1,NEND N = N - 1 EN = EN - TWO TEMPC = TEMPB TEMPB = TEMPA TEMPA = (EN*TEMPB)/X - TEMPC M = 2 - M IF (M.NE.0) THEN EM = EM - ONE ALP2EM = (EM+EM) + ALPHA IF (N.EQ.1) GO TO 130 ALPEM = (EM-ONE) + ALPHA IF (ALPEM.EQ.ZERO) ALPEM = ONE SUM = (SUM+TEMPA*ALP2EM)*ALPEM/EM END IF 120 CONTINUE END IF C--------------------------------------------------------------------- C Store B(NB). C--------------------------------------------------------------------- 130 B(N) = TEMPA IF (NEND.GE.0) THEN IF (NB.LE.1) THEN ALP2EM = ALPHA IF ((ALPHA+ONE).EQ.ONE) ALP2EM = ONE SUM = SUM + B(1)*ALP2EM GO TO 160 ELSE C--------------------------------------------------------------------- C Calculate and store B(NB-1). C--------------------------------------------------------------------- N = N - 1 EN = EN - TWO B(N) = (EN*TEMPA)/X - TEMPB IF (N.EQ.1) GO TO 150 M = 2 - M IF (M.NE.0) THEN EM = EM - ONE ALP2EM = (EM+EM) + ALPHA ALPEM = (EM-ONE) + ALPHA IF (ALPEM.EQ.ZERO) ALPEM = ONE SUM = (SUM+B(N)*ALP2EM)*ALPEM/EM END IF END IF END IF NEND = N - 2 IF (NEND.NE.0) THEN C--------------------------------------------------------------------- C Calculate via difference equation and store B(N), until N = 2. C--------------------------------------------------------------------- DO 140 L = 1,NEND N = N - 1 EN = EN - TWO B(N) = (EN*B(N+1))/X - B(N+2) M = 2 - M IF (M.NE.0) THEN EM = EM - ONE ALP2EM = (EM+EM) + ALPHA ALPEM = (EM-ONE) + ALPHA IF (ALPEM.EQ.ZERO) ALPEM = ONE SUM = (SUM+B(N)*ALP2EM)*ALPEM/EM END IF 140 CONTINUE END IF C--------------------------------------------------------------------- C Calculate B(1). C--------------------------------------------------------------------- B(1) = TWO* (ALPHA+ONE)*B(2)/X - B(3) 150 EM = EM - ONE ALP2EM = (EM+EM) + ALPHA IF (ALP2EM.EQ.ZERO) ALP2EM = ONE SUM = SUM + B(1)*ALP2EM C--------------------------------------------------------------------- C Normalize. Divide all B(N) by sum. C--------------------------------------------------------------------- 160 IF ((ALPHA+ONE).NE.ONE) SUM = SUM*FUNC(ALPHA)* + (X*HALF)** (-ALPHA) TEMPA = ENMTEN IF (SUM.GT.ONE) TEMPA = TEMPA*SUM DO 170 N = 1,NB IF (ABS(B(N)).LT.TEMPA) B(N) = ZERO B(N) = B(N)/SUM 170 CONTINUE END IF C--------------------------------------------------------------------- C Error return -- X, NB, or ALPHA is out of range. C--------------------------------------------------------------------- ELSE B(1) = ZERO NCALC = MIN(NB,0) - 1 END IF C--------------------------------------------------------------------- C Exit C--------------------------------------------------------------------- 180 RETURN C ---------- Last line of RJBESL ---------- END C END OF SUBROUTINE *** RJBESL *** C C BEGIN OF SUBROUTINE *** RYBESL *** SUBROUTINE RYBESL(X,ALPHA,NB,BY,NCALC) C C ################################################################## C C *** T. WIEDER, 05.03.1998 *** C C SOURCE FOR *** SUBROUTINE RYBESL ***: C C SUBPROGRAM LIBRARY SPECFUN C C AVAILABILITIY: PUBLIC DOMAIN C C DEVELOPER: W.J. CODY AND L. STOLTZ, APPLIED MATHEMATICS DIVISION, C ARGONNE NATIONAL LABORATORY, ARGONNE, IL 60439. C C DISTRIBUTOR: NETLIB C C ################################################################## C C---------------------------------------------------------------------- C C This routine calculates Bessel functions Y SUB(N+ALPHA) (X) C for non-negative argument X, and non-negative order N+ALPHA. C C C Explanation of variables in the calling sequence C C X - Working precision non-negative real argument for which C Y's are to be calculated. C ALPHA - Working precision fractional part of order for which C Y's are to be calculated. 0 .LE. ALPHA .LT. 1.0. C NB - Integer number of functions to be calculated, NB .GT. 0. C The first function calculated is of order ALPHA, and the C last is of order (NB - 1 + ALPHA). C BY - Working precision output vector of length NB. If the C routine terminates normally (NCALC=NB), the vector BY C contains the functions Y(ALPHA,X), ... , Y(NB-1+ALPHA,X), C If (0 .LT. NCALC .LT. NB), BY(I) contains correct function C values for I .LE. NCALC, and contains the ratios C Y(ALPHA+I-1,X)/Y(ALPHA+I-2,X) for the rest of the array. C NCALC - Integer output variable indicating possible errors. C Before using the vector BY, the user should check that C NCALC=NB, i.e., all orders have been calculated to C the desired accuracy. See error returns below. C C C******************************************************************* C******************************************************************* C C Explanation of machine-dependent constants C C beta = Radix for the floating-point system C p = Number of significant base-beta digits in the C significand of a floating-point number C minexp = Smallest representable power of beta C maxexp = Smallest power of beta that overflows C EPS = beta ** (-p) C DEL = Machine number below which sin(x)/x = 1; approximately C SQRT(EPS). C XMIN = Smallest acceptable argument for RBESY; approximately C max(2*beta**minexp,2/XINF), rounded up C XINF = Largest positive machine number; approximately C beta**maxexp C THRESH = Lower bound for use of the asymptotic form; approximately C AINT(-LOG10(EPS/2.0))+1.0 C XLARGE = Upper bound on X; approximately 1/DEL, because the sine C and cosine functions have lost about half of their C precision at that point. C C C Approximate values for some important machines are: C C beta p minexp maxexp EPS C C CRAY-1 (S.P.) 2 48 -8193 8191 3.55E-15 C Cyber 180/185 C under NOS (S.P.) 2 48 -975 1070 3.55E-15 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 2 24 -126 128 5.96E-8 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 2 53 -1022 1024 1.11D-16 C IBM 3033 (D.P.) 16 14 -65 63 1.39D-17 C VAX (S.P.) 2 24 -128 127 5.96E-8 C VAX D-Format (D.P.) 2 56 -128 127 1.39D-17 C VAX G-Format (D.P.) 2 53 -1024 1023 1.11D-16 C C C DEL XMIN XINF THRESH XLARGE C C CRAY-1 (S.P.) 5.0E-8 3.67E-2466 5.45E+2465 15.0E0 2.0E7 C Cyber 180/855 C under NOS (S.P.) 5.0E-8 6.28E-294 1.26E+322 15.0E0 2.0E7 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 1.0E-4 2.36E-38 3.40E+38 8.0E0 1.0E4 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 1.0D-8 4.46D-308 1.79D+308 16.0D0 1.0D8 C IBM 3033 (D.P.) 1.0D-8 2.77D-76 7.23D+75 17.0D0 1.0D8 C VAX (S.P.) 1.0E-4 1.18E-38 1.70E+38 8.0E0 1.0E4 C VAX D-Format (D.P.) 1.0D-9 1.18D-38 1.70D+38 17.0D0 1.0D9 C VAX G-Format (D.P.) 1.0D-8 2.23D-308 8.98D+307 16.0D0 1.0D8 C C******************************************************************* C******************************************************************* C C Error returns C C In case of an error, NCALC .NE. NB, and not all Y's are C calculated to the desired accuracy. C C NCALC .LT. -1: An argument is out of range. For example, C NB .LE. 0, IZE is not 1 or 2, or IZE=1 and ABS(X) .GE. C XMAX. In this case, BY(1) = 0.0, the remainder of the C BY-vector is not calculated, and NCALC is set to C MIN0(NB,0)-2 so that NCALC .NE. NB. C NCALC = -1: Y(ALPHA,X) .GE. XINF. The requested function C values are set to 0.0. C 1 .LT. NCALC .LT. NB: Not all requested function values could C be calculated accurately. BY(I) contains correct function C values for I .LE. NCALC, and and the remaining NB-NCALC C array elements contain 0.0. C C C Intrinsic functions required are: C C DBLE, EXP, INT, MAX, MIN, REAL, SQRT C C C Acknowledgement C C This program draws heavily on Temme's Algol program for Y(a,x) C and Y(a+1,x) and on Campbell's programs for Y_nu(x). Temme's C scheme is used for x < THRESH, and Campbell's scheme is used C in the asymptotic region. Segments of code from both sources C have been translated into Fortran 77, merged, and heavily modified. C Modifications include parameterization of machine dependencies, C use of a new approximation for ln(gamma(x)), and built-in C protection against over/underflow. C C References: "Bessel functions J_nu(x) and Y_nu(x) of real C order and real argument," Campbell, J. B., C Comp. Phy. Comm. 18, 1979, pp. 133-142. C C "On the numerical evaluation of the ordinary C Bessel function of the second kind," Temme, C N. M., J. Comput. Phys. 21, 1976, pp. 343-350. C C Latest modification: March 19, 1990 C C Modified by: W. J. Cody C Applied Mathematics Division C Argonne National Laboratory C Argonne, IL 60439 C C---------------------------------------------------------------------- CS REAL C DIMENSION BY(NB), CH(21) C *** T. WIEDER, 07.02.1999 *** C .. Scalar Arguments .. DOUBLE PRECISION ALPHA,X INTEGER NB,NCALC C .. C .. Array Arguments .. DOUBLE PRECISION BY(NB+1) C .. C .. Local Scalars .. DOUBLE PRECISION ALFA,AYE,B,C,COSMU,D,D1,D2,DDIV,DEL,DEN,DIV,DMU, + E,EIGHT,EN,EN1,ENU,EPS,EVEN,EX,F,FIVPI,G,GAMMA,H, + HALF,ODD,ONBPI,ONE,ONE5,P,PA,PA1,PI,PIBY2,PIM5,Q, + Q0,QA,QA1,R,S,SINMU,SQ2BPI,TEN9,TERM,THREE, + THRESH,TWO,TWOBYX,X2,XINF,XLARGE,XMIN,XNA,YA,YA1, + ZERO INTEGER I,K,NA C .. C .. Local Arrays .. DOUBLE PRECISION CH(21) C .. C .. Intrinsic Functions .. INTRINSIC ABS,AINT,COS,INT,LOG,MIN,SIN,SQRT C .. C .. Data statements .. C---------------------------------------------------------------------- C Mathematical constants C FIVPI = 5*PI C PIM5 = 5*PI - 15 C ONBPI = 1/PI C PIBY2 = PI/2 C SQ2BPI = SQUARE ROOT OF 2/PI C---------------------------------------------------------------------- CS DATA ZERO,HALF,ONE,TWO,THREE/0.0E0,0.5E0,1.0E0,2.0E0,3.0E0/ CS DATA EIGHT,ONE5,TEN9/8.0E0,15.0E0,1.9E1/ CS DATA FIVPI,PIBY2/1.5707963267948966192E1,1.5707963267948966192E0/ CS DATA PI,SQ2BPI/3.1415926535897932385E0,7.9788456080286535588E-1/ CS DATA PIM5,ONBPI/7.0796326794896619231E-1,3.1830988618379067154E-1/ C---------------------------------------------------------------------- C Machine-dependent constants C---------------------------------------------------------------------- CS DATA DEL,XMIN,XINF,EPS/1.0E-4,2.36E-38,3.40E38,5.96E-8/ CS DATA THRESH,XLARGE/8.0E0,1.0E4/ C *** T. WIEDER, 11.03.1998 *** C DATA DEL,XMIN,XINF,EPS/1.0D-8,4.46D-308,1.79D308,1.11D-16/ C---------------------------------------------------------------------- C Coefficients for Chebyshev polynomial expansion of C 1/gamma(1-x), abs(x) .le. .5 C---------------------------------------------------------------------- CS DATA CH/-0.67735241822398840964E-23,-0.61455180116049879894E-22, CS 1 0.29017595056104745456E-20, 0.13639417919073099464E-18, CS 2 0.23826220476859635824E-17,-0.90642907957550702534E-17, CS 3 -0.14943667065169001769E-14,-0.33919078305362211264E-13, CS 4 -0.17023776642512729175E-12, 0.91609750938768647911E-11, CS 5 0.24230957900482704055E-09, 0.17451364971382984243E-08, CS 6 -0.33126119768180852711E-07,-0.86592079961391259661E-06, CS 7 -0.49717367041957398581E-05, 0.76309597585908126618E-04, CS 8 0.12719271366545622927E-02, 0.17063050710955562222E-02, CS 9 -0.76852840844786673690E-01,-0.28387654227602353814E+00, CS A 0.92187029365045265648E+00/ DATA ZERO,HALF,ONE,TWO,THREE/0.0D0,0.5D0,1.0D0,2.0D0,3.0D0/ DATA EIGHT,ONE5,TEN9/8.0D0,15.0D0,1.9D1/ DATA FIVPI,PIBY2/1.5707963267948966192D1,1.5707963267948966192D0/ DATA PI,SQ2BPI/3.1415926535897932385D0,7.9788456080286535588D-1/ DATA PIM5,ONBPI/7.0796326794896619231D-1,3.1830988618379067154D-1/ DATA DEL,XMIN,XINF,EPS/1.0D-8,2.23D-307,1.79D308,1.11D-16/ DATA THRESH,XLARGE/16.0D0,1.0D8/ DATA CH/-0.67735241822398840964D-23,-0.61455180116049879894D-22, + 0.29017595056104745456D-20,0.13639417919073099464D-18, + 0.23826220476859635824D-17,-0.90642907957550702534D-17, + -0.14943667065169001769D-14,-0.33919078305362211264D-13, + -0.17023776642512729175D-12,0.91609750938768647911D-11, + 0.24230957900482704055D-09,0.17451364971382984243D-08, + -0.33126119768180852711D-07,-0.86592079961391259661D-06, + -0.49717367041957398581D-05,0.76309597585908126618D-04, + 0.12719271366545622927D-02,0.17063050710955562222D-02, + -0.76852840844786673690D-01,-0.28387654227602353814D+00, + 0.92187029365045265648D+00/ C .. C---------------------------------------------------------------------- EX = X ENU = ALPHA IF ((NB.GT.0) .AND. (X.GE.XMIN) .AND. (EX.LT.XLARGE) .AND. + (ENU.GE.ZERO) .AND. (ENU.LT.ONE)) THEN XNA = AINT(ENU+HALF) NA = INT(XNA) IF (NA.EQ.1) ENU = ENU - XNA IF (ENU.EQ.-HALF) THEN P = SQ2BPI/SQRT(EX) YA = P*SIN(EX) YA1 = -P*COS(EX) ELSE IF (EX.LT.THREE) THEN C---------------------------------------------------------------------- C Use Temme's scheme for small X C---------------------------------------------------------------------- B = EX*HALF D = -LOG(B) F = ENU*D E = B** (-ENU) IF (ABS(ENU).LT.DEL) THEN C = ONBPI ELSE C = ENU/SIN(ENU*PI) END IF C---------------------------------------------------------------------- C Computation of sinh(f)/f C---------------------------------------------------------------------- IF (ABS(F).LT.ONE) THEN X2 = F*F EN = TEN9 S = ONE DO 10 I = 1,9 S = S*X2/EN/ (EN-ONE) + ONE EN = EN - TWO 10 CONTINUE ELSE S = (E-ONE/E)*HALF/F END IF C---------------------------------------------------------------------- C Computation of 1/gamma(1-a) using Chebyshev polynomials C---------------------------------------------------------------------- X2 = ENU*ENU*EIGHT AYE = CH(1) EVEN = ZERO ALFA = CH(2) ODD = ZERO DO 20 I = 3,19,2 EVEN = - (AYE+AYE+EVEN) AYE = -EVEN*X2 - AYE + CH(I) ODD = - (ALFA+ALFA+ODD) ALFA = -ODD*X2 - ALFA + CH(I+1) 20 CONTINUE EVEN = (EVEN*HALF+AYE)*X2 - AYE + CH(21) ODD = (ODD+ALFA)*TWO GAMMA = ODD*ENU + EVEN C---------------------------------------------------------------------- C End of computation of 1/gamma(1-a) C---------------------------------------------------------------------- G = E*GAMMA E = (E+ONE/E)*HALF F = TWO*C* (ODD*E+EVEN*S*D) E = ENU*ENU P = G*C Q = ONBPI/G C = ENU*PIBY2 IF (ABS(C).LT.DEL) THEN R = ONE ELSE R = SIN(C)/C END IF R = PI*C*R*R C = ONE D = -B*B H = ZERO YA = F + R*Q YA1 = P EN = ZERO 30 EN = EN + ONE IF (ABS(G/ (ONE+ABS(YA)))+ABS(H/ (ONE+ABS(YA1))).GT. + EPS) THEN F = (F*EN+P+Q)/ (EN*EN-E) C = C*D/EN P = P/ (EN-ENU) Q = Q/ (EN+ENU) G = C* (F+R*Q) H = C*P - EN*G YA = YA + G YA1 = YA1 + H GO TO 30 END IF YA = -YA YA1 = -YA1/B ELSE IF (EX.LT.THRESH) THEN C---------------------------------------------------------------------- C Use Temme's scheme for moderate X C---------------------------------------------------------------------- C = (HALF-ENU)* (HALF+ENU) B = EX + EX E = (EX*ONBPI*COS(ENU*PI)/EPS) E = E*E P = ONE Q = -EX R = ONE + EX*EX S = R EN = TWO 40 IF (R*EN*EN.LT.E) THEN EN1 = EN + ONE D = (EN-ONE+C/EN)/S P = (EN+EN-P*D)/EN1 Q = (-B+Q*D)/EN1 S = P*P + Q*Q R = R*S EN = EN1 GO TO 40 END IF F = P/S P = F G = -Q/S Q = G 50 EN = EN - ONE IF (EN.GT.ZERO) THEN R = EN1* (TWO-P) - TWO S = B + EN1*Q D = (EN-ONE+C/EN)/ (R*R+S*S) P = D*R Q = D*S E = F + ONE F = P*E - G*Q G = Q*E + P*G EN1 = EN GO TO 50 END IF F = ONE + F D = F*F + G*G PA = F/D QA = -G/D D = ENU + HALF - P Q = Q + EX PA1 = (PA*Q-QA*D)/EX QA1 = (QA*Q+PA*D)/EX B = EX - PIBY2* (ENU+HALF) C = COS(B) S = SIN(B) D = SQ2BPI/SQRT(EX) YA = D* (PA*S+QA*C) YA1 = D* (QA1*S-PA1*C) ELSE C---------------------------------------------------------------------- C Use Campbell's asymptotic scheme. C---------------------------------------------------------------------- NA = 0 D1 = AINT(EX/FIVPI) I = INT(D1) DMU = ((EX-ONE5*D1)-D1*PIM5) - (ALPHA+HALF)*PIBY2 IF (I-2* (I/2).EQ.0) THEN COSMU = COS(DMU) SINMU = SIN(DMU) ELSE COSMU = -COS(DMU) SINMU = -SIN(DMU) END IF DDIV = EIGHT*EX DMU = ALPHA DEN = SQRT(EX) DO 80 K = 1,2 P = COSMU COSMU = SINMU SINMU = -P D1 = (TWO*DMU-ONE)* (TWO*DMU+ONE) D2 = ZERO DIV = DDIV P = ZERO Q = ZERO Q0 = D1/DIV TERM = Q0