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 DO 60 I = 2,20 D2 = D2 + EIGHT D1 = D1 - D2 DIV = DIV + DDIV TERM = -TERM*D1/DIV P = P + TERM D2 = D2 + EIGHT D1 = D1 - D2 DIV = DIV + DDIV TERM = TERM*D1/DIV Q = Q + TERM IF (ABS(TERM).LE.EPS) GO TO 70 60 CONTINUE 70 P = P + ONE Q = Q + Q0 IF (K.EQ.1) THEN YA = SQ2BPI* (P*COSMU-Q*SINMU)/DEN ELSE YA1 = SQ2BPI* (P*COSMU-Q*SINMU)/DEN END IF DMU = DMU + ONE 80 CONTINUE END IF IF (NA.EQ.1) THEN H = TWO* (ENU+ONE)/EX IF (H.GT.ONE) THEN IF (ABS(YA1).GT.XINF/H) THEN H = ZERO YA = ZERO END IF END IF H = H*YA1 - YA YA = YA1 YA1 = H END IF C---------------------------------------------------------------------- C Now have first one or two Y's C---------------------------------------------------------------------- BY(1) = YA BY(2) = YA1 IF (YA1.EQ.ZERO) THEN NCALC = 1 ELSE AYE = ONE + ALPHA TWOBYX = TWO/EX NCALC = 2 DO 90 I = 3,NB IF (TWOBYX.LT.ONE) THEN IF (ABS(BY(I-1))*TWOBYX.GE.XINF/AYE) GO TO 100 ELSE IF (ABS(BY(I-1)).GE.XINF/AYE/TWOBYX) GO TO 100 END IF BY(I) = TWOBYX*AYE*BY(I-1) - BY(I-2) AYE = AYE + ONE NCALC = NCALC + 1 90 CONTINUE END IF 100 DO 110 I = NCALC + 1,NB BY(I) = ZERO 110 CONTINUE ELSE BY(1) = ZERO NCALC = MIN(NB,0) - 1 END IF RETURN C---------- Last line of RYBESL ---------- END C END OF SUBROUTINE *** RYBESL *** C C BEGIN OF FUNCTION *** DLGAMA *** CS REAL FUNCTION ALGAMA(X) DOUBLE PRECISION FUNCTION DLGAMA(X) C C ################################################################## C C *** T. WIEDER, 05.03.1998 *** C C SOURCE FOR *** FUNCTION DLGAMA ***: C C SUBPROGRAM LIBRARY SPECFUN C C AVAILABILITIY: PUBLIC DOMAIN C C DEVELOPER: W.J. CODY AND L. STOLTZ, APPLED MATHEMATICS DIVISION, C ARGONNE NATIONAL LABORATORY, ARGONNE, IL 60439. C C C DISTRIBUTOR: NETLIB C C ################################################################## C C---------------------------------------------------------------------- C C This routine calculates the LOG(GAMMA) function for a positive real C argument X. Computation is based on an algorithm outlined in C references 1 and 2. The program uses rational functions that C theoretically approximate LOG(GAMMA) to at least 18 significant C decimal digits. The approximation for X > 12 is from reference C 3, while approximations for X < 12.0 are similar to those in C reference 1, but are unpublished. The accuracy achieved depends C on the arithmetic system, the compiler, the intrinsic functions, C and proper selection of the machine-dependent constants. C C C********************************************************************* C********************************************************************* C C Explanation of machine-dependent constants C C beta - radix for the floating-point representation C maxexp - the smallest positive power of beta that overflows C XBIG - largest argument for which LN(GAMMA(X)) is representable C in the machine, i.e., the solution to the equation C LN(GAMMA(XBIG)) = beta**maxexp C XINF - largest machine representable floating-point number; C approximately beta**maxexp. C EPS - The smallest positive floating-point number such that C 1.0+EPS .GT. 1.0 C FRTBIG - Rough estimate of the fourth root of XBIG C C C Approximate values for some important machines are: C C beta maxexp XBIG C C CRAY-1 (S.P.) 2 8191 9.62E+2461 C Cyber 180/855 C under NOS (S.P.) 2 1070 1.72E+319 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 2 128 4.08E+36 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 2 1024 2.55D+305 C IBM 3033 (D.P.) 16 63 4.29D+73 C VAX D-Format (D.P.) 2 127 2.05D+36 C VAX G-Format (D.P.) 2 1023 1.28D+305 C C C XINF EPS FRTBIG C C CRAY-1 (S.P.) 5.45E+2465 7.11E-15 3.13E+615 C Cyber 180/855 C under NOS (S.P.) 1.26E+322 3.55E-15 6.44E+79 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 3.40E+38 1.19E-7 1.42E+9 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 1.79D+308 2.22D-16 2.25D+76 C IBM 3033 (D.P.) 7.23D+75 2.22D-16 2.56D+18 C VAX D-Format (D.P.) 1.70D+38 1.39D-17 1.20D+9 C VAX G-Format (D.P.) 8.98D+307 1.11D-16 1.89D+76 C C************************************************************** C************************************************************** C C Error returns C C The program returns the value XINF for X .LE. 0.0 or when C overflow would occur. The computation is believed to C be free of underflow and overflow. C C C Intrinsic functions required are: C C LOG C C C References: C C 1) W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for C the Natural Logarithm of the Gamma Function,' Math. Comp. 21, C 1967, pp. 198-203. C C 2) K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May, C 1969. C C 3) Hart, Et. Al., Computer Approximations, Wiley and sons, New C York, 1968. C C C Authors: W. J. Cody and L. Stoltz C Argonne National Laboratory C C Latest modification: June 16, 1988 C C---------------------------------------------------------------------- CS REAL C .. Scalar Arguments .. DOUBLE PRECISION X C .. C .. Local Scalars .. DOUBLE PRECISION CORR,D1,D2,D4,EPS,FOUR,FRTBIG,HALF,ONE,PNT68,RES, + SQRTPI,THRHAL,TWELVE,TWO,XBIG,XDEN,XINF,XM1,XM2, + XM4,XNUM,Y,YSQ,ZERO INTEGER I C .. C .. Local Arrays .. DOUBLE PRECISION C(7),P1(8),P2(8),P4(8),Q1(8),Q2(8),Q4(8) C .. C .. Intrinsic Functions .. INTRINSIC LOG C .. C .. Data statements .. C---------------------------------------------------------------------- C Mathematical constants C---------------------------------------------------------------------- CS DATA ONE,HALF,TWELVE,ZERO/1.0E0,0.5E0,12.0E0,0.0E0/, CS 1 FOUR,THRHAL,TWO,PNT68/4.0E0,1.5E0,2.0E0,0.6796875E0/, CS 2 SQRTPI/0.9189385332046727417803297E0/ C---------------------------------------------------------------------- C Machine dependent parameters C---------------------------------------------------------------------- CS DATA XBIG,XINF,EPS,FRTBIG/4.08E36,3.401E38,1.19E-7,1.42E9/ C---------------------------------------------------------------------- C Numerator and denominator coefficients for rational minimax C approximation over (0.5,1.5). C---------------------------------------------------------------------- CS DATA D1/-5.772156649015328605195174E-1/ CS DATA P1/4.945235359296727046734888E0,2.018112620856775083915565E2, CS 1 2.290838373831346393026739E3,1.131967205903380828685045E4, CS 2 2.855724635671635335736389E4,3.848496228443793359990269E4, CS 3 2.637748787624195437963534E4,7.225813979700288197698961E3/ CS DATA Q1/6.748212550303777196073036E1,1.113332393857199323513008E3, CS 1 7.738757056935398733233834E3,2.763987074403340708898585E4, CS 2 5.499310206226157329794414E4,6.161122180066002127833352E4, CS 3 3.635127591501940507276287E4,8.785536302431013170870835E3/ C---------------------------------------------------------------------- C Numerator and denominator coefficients for rational minimax C Approximation over (1.5,4.0). C---------------------------------------------------------------------- CS DATA D2/4.227843350984671393993777E-1/ CS DATA P2/4.974607845568932035012064E0,5.424138599891070494101986E2, CS 1 1.550693864978364947665077E4,1.847932904445632425417223E5, CS 2 1.088204769468828767498470E6,3.338152967987029735917223E6, CS 3 5.106661678927352456275255E6,3.074109054850539556250927E6/ CS DATA Q2/1.830328399370592604055942E2,7.765049321445005871323047E3, CS 1 1.331903827966074194402448E5,1.136705821321969608938755E6, CS 2 5.267964117437946917577538E6,1.346701454311101692290052E7, CS 3 1.782736530353274213975932E7,9.533095591844353613395747E6/ C---------------------------------------------------------------------- C Numerator and denominator coefficients for rational minimax C Approximation over (4.0,12.0). C---------------------------------------------------------------------- CS DATA D4/1.791759469228055000094023E0/ CS DATA P4/1.474502166059939948905062E4,2.426813369486704502836312E6, CS 1 1.214755574045093227939592E8,2.663432449630976949898078E9, CS 2 2.940378956634553899906876E10,1.702665737765398868392998E11, CS 3 4.926125793377430887588120E11,5.606251856223951465078242E11/ CS DATA Q4/2.690530175870899333379843E3,6.393885654300092398984238E5, CS 2 4.135599930241388052042842E7,1.120872109616147941376570E9, CS 3 1.488613728678813811542398E10,1.016803586272438228077304E11, CS 4 3.417476345507377132798597E11,4.463158187419713286462081E11/ C---------------------------------------------------------------------- C Coefficients for minimax approximation over (12, INF). C---------------------------------------------------------------------- CS DATA C/-1.910444077728E-03,8.4171387781295E-04, CS 1 -5.952379913043012E-04,7.93650793500350248E-04, CS 2 -2.777777777777681622553E-03,8.333333333333333331554247E-02, CS 3 5.7083835261E-03/ DATA ONE,HALF,TWELVE,ZERO/1.0D0,0.5D0,12.0D0,0.0D0/,FOUR,THRHAL, + TWO,PNT68/4.0D0,1.5D0,2.0D0,0.6796875D0/, + SQRTPI/0.9189385332046727417803297D0/ DATA XBIG,XINF,EPS,FRTBIG/2.55D305,1.79D308,2.22D-16,2.25D76/ DATA D1/-5.772156649015328605195174D-1/ DATA P1/4.945235359296727046734888D0,2.018112620856775083915565D2, + 2.290838373831346393026739D3,1.131967205903380828685045D4, + 2.855724635671635335736389D4,3.848496228443793359990269D4, + 2.637748787624195437963534D4,7.225813979700288197698961D3/ DATA Q1/6.748212550303777196073036D1,1.113332393857199323513008D3, + 7.738757056935398733233834D3,2.763987074403340708898585D4, + 5.499310206226157329794414D4,6.161122180066002127833352D4, + 3.635127591501940507276287D4,8.785536302431013170870835D3/ DATA D2/4.227843350984671393993777D-1/ DATA P2/4.974607845568932035012064D0,5.424138599891070494101986D2, + 1.550693864978364947665077D4,1.847932904445632425417223D5, + 1.088204769468828767498470D6,3.338152967987029735917223D6, + 5.106661678927352456275255D6,3.074109054850539556250927D6/ DATA Q2/1.830328399370592604055942D2,7.765049321445005871323047D3, + 1.331903827966074194402448D5,1.136705821321969608938755D6, + 5.267964117437946917577538D6,1.346701454311101692290052D7, + 1.782736530353274213975932D7,9.533095591844353613395747D6/ DATA D4/1.791759469228055000094023D0/ DATA P4/1.474502166059939948905062D4,2.426813369486704502836312D6, + 1.214755574045093227939592D8,2.663432449630976949898078D9, + 2.940378956634553899906876D10,1.702665737765398868392998D11, + 4.926125793377430887588120D11,5.606251856223951465078242D11/ DATA Q4/2.690530175870899333379843D3,6.393885654300092398984238D5, + 4.135599930241388052042842D7,1.120872109616147941376570D9, + 1.488613728678813811542398D10,1.016803586272438228077304D11, + 3.417476345507377132798597D11,4.463158187419713286462081D11/ DATA C/-1.910444077728D-03,8.4171387781295D-04, + -5.952379913043012D-04,7.93650793500350248D-04, + -2.777777777777681622553D-03,8.333333333333333331554247D-02, + 5.7083835261D-03/ C .. C---------------------------------------------------------------------- Y = X IF ((Y.GT.ZERO) .AND. (Y.LE.XBIG)) THEN IF (Y.LE.EPS) THEN RES = -LOG(Y) ELSE IF (Y.LE.THRHAL) THEN C---------------------------------------------------------------------- C EPS .LT. X .LE. 1.5 C---------------------------------------------------------------------- IF (Y.LT.PNT68) THEN CORR = -LOG(Y) XM1 = Y ELSE CORR = ZERO XM1 = (Y-HALF) - HALF END IF IF ((Y.LE.HALF) .OR. (Y.GE.PNT68)) THEN XDEN = ONE XNUM = ZERO DO 10 I = 1,8 XNUM = XNUM*XM1 + P1(I) XDEN = XDEN*XM1 + Q1(I) 10 CONTINUE RES = CORR + (XM1* (D1+XM1* (XNUM/XDEN))) ELSE XM2 = (Y-HALF) - HALF XDEN = ONE XNUM = ZERO DO 20 I = 1,8 XNUM = XNUM*XM2 + P2(I) XDEN = XDEN*XM2 + Q2(I) 20 CONTINUE RES = CORR + XM2* (D2+XM2* (XNUM/XDEN)) END IF ELSE IF (Y.LE.FOUR) THEN C---------------------------------------------------------------------- C 1.5 .LT. X .LE. 4.0 C---------------------------------------------------------------------- XM2 = Y - TWO XDEN = ONE XNUM = ZERO DO 30 I = 1,8 XNUM = XNUM*XM2 + P2(I) XDEN = XDEN*XM2 + Q2(I) 30 CONTINUE RES = XM2* (D2+XM2* (XNUM/XDEN)) ELSE IF (Y.LE.TWELVE) THEN C---------------------------------------------------------------------- C 4.0 .LT. X .LE. 12.0 C---------------------------------------------------------------------- XM4 = Y - FOUR XDEN = -ONE XNUM = ZERO DO 40 I = 1,8 XNUM = XNUM*XM4 + P4(I) XDEN = XDEN*XM4 + Q4(I) 40 CONTINUE RES = D4 + XM4* (XNUM/XDEN) ELSE C---------------------------------------------------------------------- C Evaluate for argument .GE. 12.0, C---------------------------------------------------------------------- RES = ZERO IF (Y.LE.FRTBIG) THEN RES = C(7) YSQ = Y*Y DO 50 I = 1,6 RES = RES/YSQ + C(I) 50 CONTINUE END IF RES = RES/Y CORR = LOG(Y) RES = RES + SQRTPI - HALF*CORR RES = RES + Y* (CORR-ONE) END IF ELSE C---------------------------------------------------------------------- C Return for bad arguments C---------------------------------------------------------------------- RES = XINF END IF C---------------------------------------------------------------------- C Final adjustments and return C---------------------------------------------------------------------- CS ALGAMA = RES DLGAMA = RES RETURN C ---------- Last line of DLGAMA ---------- END C END OF FUNCTION *** DLGAMA *** C C BEGIN OF SUBROUTINE *** INTHP *** SUBROUTINE INTHP(A,B,D,F,M,P,EPS,INF,QUADR) C C ################################################################## C C *** T. WIEDER, 05.03.1998 *** C C SOURCE FOR SUBROUTINE *** INTHP ***: C C SUBPROGRAM LIBRARY TOMS C C AVAILABILITIY: PUBLIC DOMAIN C C DEVELOPER: K. SIKORSKY, F. STENGER, AND J. SCHWING. C C DISTRIBUTOR: TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 10 (1984), PP. 152 - 160. C C ################################################################## C C ALGORITHM 614 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.10, NO. 2, C JUN., 1984, P. 152-160. C C THIS SUBROUTINE COMPUTES INTEGRAL OF FUNCTIONS WHICH C MAY HAVE SINGULARITIES AT ONE OR BOTH END-POINTS OF AN C INTERVAL (A,B), SEE [1, 2]. IT CONTAINS FOUR DIFFERENT C QUADRATURE ROUTINES: ONE OVER A FINITE INTERVAL (A,B), C TWO OVER (A,+INFINITY), AND ONE OVER (-INFINITY,+INFINITY). C OF THE TWO FORMULAS OVER (A,+INFINITY), THE FIRST (INF=2 C BELOW) IS MORE SUITED TO NON-OSCILLATORY INTEGRANDS, WHILE C THE SECOND (INF=3) IS MORE SUITED TO OSCILLATORY INTEGRANDS. C THE USER SUPPLIES THE INTEGRAND FUNCTION, HE SPECIFIES THE C INTERVAL, AS WELL AS THE RELATIVE ERROR TO WHICH THE INTE- C GRAL IS TO BE EVALUATED. C THE FORMULAS ARE OPTIMAL IN CERTAIN HARDY SPACES H(P,DD), C SEE [1, 2]. HERE DD IS AN OPEN DOMAIN IN THE COMPLEX PLANE, C A AND B BELONG TO THE BOUNDARY OF DD AND H(P,DD), P.GT.1, IS C THE SET OF ALL ANALYTIC FUNCTONS IN DD WHOSE P-TH NORM DEFI- C NED AS IN [2] IS FINITE. C IF THE USER IS UNABLE TO SPECIFY THE PARAMETERS P AND D C OF THE SPACE H(P,DD) TO WHICH HIS INTEGRAND BELONGS, THE C ALGORITHM TERMINATES ACCORDING TO A HEURISTIC CRITERION, SEE C [2] AND COMMENTS TO EPS. C IF THE USER CAN SPECIFY THE PARAMETERS P AND D OF THE C SPACE H(P,DD) TO WHICH HIS INTEGRAND BELONGS, THE ALGORITHM C TERMINATES WITH AN ANSWER HAVING A GUARANTEED ACCURACY ( DE- C TEMINISTIC CRITERION, SEE [1, 2] AND COMMENTS TO EPS). C C C INPUT PARAMETERS C C C A = LOWER LIMIT OF INTEGRATION (SEE COMMENTS TO INF). C C B = UPPER LIMIT OF INTEGRATION (SEE COMMENTS TO INF). C C D = A PARAMETER OF THE CLASS H(P,DD) (SEE COMMENTS TO C INF). C C USER SETS D: C C HEURISTIC TERMINATION C = ANY REAL NUMBER C C DETERMINISTIC TERMINATION C = A NUMBER IN THE RANGE 0.LT.D.LE.PI/2. C C F = A NAME OF AN EXTERNAL INTEGRAND FUNCTION TO BE C SUPPLIED BY THE USER. F(X) COMPUTES THE VALUE OF C A FUNCTION F AT A POINT X. THE STATEMENT C ...EXTERNAL F... MUST APPEAR IN THE MAIN PROGRAM. C C M = MAXIMAL NUMBER OF FUNCTION EVALUATIONS ALLOWED IN C THE COMPUTATIONS, M.GE.3.( ALTERED ON EXIT ). C C P = 0, 1, .GT.1 A PARAMETER OF THE CLASS H(P,DD). C C USER SETS P: C = 0 - HEURISTIC TERMINATION. C = 1 - DETERMINISTIC TERMINATION WITH THE INFINITY C NORM. C .GT.1 -DETERMINISTIC TERMINATION WITH THE P-TH NORM. C C EPS = A REAL NUMBER - THE RELATIVE ERROR BOUND - SEE C REMARKS BELOW. ( ALTERED ON EXIT ). C C C INF = 1, 2, 3, 4 - INFORMATION PARAMETER. ( ALTERED ON EXIT ). C C = 1 SIGNIFIES AN INFINITE INTERVAL (A,B)=REAL LINE, C A AND B ANY NUMBERS. C DETERMINISTIC TERMINATION - C DD=STRIP(Z:ABS(IM(Z)).LT.D). C C = 2 SIGNIFIES A SEMI-INFINITE INTERVAL (A, +INFINITY) C USER SUPPLIES A, B ANY NUMBER. C QUADRATURE SUITED TO NON-OSCILLATORY INTEGRANDS. C DETERMINISTIC TERMINATION - C DD=SECTOR(Z:ABS(ARG(Z-A)).LT.D). C C = 3 SIGNIFIES A SEMI INFINITE INTERVAL (A,+INFINITY) C USER SUPPLIES A, B ANY NUMBER. C QUADRATURE SUITED TO OSCILLATORY INTEGRANDS. C DETERMINISTIC TERMINATION - C DD=REGION(Z:ABS(ARG(SINH(Z-A))).LT.D). C C = 4 SIGNIFIES A FINITE INTERVAL (A,B). C USER SUPPLIES A AND B. C DETERMINISTIC TERMINATION - C DD=LENS REGION(Z:ABS(ARG((Z-A)/(B-Z))).LT.D). C C C OUTPUT PARAMETERS C C C M = THE NUMBER OF FUNCTION EVALUATIONS USED IN THE C QUADRATURE. C C EPS = THE RELATIVE ERROR BOUND (SEE REMARKS BELOW). C C DETERMINISTIC TERMINATION C C = THE RELATIVE ERROR REXA BOUND, I.E., C REXA(F,M(OUTPUT)) .LE. EPS. C C HEURISTIC TERMINATION C C = MAX(EPS(INPUT),MACHEP). C C INF = 0, 1 - DETERMINISTIC TERMINATION C C = 0 COMPUTED QUADRATURE QCOM(F,M(EPS)), SEE REMARKS C BELOW. C C = 1 COMPUTED QUADRATURE QCOM(F,M1), SEE REMARKS C BELOW. C C INF = 2, 3, 4 - HEURISTIC TERMINATION. C C = 2 INTEGRATION COMPLETED WITH EPS=MAX(EPS(INPUT), C MACHEP). WE CAN EXPECT THE RELATIVE ERROR C REXA TO BE OF THE ORDER OF EPS (FOR SOME P.GE.1). C C = 3 INTEGRATION NOT COMPLETED. ATTEMPT TO EXCEED THE C MAXIMAL ALLOWED NUMBER OF FUNCTION EVALUATIONS M. C TRUNCATION CONDITIONS (SEE [2]) SATISFIED. QUADR C SET TO BE EQUAL TO THE LAST TRAPEZOIDAL APPRO- C XIMATION. IT IS LIKELY THAT QUADR APPROXIMATES THE C INTEGRAL IF M IS LARGE. C C = 4 INTEGRATION NOT COMPLETED. ATTEMPT TO EXCEED THE C MAXIMAL ALLOWED NUMBER OF FUNCTION EVALUATIONS M. C TRUNCATION CONDITIONS (SEE [2]) NOT SATISFIED. C QUADR SET TO BE EQUAL TO THE COMPUTED TRAPEZOIDAL C APPROXIMATION. IT IS UNLIKELY THAT QUADR APPROXIMATES C THE INTEGRAL. C C INF = 10, 11, 12, 13 - INCORRECT INPUT C C = 10 M.LT.3. C C = 11 P DOES NOT SATISFY P=0, P=1 OR P.GT.1 OR IN THE C CASE OF DETERMINISTIC TERMINATION D DOES NOT C SATISFY 0.LT.D.LE.PI/2. C C = 12 A.GE.B IN CASE OF A FINITE INTERVAL. C C = 13 INF NOT EQUAL TO 1, 2, 3, OR 4. C C C QUADR = THE COMPUTED VALUE OF QUADRATURE. C C C REMARKS: C C LET QEXA(F,M) ( QCOM(F,M) ) BE THE EXACT (COMPUTED) C VALUE OF THE QUADRATURE WITH M FUNCTION EVALUATIONS, C AND LET REXA(F,M) ( RCOM(F,M) ) BE THE RELATIVE ERROR C OF QEXA (QCOM) ,I.E., C REXA(F,M)=ABS(INTEGRAL(F)-QEXA(F,M))/NORM(F), C RCOM(F,M)=ABS(INTEGRAL(F)-QCOM(F,M))/NORM(F), C WITH THE NOTATION 0/0=0. C DUE TO THE ROUNDOFF ONE CANNOT EXPECT THE ERROR C RCOM TO BE LESS THAN THE RELATIVE MACHINE PRECISION C MACHEP. THEREFORE THE INPUT VALUE OF EPS IS CHANGED C ACCORDING TO THE FORMULA C EPS=MAX(EPS,MACHEP). C C DETERMINISTIC TERMINATON CASE C C THE NUMBER OF FUNCTON EVALUATIONS M(EPS) IS COMPUTED C SO THAT THE ERROR REXA IS NO GREATER THAN EPS,I.E., C C (*) REXA(F,M(EPS)) .LE. EPS . C C IF M(EPS).LE.M THEN THE QUADRATURE QCOM(F,M(EPS)) IS COM- C PUTED. OTHERWISE, WHICH MEANS THAT EPS IS TOO SMALL WITH C RESPECT TO M, THE QUADRATURE QCOM(F,M1) IS COMPUTED, WHERE C M1=2*INT((M-1)/2)+1. IN THIS CASE EPS IS CHANGED TO THE C SMALLEST NUMBER FOR WHICH THE ESTIMATE (*) HOLDS WITH C M(EPS)=M1 FUNCTION EVALUATIONS. C C HEURISTIC TERMINATION CASE C C WE CAN EXPECT THE RELATIVE ERROR REXA TO BE OF THE C ORDER OF EPS, SEE [2]. IF EPS IS TOO SMALL WITH RESPECT C TO M THEN THE QUADRATURE QCOM(F,M) IS COMPUTED. C C ROUNDOFF ERRORS C C IN BOTH DETERMINISTIC AND HEURISTIC CASES THE ROUND- C OFF ERROR C ROFF=ABS(QEXA(F,M)-QCOM(F,M)) C CAN BE ESTIMATED BY C C (**) ROFF .LE. 3*C1*R*MACHEP, C C WHERE R=QCOM(ABS(F),M)+(1+2*C2)/3*SUM(W(I),I=1,2,...M) C AND C1 IS OF THE ORDER OF UNITY, C1=1/(1-3*MACHEP), W(I) C ARE THE WEIGHTS OF THE QUADRATURE, SEE [2], AND C2 IS C A CONSTANT ESTIMATING THE ACCURACY OF COMPUTING FUNCTION C VALUES, I.E., C ABS(EXACT(F(X))-COMPUTED(F(X))).LE.C2*MACHEP. C IF THE INTEGRAND VALUES ARE COMPUTED INACCURATELY, I.E., C C2 IS LARGE, THEN THE ESTIMATE (**) IS LARGE AND ONE CAN C EXPECT THE ACTUAL ERROR ROFF TO BE LARGE. NUMERICAL TESTS C INDICATE THAT THIS HAPPENS ESPECIALLY WHEN THE INTEGRAND C IS EVALUATED INACCURATELY NEAR A SINGULARITY. THE WAYS OF C CIRCUMVENTING SUCH PITFALLS ARE EXPLAINED IN [2]. C C REFERENCES: C C [1] SIKORSKI,K., OPTIMAL QUADRATURE ALGORITHMS IN HP C SPACES, NUM. MATH., 39, 405-410 (1982). C [2] SIKORSKI,K., STENGER,F., OPTIMAL QUADRATURES IN C HP SPACES, ACM TOMS. C C C IMPLICIT NONE C C .. Scalar Arguments .. DOUBLE PRECISION A,B,D,EPS,P,QUADR INTEGER INF,M C .. C .. Function Arguments .. DOUBLE PRECISION F EXTERNAL F C .. C .. Local Scalars .. DOUBLE PRECISION ALFA,BA,C,C0,COR,E1,EPS3,EXPH,EXPH0,H,H0,H1,PI,S, + S1,SQ2,SR,SUM,SUM1,SUM2,T,U,V,V0,V1,V2,W,W1,W2, + W3,W4 INTEGER I,I1,K,L,L1,M1,M2,N,N1 LOGICAL INF1,INF2 C .. C .. External Functions .. DOUBLE PRECISION D1MACH EXTERNAL D1MACH C .. C .. Intrinsic Functions .. INTRINSIC ABS,DBLE,DLOG,DSQRT,DTAN,EXP,FLOAT,INT,SQRT C .. PI = 4.0D0*DTAN(1.0D0) C C CHECK THE INPUT DATA C IF (INF.NE.1 .AND. INF.NE.2 .AND. INF.NE.3 .AND. + INF.NE.4) GO TO 300 IF (M.LT.3) GO TO 270 IF (P.LT.1.0D0 .AND. P.NE.0.0D0) GO TO 280 IF (P.GE.1.0D0 .AND. (D.LE.0.D0.OR.D.GT.PI/2.D0)) GO TO 280 IF (INF.EQ.4 .AND. A.GE.B) GO TO 290 C SQ2 = DSQRT(2.0D0) I1 = INF - 2 BA = B - A N1 = 0 C C COMPUTE THE RELATIVE MACHINE PRECISION AND CHECK C THE VALUE OF EPS. CAUTION...THIS LOOP MAY NOT WORK ON A C MACHINE THAT HAS AN ACCURATED ARITHMETIC PROCESS COMPARED C TO THE STORAGE PRECISION. THE VALUE OF U MAY NEED TO BE C SIMPLY DEFINED AS THE RELATIVE ACCURACY OF STORAGE PRECISION. C U = 1.0D0 10 U = U/10.0D0 T = 1.0D0 + U IF (1.0D0.NE.T) GO TO 10 U = U*10.0D0 IF (EPS.LT.U) EPS = U C IF (P.EQ.0.D0) GO TO 40 C C SET UP DATA FOR THE DETERMINISTIC TERMINATION C IF (P.EQ.1.0D0) ALFA = 1.0D0 IF (P.GT.1.0D0) ALFA = (P-1.0D0)/P C = 2.0D0*PI/ (1.0D0-1.D0/EXP(PI*DSQRT(ALFA))) + 4.0D0**ALFA/ALFA W = DLOG(C/EPS) W1 = 1.0D0/ (PI*PI*ALFA)*W*W N = INT(W1) IF (W1.GT.FLOAT(N)) N = N + 1 IF (W1.EQ.0.D0) N = 1 N1 = 2*N + 1 SR = DSQRT(ALFA*DBLE(N)) IF (N1.LE.M) GO TO 20 C C EPS TOO SMALL WITH RESPECT TO M. COMPUTE THE NEW EPS C GUARANTEED BY THE VALUE OF M. C N1 = 1 N = INT(FLOAT((M-1)/2)) SR = DSQRT(ALFA*DBLE(N)) M = 2*N + 1 EPS = C/EXP(PI*SR) GO TO 30 C 20 M = N1 N1 = 0 30 H = 2.0D0*D/SR SUM2 = 0.0D0 L1 = N K = N INF1 = .FALSE. INF2 = .FALSE. H0 = H GO TO 50 C C SET UP DATA FOR THE HEURISTIC TERMINATION C 40 H = 1.0D0 H0 = 1.0D0 EPS3 = EPS/3.0D0 SR = SQRT(EPS) V1 = EPS*10.0D0 V2 = V1 M1 = M - 1 N = INT(FLOAT(M1/2)) M2 = N L1 = 0 INF1 = .TRUE. INF2 = .FALSE. C C INITIALIZE THE QUADRATURE C 50 I = 0 IF (INF.EQ.1) SUM = F(0.0D0) IF (INF.EQ.2) SUM = F(A+1.0D0) IF (INF.EQ.3) SUM = F(A+DLOG(1.0D0+SQ2))/SQ2 IF (INF.EQ.4) SUM = F((A+B)/2.0D0)/4.0D0*BA C C COMPUTE WEIGHTS, NODES AND FUNCTION VALUES C 60 EXPH = EXP(H) EXPH0 = EXP(H0) H1 = H0 E1 = EXPH0 U = 0.0D0 COR = 0.0D0 C C *** T. WIEDER, 24.02.1998 *** C70 IF (I1) 80,90,100 70 IF (I1.LT.0) THEN GO TO 80 ELSE IF (I1.EQ.0) THEN GO TO 90 ELSE IF (I1.GT.0) THEN GO TO 100 END IF C 80 V = F(H1) H1 = H1 + H GO TO 150 C 90 V = E1*F(A+E1) E1 = E1*EXPH GO TO 150 C 100 IF (INF.EQ.4) GO TO 140 W1 = DSQRT(E1+1.0D0/E1) W2 = DSQRT(E1) IF (E1.LT.0.1D0) GO TO 110 S = DLOG(E1+W1*W2) GO TO 130 110 W3 = E1 W4 = E1*E1 C0 = 1.0D0 S = E1 S1 = E1 T = 0.0D0 120 C0 = -C0* (0.5D0+T)* (2.0D0*T+1.D0)/ (2.0D0*T+3.0D0)/ (T+1.0D0) T = T + 1.0D0 W3 = W3*W4 S = S + C0*W3 IF (S.EQ.S1) GO TO 130 S1 = S GO TO 120 130 V = W2/W1*F(A+S) C *** 05.02.1998 *** IF (ABS(E1).GT.D1MACH(2)/10.0D0) THEN E1 = E1/10.0D0 E1 = E1*EXPH ELSE E1 = E1*EXPH END IF C *** 05.02.1998 *** GO TO 150 C 140 W1 = E1 + 1.0D0 V = E1/W1/W1*F((A+B*E1)/W1)*BA E1 = E1*EXPH C C SUMMATION ALGORITHM C 150 I = I + 1 SUM1 = U + V IF (ABS(U).LT.ABS(V)) GO TO 160 COR = V - (SUM1-U) + COR GO TO 170 160 COR = U - (SUM1-V) + COR 170 U = SUM1 IF (I.LT.L1) GO TO 70 C C SWITCH TO CHECK TRUNCATION CONDITION ( HEURISTIC C TERMINATION) C IF (INF1) GO TO 190 C C SWITCH TO COMPUTE THE MIDORDINATE APPROXIMATION C ( HEURISTIC TERMINATION ) OR TO STOP ( DETERMINIS- C TIC TERMINATION) C IF (INF2) GO TO 210 C C SET UP PARAMETERS TO CONTINUE SUMMATION C L1 = K 180 INF2 = .TRUE. I = 0 EXPH = 1.D0/EXPH H0 = -H0 E1 = 1.0D0/EXPH0 H1 = H0 H = -H GO TO 70 C C TRUNCATION CONDITION C 190 V0 = V1 V1 = V2 V2 = ABS(V) IF (V0+V1+V2.LE.EPS3) GO TO 200 IF (I.LT.M2) GO TO 70 N1 = 5 200 IF (INF2) K = I IF (.NOT.INF2) L = I V1 = 10.D0*EPS V2 = V1 M2 = M1 - L IF (.NOT.INF2) GO TO 180 C C N1=5 - TRUNCATION CONDITION NOT SATISFIED C IF (N1.EQ.5) GO TO 260 C C TRUNCATION CONDITION SATISFIED, SUM2=TRAPEZOIDAL C APPROXIMATION C SUM2 = SUM1 + COR + SUM M2 = 2* (K+L) C C CHECK THE NUMBER OF FUNCTION EVALUATIONS C IF (M2.GT.M1) GO TO 240 C C INITIALIZE ITERATION C INF1 = .FALSE. INF2 = .FALSE. L1 = L I = 0 H = -H H0 = H/2.D0 GO TO 60 C C P.GE.1 = DETERMINISTIC TERMINATION C 210 IF (P.GE.1.D0) GO TO 220 C C COMPUTE THE MIDORDINATE APPROXIMATION SUM1 C H = -H SUM1 = (SUM1+COR)*H W1 = (SUM1+SUM2)/2.0D0 C C TERMINATION CONDITION C IF (ABS(SUM1-SUM2).LE.SR) GO TO 230 C C SET UP DATA FOR THE NEXT ITERATION C M2 = 2*M2 IF (M2.GT.M1) GO TO 250 I = 0 K = 2*K L = 2*L L1 = L H = H/2.0D0 H0 = H/2.0D0 SUM2 = W1 INF2 = .FALSE. GO TO 60 C C FINAL RESULTS C 220 QUADR = -H* (SUM1+COR+SUM) INF = N1 RETURN C 230 QUADR = W1 INF = 2 M = M2 + 1 RETURN C 240 QUADR = SUM2 INF = 3 M = K + L + 1 RETURN C 250 QUADR = W1 INF = 3 M = M2/2 + 1 RETURN C 260 QUADR = U + COR + SUM INF = 4 M = K + L + 1 RETURN C 270 INF = 10 RETURN C 280 INF = 11 RETURN C 290 INF = 12 RETURN C 300 INF = 13 RETURN C END C END OF SUBROUTINE *** INTHP *** C C ################################################################## C ################################################################## C C BEGIN OF SUBROUTINE *** OSCINT *** SUBROUTINE OSCINT(AZERO,PERIOD,RFIRST,EPS,NQUAD,NDIM1,NDIM2,GAUSS, + HFUN,GPER,WORK,SAVPER,WEIGHT,ABSCIS,QLIST, + RESULT,ISTATE) C C ################################################################## C C *** T. WIEDER, 05.03.1998 *** C C SOURCE FOR SUBROUTINE *** OSCINT ***: C C SUBPROGRAM LIBRARY TOMS C C AVAILABILITIY: PUBLIC DOMAIN C C DEVELOPER: JAMES LYNESS AND GWENDOLEN HINES. C C DISTRIBUTOR: TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 12 (1986), PP. 24 - 25. C C ################################################################## C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C THIS ROUTINE PROVIDES AN APPROXIMATION **RESULT** TO THE C INTEGRAL OF AN OVERALL INTEGRAND **HFUN(X)*GPER(X)** OVER A C SEMI-INFINITE INTERVAL (AZERO,INFINITY). THE OVERALL INTEGRAND C FUNCTION SHOULD BE ULTIMATELY OSCILLATING IN SIGN AND PERIODIC C WITH USER-PROVIDED **PERIOD**. C C THE ROUTINE CONSIDERS A SEQUENCE OF INTERVALS. ALL THE C THE INTERVALS, EXCEPT THE FIRST, ARE OF LENGTH 0.5*PERIOD. C IT USES A QUADRATURE RULE TO APPROXIMATE THE INTEGRAL OVER C EACH INTERVAL. THESE APPROXIMATIONS ARE STORED IN C QLIST(J) J = 1,2,3,... . THE VALUES OF QLIST ULTIMATELY C OSCILLATE IN SIGN. THEN THE ROUTINE USES A SERIES ACCELERATION C TECHNIQUE BASED ON THE EULER TRANSFORMATION TO SUM THIS SEQUENCE. C C C INPUT PARAMETERS: C C AZERO LOWER LIMIT OF INTEGRATION C C PERIOD LEAST POSITIVE PERIOD OR ULTIMATE PERIOD OF C INTEGRAND FUNCTION C C RFIRST RIGHT-HAND ENDPOINT OF FIRST INTERVAL. SUGGESTED VALUE C FOR STRAIGHTFORWARD PROBLEMS IS AZERO. C SEE NOTE 1 BELOW C C EPS THE REQUESTED ACCURACY C C NQUAD INTEGER. NQ = ABS(NQUAD) IS THE NUMBER OFABSCISSAS TO C BE USED BY THE QUADRATURE RULE IN EACH INTERVAL. C NQUAD > 1 TRAPEZOIDAL RULE IS USED C NQUAD < 0 RULE SPECIFIED BY GAUSS BELOW IS USED C C NDIM1,NDIM2 DIMENSIONS OF THE OUTPUT ARRAY **WORK**. C NDIM1 = 10 IMPLIES NMAX = 100 C NDIM1 .NE. 10 IMPLIES NMAX = MIN(100,NDIM1) C NMAX IS A PHYSICAL LIMIT. IF CALCULATION IS NOT C COMPLETE AFTER NMAX INTERVALS HAVE BEEN CONSIDERED C IT IS THEN ABANDONED. (SEE NOTE 5 BELOW) C NDIM2 SUGGESTED VALUE IS 15. NORMALLY SHOULD EXCEED 4 C C GAUSS NAME OF A SUBROUTINE WHICH PROVIDES WEIGHTS C AND ABSCISSAS. WHEN NQUAD > 0 THIS IS NOT CALLED. C SEE NOTE 2 BELOW FOR USE WHEN NQUAD < 0. C C HFUN NAME OF A FUNCTION SUBROUTINE. THIS IS ONE FACTOR OF C THE INTEGRAND FUNCTION. THE OTHER FACTOR IS GPER BELOW C C GPER NAME OF A FUNCTION SUBROUTINE. THIS IS THE C CO-FACTOR OF HFUN IN THE INTEGRAND FUNCTION. C GPER MUST BE PERIODIC WITH PERIOD COINCIDING WITH THE C INPUT PARAMETER, PERIOD ABOVE. (SEE NOTE 3 BELOW) C C C OUTPUT PARAMETERS: C C WORK(NDIM1,NDIM2) HOLDS COMPLETED FINITE AVERAGE TABLE C C SAVPER( ) ARRAY TO SAVE FUNCTION VALUES OF GPER. C DIMENSION NOT LESS THAN 2*IABS(NQUAD) C C C WEIGHT( ) ARRAY TO HOLD WEIGHTS FOR GAUSSIAN QUADRATURE C DIMENSION NOT LESS THAN MAX(1,-NQUAD) C C ABSCIS( ) ARRAY TO HOLD ABSCISSAS FOR GAUSSIAN QUADRATURE C DIMENSION NOT LESS THAN MAX(1,-NQUAD) C C QLIST(100) ARRAY TO HOLD THE RESULT OF EACH QUADRATURE C C RESULT OVERALL RESULT OF THE INTEGRATION C C ISTATE(6) VECTOR OF INTEGERS GIVES STATUS OF RESULT C C ISTATE(1) AN INDICATOR WARNING ABOUT SUSPICIOUS RESULTS - C ZERO: THE RUN WAS APPARENTLY SUCCESSFUL. C POSITIVE: THE RUN WAS APPARENTLY SUCCESSFUL , BUT C THERE ARE UNSATISFACTORY FEATURES OF C POSSIBLE INTEREST TO THE SOPHISTICATED C USER C NEGATIVE: UNSUCCESSFUL RUN - DISREGARD THE RESULT C SEE NOTES BELOW FOR COMPLETE SPECIFICATION C ISTATE(2) LSIGCH - INDICATES LAST INTERVAL IN WHICH THE SIGN C OF THE INTEGRAL COINCIDED WITH THAT OF C THE INTEGRAL OVER THE PREVIOUS INTERVAL C SEE NOTE 4 BELOW. C ISTATE(3) NOW - THE NUMBER OF FINITE INTEGRALS, QLIST(Q), C EVALUATED IN THE CALCULATION. C ISTATE(4) NCOL - THE COLUMN OF THE FINITE AVERAGE TABLE, C (WORK), ON WHICH THE RESULT IS BASED. C ISTATE(5) NROW - THE ROW OF THE FINITE AVERAGE TABLE, C (WORK), ON WHICH THE RESULT IS BASED. C ISTATE(6) NCOUNT - THE NUMBER OF CALLS TO FUNCTION HFUN. C C C NOTE 1 C C RFIRST C THIS ALLOWS THE USER TO LOCATE HIS SUBDIVISION. FOR CAUTIOUS C RUNNING, ARRANGE RFIRST TO COINCIDE WITH AN "ULTIMATE ZERO". FOR C SLIGHTLY MORE ADVENTUROUS BUT LESS RELIABLE RUNNING, ARRANGE C RFIRST TO COINCIDE WITH AN "ULTIMATE" PEAK. OTHERWISE, SET C RFIRST < AZERO, IN WHICH CASE, OSCINT USES AZERO INSTEAD OF RFIRST. C C C NOTE 2 C C GAUSS C THIS IS THE NAME OF A USER-PROVIDED SUBROUTINE OF THE FORM C GAUSS(ITYPE,A,B,C,D,N,WEIGHT(N),ABS(N)CIS,IFAIL) C IT IS CALLED ONLY WHEN NQUAD IS NEGATIVE WITH N = -NQUAD. IT C RETURNS A SET OF WEIGHTS AND ABSCISSAS, SUITABLE FOR INTEGRATION C OVER THE INTERVAL (-1,1). C THE FIRST FIVE INPUT PARAMETERS OF GAUSS SHOULD BE IGNORED. C IFAIL MAY BE USED FOR A WARNING MESSAGE. IF GAUSS RETURNS C IFAIL .NE. 0, OSCINT ABORTS, SETTING ISTATE(1) = -4000. C THE USER WHO HAS THE NAG LIBRARY AVAILABLE MAY SET GAUSS = D01BCF. C C C NOTE 3 C C HFUN AND GPER C ONE MAY ALWAYS USE HFUN FOR THE INTEGRAND FUNCTION, F(X), AND C CODE GPER TO RETURN THE VALUE 1.0D0. HOWEVER, WHEN F(X) HAS A C PERIODIC FACTOR, J(X), REPETITIVE EVALUATION OF J(X) IN EACH C INTERVAL MAY BE AVOIDED BY SETTING GPER = J(X) AND HFUN = H(X). C THE ROUTINE CHECKS THAT GPER IS INDEED PERIODIC BY MAKING SOME C EVALUATIONS IN THE THIRD, FOURTH, AND FIFTH INTERVALS. IF IT C FINDS GPER IS NOT PERIODIC, IT TERMINATES WITH ISTATE(1) = -5000. C IF PERIOD <= 10**-5 , THE ROUTINE TERMINATES WITH ISTATE(1) = -3000. C C C NOTE 4 C C TERMINATION C IN THIS NOTE, THE "NORMAL SIGN CHANGE PATTERN" IS ONE IN WHICH C SUCCESSIVE VALUES OF QLIST(Q) OSCILLATE IN SIGN. THE "GRACE C PERIOD" IS Q<=10 WHEN THE NORMAL SIGN CHANGE PATTERN IS NOT C INSISTED ON. LSIGCH IS THE HIGHEST VALUE OF Q FOR WHICH C QLIST(Q)*QLIST(Q-1) IS POSITIVE. TERMINATION COMES ABOUT C AFTER CALCULATING QLIST(NOW), WHEN EITHER C C (A) AN APPROXIMATION OF SUFFICIENT ACCURACY IS CURRENTLY AVAIL- C ABLE (THE ROUTINE SETS ISTATE(1) = MAX(0,4-(NOW-LSIGCH)) C OR C C (B) THE NORMAL SIGN CHANGE PATTERN IS VIOLATED AFTER THE GRACE C PERIOD. (THE ROUTINE SETS ISTATE(1) = -200) C OR C C (C) THE USER SET LIMIT, NMAX, OF INTERVALS HAVE BEEN CALCULATED C I.E. NOW = NMAX. (THE ROUTINE SETS ISTATE(1) = -100) C C THE ROUTINE CHECKS (A), THEN (B), AND THEN (C). TERMINATION C UNDER (A) MAY OCCUR WITHOUT ANY NORMAL SIGN CHANGE PATTERN C EMERGING. THIS MAY BE DUE TO MISUSE OF THE ROUTINE, BUT THE C PROBLEM IS SMALL ENOUGH TO BE CORRECTLY HANDLED. IN THIS CASE, C ISTATE(1) = MAX(0,4-(NOW-LSIGCH)) MAY BE A SMALL POSITIVE INTEGER. C C C NOTE 5 C C VARIABLE DIMENSIONS AND STORAGE ECONOMY C THE PROGRAM WHICH CALLS OSCINT HAS TO PROVIDE NUMERICAL VALUES C OF THE INPUT PARAMETERS NQUAD, NDIM1 AND NDIM2. IT MUST ALSO C INCLUDE A DIMENSION STATEMENT IN WHICH THE FIRST FIVE AND THE C LAST OUTPUT PARAMETERS OF OSCINT ARE DIMENSIONED. NOTE THAT THE C DIMENSIONS OF WEIGHT AND ABSCIS ARE BOTH AT LEAST 1 WHEN NQUAD C IS POSITIVE AND AT LEAST -NQUAD OTHERWISE. HOWEVER, THE C DIMENSION OF SAVPER IS AT LEAST 2*ABS(NQUAD) REGARDLESS OF THE C SIGN OF NQUAD. C GENERALLY, THESE VARIABLE DIMENSION STATEMENTS ALLOW ECONOMIC C USE OF STORAGE. THERE IS ONE FURTHER STORAGE SAVING FEATURE. C WHEN NDIM1 = 10, THE ROUTINE CARRIES OUT THE SAME CALCULATION C AS IT WOULD IF NDIM1 = 100. HOWEVER, INSTEAD OF USING A WORK C ARRAY OF DIMENSION (100,NDIM2), IT USES A WORK ARRAY OF DIMENSION C (10,NDIM2) AND OVERWRITES IT, AS AND WHEN NECESSARY, AS THE C CALCULATION PROCEEDS. IN THIS CASE, OSCINT BEHAVES AS IF NDIM2 WERE C REPLACED BY MIN(20,NDIM2). C SOME OBVIOUS ERRORS IN DIMENSIONING, SUCH AS NEGATIVE C DIMENSIONS, CAUSE THE ROUTINE TO TERMINATE WITH C ISTATE(1) = -6000 C HOWEVER, INADEQUATE DIMENSIONING IN THE CALLING PROGRAM MAY GO C UNDETECTED AND MAY LEAD TO RANDOM OR CHAOTIC OUTPUT. C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C IMPLICIT NONE C *** T. WIEDER, 15.04.1998 *** C DOUBLE PRECISION SAVPER(*),WEIGHT(*),ABSCIS(*) C C .. Parameters .. INTEGER NDIMMA PARAMETER (NDIMMA=1000) C .. C .. Scalar Arguments .. DOUBLE PRECISION AZERO,EPS,PERIOD,RESULT,RFIRST INTEGER NDIM1,NDIM2,NQUAD C .. C .. Array Arguments .. DOUBLE PRECISION ABSCIS(36),QLIST(0:NDIMMA),SAVPER(0:256), + WEIGHT(36),WORK(0:NDIM1,0:NDIM2) INTEGER ISTATE(6) C .. C .. Function Arguments .. DOUBLE PRECISION GPER,HFUN EXTERNAL GPER,HFUN C .. C .. Subroutine Arguments .. EXTERNAL GAUSS C .. C .. Local Scalars .. DOUBLE PRECISION CURR,HASPER,PREV,S,WMIN INTEGER I,II,ISUB,J,JP,K,MJ,MJ1,MJ2,MJK,MJM2,NMAX,NML,NOW,NOWJP, + NROUND,NROW C .. C .. External Functions .. DOUBLE PRECISION QRULE EXTERNAL QRULE C .. C .. Intrinsic Functions .. INTRINSIC ABS,DABS,MAX,MIN,MOD C .. NROUND = 0 HASPER = 0.5D0*PERIOD WMIN = 1.0D0 DO 10 I = 1,6 ISTATE(I) = 0 10 CONTINUE IF (NQUAD.EQ.0 .OR. NQUAD.EQ.1 .OR. NDIM1.LT.1 .OR. + NDIM2.LT.1) ISTATE(1) = -6000 IF (PERIOD.LT.10.0D-5) ISTATE(1) = -5000 IF (ISTATE(1).LT.0) GO TO 120 C C *** T. WIEDER, 18.02.1998 *** IF (NDIM1.EQ.10) THEN NMAX = 100 ELSE C NMAX = MIN(100,NDIM1) NMAX = MIN(NDIMMA,NDIM1) END IF C JP = 0 S = 0.0D0 C C LOOP TO CONSTRUCT TABLE C 20 IF (JP.EQ.NMAX) ISTATE(1) = -100 IF (ISTATE(1).NE.0) GO TO 40 C K = 0 JP = JP + 1 J = MIN(JP,NDIM2) IF (NDIM1.EQ.10) J = MIN(J,20) 30 IF (K.LT.J) THEN K = K + 1 IF (NDIM1.EQ.10) THEN MJ = MOD(JP,10) MJM2 = MOD(JP-2,10) MJ1 = MOD(JP-K+1,10) MJ2 = MOD(JP-K+2,10) MJK = MOD(JP-K,10) IF (MJ.EQ.0) MJ = 10 IF (MJM2.EQ.0) MJM2 = 10 IF (MJ1.EQ.0) MJ1 = 10 IF (MJ2.EQ.0) MJ2 = 10 IF (MJK.EQ.0) MJK = 10 C ELSE MJ = JP MJM2 = JP - 2 MJ1 = JP - K + 1 MJ2 = JP - K + 2 MJK = JP - K END IF C IF (K.EQ.1) THEN C C *** AMENDMENT SUGGESTED BY DR. T.R. HOPKINS, 15.04.1998 *** C WORK(MJ,K)=QRULE(JP-1,AZERO,HASPER,RFIRST,NQUAD,NDIM1,NDIM2,GAUSS C 1HFUN,GPER,SAVPER,WEIGHT,ABSCIS,QLIST,ISTATE) WORK(MJ,K) = QRULE(JP-1,AZERO,HASPER,RFIRST,NQUAD,GAUSS, + HFUN,GPER,SAVPER,WEIGHT,ABSCIS,QLIST,ISTATE) C ELSE WORK(MJ1,K) = (WORK(MJ1,K-1)+WORK(MJ2,K-1))/2.0D0 IF (DABS(WORK(MJ1,K)).LT.WMIN) THEN WMIN = DABS(WORK(MJ1,K)) NOW = K NROW = MJ1 NOWJP = JP END IF C CURR = ABS(WORK(MJ1,K)) PREV = ABS(WORK(MJK,K)) IF (JP.NE.K) THEN IF ((CURR.LT.EPS) .AND. (PREV.LT.EPS)) THEN NOW = K NROW = MJ1 NOWJP = JP GO TO 40 C END IF C END IF C END IF C GO TO 30 C END IF C GO TO 20 C 40 DO 50 I = 1,NOWJP - NOW S = S + QLIST(I) 50 CONTINUE IF (NDIM1.NE.10) GO TO 100 NROUND = NOW - 10 IF (NROUND.LE.0) NROUND = 0 IF (NROUND.LE.0) GO TO 100 C DO 60 I = NOWJP - NOW + 1,NOWJP - 10 S = S + QLIST(I) 60 CONTINUE ISUB = NOWJP - 9 70 IF (ISUB.GT.10) ISUB = ISUB - 10 IF (ISUB.GT.10) GO TO 70 DO 80 J = 1,NROUND S = S + WORK(ISUB,J)/2.0D0 80 CONTINUE DO 90 I = NROW,NROW + NROUND - 1 II = I IF (I.GT.10) II = I - 10 S = S - WORK(II,NROUND+1) 90 CONTINUE C 100 DO 110 J = NROUND + 1,NOW - 1 S = S + WORK(NROW,J)/2.0D0 110 CONTINUE S = S + WORK(NROW,NOW) RESULT = S ISTATE(3) = JP ISTATE(4) = NOW ISTATE(5) = NROW NML = ISTATE(3) - ISTATE(2) IF (NML.LT.4 .AND. ISTATE(1).EQ.0) ISTATE(1) = MAX(0,4-NML) 120 RETURN END C END OF SUBROUTINE *** OSCINT *** C C BEGIN OF FUNCTION *** QRULE *** C C *** AMENDMENT SUGGESTED BY DR. T.R. HOPKINS, 15.04.1998 *** C DOUBLE PRECISION FUNCTION QRULE (J,AZERO,HASPER,RFIRST,NQUAD, C 1NDIM1,NDIM2,GAUSS,HFUN,GPER,SAVPER,WEIGHT,ABSCIS,QLIST,ISTATE) DOUBLE PRECISION FUNCTION QRULE(J,AZERO,HASPER,RFIRST,NQUAD,GAUSS, + HFUN,GPER,SAVPER,WEIGHT,ABSCIS, + QLIST,ISTATE) C C ################################################################## C C *** T. WIEDER, 05.03.1998 *** C C SOURCE FOR *** FUNCTION QRULE ***: C C SUBPROGRAM LIBRARY TOMS C C AVAILABILITIY: PUBLIC DOMAIN C C DEVELOPER: JAMES LYNESS AND GWENDOLEN HINES. C C DISTRIBUTOR: TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 12 (1986), PP. 24 - 25. C C ################################################################## C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C THIS ROUTINE EVALUATES THE INTEGRAL OF THE FUNCTION, C HFUN(X)*GPER(X) OVER THE INTERVAL (A,B) WHERE: C C GENERALLY (J>0) A = PREVIOUS VALUE OF B C B = A + HASPER C (J=0) A = AZERO C B = RFIRST WHEN RFIRST > AZERO C B = AZERO + HASPER WHEN RFIRST <= AZERO C C C INPUT PARAMETERS: C C J DEFINES WHICH TERM, QLIST(J), OF C THE SERIES IS BEING EVALUATED C C HASPER HALF THE PERIOD C C C OTHER INPUT AND OUTPUT PARAMETERS ARE IDENTICAL TO THOSE IN C OSCINT (DESCRIBED ABOVE). C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C *** T. WIEDER, 18.02.1998 *** C IMPLICIT NONE C ORIGINAL VALUE FOR *** IGRACE *** = 9 C ORIGINAL VALUE FOR *** NDIMMA *** = 100 C *** T. WIEDER, 15.04.1998 *** C DOUBLE PRECISION SAVPER(*),WEIGHT(*),ABSCIS(*) C C .. Parameters .. INTEGER IGRACE,NDIMMA PARAMETER (IGRACE=99,NDIMMA=1000) C .. C .. Scalar Arguments .. DOUBLE PRECISION AZERO,HASPER,RFIRST INTEGER J,NQUAD C .. C .. Array Arguments .. DOUBLE PRECISION ABSCIS(36),QLIST(0:NDIMMA),SAVPER(0:256), + WEIGHT(36) INTEGER ISTATE(6) C .. C .. Function Arguments .. DOUBLE PRECISION GPER,HFUN EXTERNAL GPER,HFUN C .. C .. Subroutine Arguments .. EXTERNAL GAUSS C .. C .. Local Scalars .. DOUBLE PRECISION A,AA,B,BB,CC,DD,DIFF,FUN,GMAX,TENDPT,WSUM,WT,XI,Y INTEGER I,IELM,IFAIL,INDEX,ITYPE,NPTS C .. C .. Intrinsic Functions .. INTRINSIC ABS,DABS,MOD C .. C .. Save statement .. SAVE B,TENDPT C .. NPTS = ABS(NQUAD) C IF (J.EQ.0) THEN INDEX = 0 DO 10 I = 1,NDIMMA QLIST(I) = 0.0D0 10 CONTINUE C CALL GAUSS IF NQUAD < 0 IF (NQUAD.LT.0) THEN ITYPE = 0 AA = -1.0D0 BB = 1.0D0 CC = 0.0D0 DD = 0.0D0 IFAIL = 0 CALL GAUSS(ITYPE,AA,BB,CC,DD,NPTS,WEIGHT,ABSCIS,IFAIL) IF (IFAIL.NE.0) THEN ISTATE(1) = -4000 GO TO 40 C END IF C END IF C END IF C QRULE = 0.0D0 C C SET INTERVAL ENDPOINTS, A AND B IF (J.NE.0) THEN A = B C ELSE A = AZERO END IF C IF (RFIRST.LE.AZERO .OR. J.NE.0) THEN B = A + HASPER C ELSE B = RFIRST END IF C C LOOP ON ABSCISSAS FOR INTERVAL #J - STARTS TO CALCULATE QRULE DO 30 I = 1,NPTS XI = I C C CALCULATE ABSCISSA Y IF (NQUAD.LT.0) THEN Y = (B-A)*ABSCIS(I)/2.0D0 + (B+A)/2.0D0 WT = WEIGHT(I) C ELSE Y = ((XI-1)*B+A* (NPTS-I))/ (NPTS-1) WT = 1.0D0 END IF C C IELM IS LOCATION IN SAVPER FOR GPER FUNCTION VALUES C J=0: IGNORES. J=1,2: WHERE TO PUT VALUE. J>2: WHERE TO GET C VALUE FROM IELM = NPTS*MOD(J-1,2) + I C IF (J.GT.0 .AND. I.EQ.1 .AND. NQUAD.GT.0) GO TO 20 IF (J.EQ.0) THEN FUN = HFUN(Y)*GPER(Y) C ELSE IF (J.LT.3) THEN SAVPER(IELM) = GPER(Y) C CHECK FOR CONSTANT FUNCTION IF (SAVPER(IELM).EQ.SAVPER(1)) INDEX = INDEX + 1 C CODE MODIFIED BY DR. T.R. HOPKINS, 28.04.1998: IF (IELM.EQ.1) THEN GMAX = DABS(SAVPER(1)) ELSE IF (DABS(SAVPER(IELM)).GT.DABS(SAVPER(IELM-1))) THEN GMAX = DABS(SAVPER(IELM)) END IF END IF C FUN = HFUN(Y)*SAVPER(IELM) C ELSE IF (J.LT.5) THEN C CHECK THAT GPER IS PERIODIC WITH PERIOD=PERIOD DIFF = DABS(SAVPER(IELM)-GPER(Y)) C C *** T. WIEDER, 19.02.1998 *** C IF (DIFF.LT.GMAX*1.0D-5) THEN IF (DIFF.LT.GMAX*1.0D1) THEN FUN = HFUN(Y)*SAVPER(IELM) C ELSE ISTATE(1) = -3000 END IF C ELSE IF (INDEX.EQ.2*NPTS) THEN FUN = HFUN(Y) C ELSE FUN = HFUN(Y)*SAVPER(IELM) END IF C END IF C ISTATE(6) = ISTATE(6) + 1 C 20 IF (NQUAD.GT.0) THEN IF (I.EQ.1) THEN IF (J.EQ.0) FUN = FUN/2.0D0 IF (J.GT.0) FUN = TENDPT END IF C IF (I.EQ.NPTS) THEN FUN = FUN/2.0D0 TENDPT = FUN END IF C QRULE = QRULE + FUN C ELSE QRULE = QRULE + WT*FUN END IF C 30 CONTINUE C LOOP ON ABSCISSA FOR INTERVAL #J ENDS IF (INDEX.EQ.2*NPTS .AND. SAVPER(1).NE.1) THEN QRULE = QRULE*SAVPER(1) END IF C IF (NQUAD.GT.0) WSUM = NPTS - 1 IF (NQUAD.LT.0) WSUM = 2.0D0 QRULE = QRULE* (B-A)/WSUM QLIST(J+1) = QRULE C IF (QRULE*QLIST(J).GT.0) ISTATE(2) = J C C *** T. WIEDER, 18.02.1998 *** C IF (ISTATE(2).GT.9) ISTATE(1) = -200 IF (ISTATE(2).GT.IGRACE) ISTATE(1) = -200 C 40 RETURN END C END OF SUBROUTINE *** QRULE *** C C BEGIN OF SUBROUTINE *** G5AND9 *** SUBROUTINE G5AND9(ITYPE,A,B,C,D,NQUAD,WEIGHT,ABSCIS,IERR) C C ################################################################## C C *** T. WIEDER, 05.03.1998 *** C C SOURCE FOR SUBROUTINE *** G5AND9 ***: C C SUBPROGRAM LIBRARY TOMS C C AVAILABILITIY: PUBLIC DOMAIN C C DEVELOPER: JAMES LYNESS AND GWENDOLEN HINES. C C DISTRIBUTOR: TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 12 (1986), PP. 24 - 25. C C ################################################################## C C THIS IS AN EXTRACT FROM A QUADRATURE ROUTINE CONSTRUCTED ONLY FOR C USE IN A DRIVER WHICH ILLUSTRATES OSCINT. A NAG LIBRARY C SUBSCRIBER MAY REPLACE THIS BY D01BCF FOR GENERAL USE. C C IMPLICIT NONE C C C .. Scalar Arguments .. DOUBLE PRECISION A,B,C,D INTEGER IERR,ITYPE,NQUAD C .. C .. Array Arguments .. DOUBLE PRECISION ABSCIS(NQUAD),WEIGHT(NQUAD) C .. C .. Local Scalars .. INTEGER I C .. C .. Local Arrays .. DOUBLE PRECISION ABSC5(5),ABSC9(9),WT5(5),WT9(9) C .. C .. Data statements .. C THE FOLLOWING DATA ARE WEIGHTS AND ABSCISSAS OF THE FIVE POINT C AND THE NINE POINT GAUSS-LEGENDRE QUADRATURE RULES RESPECTIVELY. C NORMALISED TO THE INTERVAL (-1,1). DATA WT5/.2369268850561891D0,.4786286704993665D0, + .5688888888888889D0,.4786286704993665D0,.2369268850561891D0/ DATA ABSC5/-.9061798459386640D0,-.5384693101056831D0,0.0D0, + .5384693101056830D0,.9061798459386640D0/ DATA WT9/.0812743883615745D0,.1806481606948574D0, + .2606106964029355D0,.3123470770400029D0,.3302393550012598D0, + .3123470770400029D0,.2606106964029356D0,.1806481606948575D0, + .8127438836157467D-01/ DATA ABSC9/-.9681602395076261D0,-.8360311073266358D0, + -.6133714327005904D0,-.3242534234038090D0,0.0D0, + .3242534234038087D0,.6133714327005902D0,.8360311073266357D0, + .9681602395076260D0/ C .. C C *** T. WIEDER, 16.04.1998 *** C THIS IS A DIRTY TRICK TO PREVENT COMPILER FROM WARNING MESSAGE: A = A B = B C = C D = D ITYPE = ITYPE C IERR = 59 C IF (NQUAD.EQ.5) THEN IERR = 0 DO 10 I = 1,5 ABSCIS(I) = ABSC5(I) WEIGHT(I) = WT5(I) 10 CONTINUE END IF C IF (NQUAD.EQ.9) THEN IERR = 0 DO 20 I = 1,9 ABSCIS(I) = ABSC9(I) WEIGHT(I) = WT9(I) 20 CONTINUE END IF C RETURN END C END OF SUBROUTINE *** G5AND9 *** C BEGIN OF SUBROUTINE *** TBESSE *** SUBROUTINE TBESSE C C ################################################################# C C FIRST VERSION: 06.02.1998 C PREVIOUS VERSION: 06.02.1998 C LATEST VERSION: 05.03.1998 C C PROGRAMMED BY: C THOMAS WIEDER C E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE C C PURPOSE: C C TEST SUBROUTINE *** RJBESL ***. C C ################################################################# C C IMPLICIT NONE C C .. Local Scalars .. DOUBLE PRECISION ALPHA,NU,X INTEGER N,NB,NCALC C .. C .. Local Arrays .. DOUBLE PRECISION RESULT(101) C .. C .. External Subroutines .. EXTERNAL RJBESL C .. C .. Intrinsic Functions .. INTRINSIC INT integer i1mach, nout, nin nout = i1mach(2) nin = i1mach(1) C .. OPEN (90,FILE='TBESSE.out',STATUS='UNKNOWN',ACCESS='SEQUENTIAL', + ERR=30) WRITE(NOUT,FMT=9000) READ(NIN,FMT=*,ERR=20) X,NU,ALPHA C NB = INT(NU) + 1 WRITE (90,FMT=9010) WRITE(NOUT,FMT=9010) C CALL RJBESL(X,ALPHA,NB,RESULT,NCALC) C WRITE(NOUT,FMT=9020) NCALC WRITE(NOUT,FMT=9030) NB IF (NCALC.NE.NB) THEN WRITE(NOUT,FMT=9040) END IF DO 10 N = 1,NB WRITE (90,FMT=9050) N - 1 + ALPHA,X,RESULT(N) WRITE(NOUT,FMT=9050) N - 1 + ALPHA,X,RESULT(N) 10 CONTINUE C CLOSE (90) GO TO 40 20 WRITE(NOUT,FMT=9060) GO TO 40 30 WRITE(NOUT,FMT=9070) 40 STOP 9000 FORMAT (' MESSAGE TBESSE:',/,' GIVE X, NU, ALPHA',:,/) 9010 FORMAT ('# MESSAGE TBESSE: NU+ALPHA, X, BESSEL(NU+ALPHA,X)') 9020 FORMAT ('# MESSAGE TBESSE: NCALC =',I6) 9030 FORMAT ('# MESSAGE TBESSE: NB =',I6) 9040 FORMAT ('# MESSAGE TBESSE: NCALC NOT EQUAL NB!') 9050 FORMAT (D17.9,2X,D17.9,2X,D17.9) 9060 FORMAT (/,' MESSAGE TBESSE: ERROR DURING INPUT!') 9070 FORMAT (/,'MESSAGE TBESSE: ERROR ON OPENING FILE!') END C END OF SUBROUTINE *** TBESSE *** C C BEGIN OF SUBROUTINE *** TINTHP *** SUBROUTINE TINTHP C C ################################################################# C C FIRST VERSION: 05.02.1998 C PREVIOUS VERSION: 05.02.1998 C LATEST VERSION: 06.02.1998 C C PROGRAMMED BY: C THOMAS WIEDER C E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE C C PURPOSE: C TEST SUBROUTINE *** INTHP ***. C C ################################################################# C C IMPLICIT NONE C C .. Scalars in Common .. DOUBLE PRECISION EXACT,XSPLIT INTEGER IINT C .. C .. Local Scalars .. DOUBLE PRECISION A,B,DCLASS,DX,EPS,PCLASS,RESULT,X INTEGER I,INF,MCALL C .. C .. External Functions .. DOUBLE PRECISION TFUNC EXTERNAL TFUNC C .. C .. External Subroutines .. EXTERNAL INTHP C .. C .. Common blocks .. COMMON /LOSUNG/XSPLIT,EXACT,IINT integer i1mach, nout, nin nout = i1mach(2) nin = i1mach(1) C .. OPEN (81,FILE='test.plo',STATUS='unknown',ACCESS='sequential') C 10 WRITE(NOUT,FMT=9000) READ(NIN,FMT=*,ERR=10) IINT IF (IINT.EQ.0) GO TO 30 C X = 0.0D0 DX = 1.0D0 DO 20 I = 1,1000 WRITE (81,FMT=*) X,TFUNC(X) X = X + DX 20 CONTINUE C C DO AN INDEFINITE INTEGRAL FROM *** BOUND *** TO INFINITY: A = 0.0D0 B = 1.0D44 MCALL = 20000 C C *** 05.02.1998 *** C GOOD VALUE FOR *** EPS *** = 1.0D-3 EPS = 1.0D-03 C *** 05.02.1998 *** C C *** 05.02.1998 *** C GOOD VALUE FOR *** DCLASS *** = 10.0D0 C ALL OTHER VALUES RESULT IN FAILORE! DCLASS = 10.0D0 C *** 05.02.1998 *** C C *** 05.02.1998 *** C GOOD VALUE FOR *** PCLASS *** = 0.0D0 C ALL OTHER VALUES RESULT IN FAILORE! PCLASS = 0.0D0 C *** 05.02.1998 *** C C *** 05.02.1998 *** C TAKE *** INF = 2 *** FOR NON-OSCILLATING INTEGRALS! C TAKE *** INF = 3 *** FOR OSCILLATING INTEGRALS! INF = 3 C *** 05.02.1998 *** C CALL INTHP(A,B,DCLASS,TFUNC,MCALL,PCLASS,EPS,INF,RESULT) C WRITE(NOUT,FMT=*) 'RESULT =',RESULT WRITE(NOUT,FMT=*) 'NUMBER OF FUNCTION CALLS MACLL =',MCALL WRITE(NOUT,FMT=*) 'INFORMATION FROM INTHP =',INF WRITE(NOUT,FMT=*) 'EXACT SOLUTION:',EXACT C GO TO 10 C 30 RETURN 9000 FORMAT (1X,' GIVE IDENTIFICATION NUMBER OF INTEGRAND!',/, + ' (1,2,3 4,5,6,7,8,9, 0 = STOP):',/) END SHAR_EOF fi # end of overwriting check if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << SHAR_EOF > 'src.f' C BEGIN OF SUBROUTINE *** DHINIT *** SUBROUTINE DHINIT(LUININ,LUOUIN,LUERIN) C C ################################################################# C C FIRST VERSION: 10.10.1998 C PREVIOUS VERSION: 10.12.1998 C LATEST VERSION: 03.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 PURPOSE: C C STORE LOGICAL FILE NUMBERS IN COMMON BLOCK /INOUTE/ C TO TRANSFER THESE NUMBERS TO ALL OTHER SUBROUTINES C VIA COMMON BLOCK /INOUTE/. C C ################################################################# C C IMPLICIT NONE C C ----------------------------------------------------------------- C C INITIALIZE SOME VARAIBLES: C C .. Scalar Arguments .. INTEGER LUERIN,LUININ,LUOUIN C .. C .. Scalars in Common .. INTEGER IDAT,IREAD,LUERR,LUIN,LUOUT CHARACTER*80 FILEER,FILEIN,FILEOU C .. C .. Common blocks .. COMMON /INOUTE/LUIN,LUOUT,LUERR COMMON /LOGIK/IDAT,IREAD COMMON /NAMEN/FILEIN,FILEOU,FILEER C .. IREAD = 1 C C SAVE INPUT: C LUIN = LUININ LUERR = LUERIN LUOUT = LUOUIN C RETURN END C END OF SUBROUTINE *** DHINIT *** C C BEGIN OF SUBROUTINE *** HANKEL *** SUBROUTINE HANKEL(XI,NU,IVERBO,IERR1,HT) C C ################################################################## C C FIRST VERSION: 04.01.1998 C PREVIOUS VERSION: 14.10.1998 C LATEST VERSION: 12.12.1998 C C PROGRAMMED BY: C THOMAS WIEDER C TECHNICAL UNIVERSITY DARMSTADT C GERMANY C E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE C HTTP://WWW.UNI-KASSEL.DE/~WIEDER C C PURPOSE: C C SUBROUTINE *** HANKEL *** CALCULATES THE HANKEL TRANSFROM C OF A REAL-VALUED FUNCTION *** FUN ***. 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 FORM OF ** FUN *** IS: Y = FUN(X). C C ------------------------------------------------------------------ C C INPUT VARIABLES: C C ------------------------------------------------------------------ C C *** IDATIN ***: C C FOR *** IDAT=1 *** THE USER HAS TO PROVIDE A C SUBROUTINE *** DHFUNC *** C WHICH GIVES THE FUNCTION *** FUN ***: Y = FUN(X). C C THE FORM OF SUBROUTINE *** DHFUNC *** IS: C C SUBROUTINE DHFUNC(NU,X,Y,XLIMIT,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 C OUPUT VARIABLES: Y,XLIMIT,IERR 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 ------------------------------------------------------------------ C C FOR *** IDAT=2 *** THE USER HAS TO PROVIDE AN INPUT FILE C WHICH CONTAINS *** FUN *** AS A TABLE OF (X,Y)-PAIRS. C C THE FORM OF THE INPUT FILE IS: C C X_1 Y_1 C ... C X_N Y_N C ... C X_NEND Y_NEND C C THESE DATA WILL BE READ BY A CALL TO SUBROUTINE *** FUNK2 *** C PRIOR TO THE CALL TO SUBROUTINE *** HANKEL *** !!! !!! !!! C C ------------------------------------------------------------------ C C FOR *** IDAT=3 *** THE USER HAS TO PASS FUNCTION *** FUN *** C IN FORM OF TWO ARRAYS *** XDAT(N) ***, **** YDAT(N) *** C VIA A CALL TO SUBROUTINE *** FUNK3 *** C PRIOR TO THE CALL TO SUBROUTINE *** HANKEL *** !!! !!! !!! C C THE ARRAY *** XDAT(N) ***, *** YDAT(N) *** CONTAIN THE VALUES C DATA X_1,Y_1, ..., X_N,Y_N, ..., X_NEND,Y_NEND. C C *** XI ***: C C POSITION AT WHICH THE HANKEL TRANSFORM HAS TO BE EVALUATED. C *** XI *** IS A POSITIVE REAL VARIABLE (IN DOUBLE PRECISION). C C *** NU ***: C C *** NU *** IS THE ORDER OF THE HANKEL TRANSFORM. C *** NU *** MAY BE A POSITIVE OR NEGATIVE REAL VALUE C IN DOUBLE PRECISION. C HOWEVER, A MAXIMUM MAGNITUDE HAS BEEN SET FOR *** NU *** C WITHIN SUBROUTINE *** HANKEL ***. YOU MAY ALTER THIS SETTING. C C ------------------------------------------------------------------ C C *** IVERBO ***: C C A SWITCH BY WHICH THE USER DECIDES HOW MUCH INTERMEDIATE C OUTPUT WILL BE GENERATED BY SUBROUTINE *** HANKEL ***. C PLEASE SET *** IVERBO *** EQUAL TO 0, 1 OR 2. C C ------------------------------------------------------------------ C C *** LUOUT *** = LOGICAL NUMBER OF OUTPUT FILE C *** LUERR *** = LOGICAL NUMBER OF ERROR LOG FILE C C ------------------------------------------------------------------ C C OUTPUT VARIBLES: C C ------------------------------------------------------------------ C C *** IERR1 ***: C C *** IERR1 *** IS AN ERROR FLAG. C *** IRERR1 *** = 0 --> NO ERROR HAS BEEN DETECTED, C SUCCESSFULL CALCULATION. C *** IERR1 *** NOT EQUAL TO 0 --> AN ARROR OCCURED, C PLEASE CONSULT THE OUTPUT FILES! C C ------------------------------------------------------------------ C C *** HT ***: C C *** HT *** IS THE VALUE OF THE HANKEL TRANSFORM. C *** HT *** IS A POSITIVE OR NEGATIVE REAL VARIABLE C (IN DOUBLE PRECISION). C C ------------------------------------------------------------------ C C LITERATURE: C C HANDBUCH DER PHYSIK, S. FLUEGGE (HRSG.), C SPRINGER VERLAG, BERLIN, 1955, BAND II C I.N. SNEDDON: FUNCTIONAL ANALYSIS, PP. 198-348, C SECTION IV: THE HANKEL TRANSFORM, PP. 298-307. C C A. ERDELYI, W. MAGNUS, F. OBERHETTINGER, F.G. TRICOMI: C TABLES OF INTEGRAL TRANSFORMS, VOL. II, C McGRAW-HILL BOOK COMPANY, NEW YORK, 1954. C CHAPTER VIII: HANKEL TRANSFORMS. C C ------------------------------------------------------------------ C C SUCCESSFULLY COMPILED BY: C C 1.) AIX FORTRAN XLF (UNIX) C 2.) NAG FTN90 (DOS) C 3.) F90 NAG (UNIX) C 4.) DIGITAL VISUAL FORTRAN (WINDOWS NT) C C ------------------------------------------------------------------ C C CALLING SEQUENCE: C C CALLS SUBROUTINE *** FUNCHE ***. C CALLS SUBROUTINE *** HNKLNM ***. C CALLS SUBROUTINE *** DHINFO ***. C C ------------------------------------------------------------------ C C MOST IMPORTANT VARIABLES: C C XI = PARAMETER OF THE HANKEL TRANSFORM C = POSITION AT WHICH THE HANKEL TRANSFORM *** HT *** C OF FUNCTION *** FUN *** HAS TO BE EVALUATED. C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C XI IS AN INPUT VARIABLE TO SUBROUTINE *** HANKEL ***. C XI IS A POSITIVE, REAL-VALUED VARIABLE! C XI > 0 C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C NU = ORDER OF HANKEL TRANSFORM C = ORDER OF BESSEL FUNCTION *** DHJNUX ***. C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C NU IS AN INPUT VARIABLE TO SUBROUTINE *** HANKEL ***. C NU IS A REAL-VALUED VARIABLE! C NU MAY BE POSITIVE, ZERO, OR NEGATIVE. C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C HT IS AN OUTPUT VARIABLE FROM SUBROUTINE *** HANKEL ***. C HT IS A REAL-VALUED VARIBALE! C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C IDAT EQUAL TO 1 --> FUNCTION *** FUN *** IS GIVEN ANALYTICALLY. C IDAT EQUAL TO 2 --> FUNCTION *** FUN *** IS GIVEN AS TABULATED C DATA WITHIN A FILE. C IDAT EQUAL TO 3 --> FUNCTION *** FUN *** IS PASSED AS TABULATED C DATA VIA CALL. C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C IDAT IS AN INPUT VARIABLE TO SUBROUTINE *** HANKEL ***. C IDAT IS A INTEGER-VALUED VARIABLE! C IDAT MAY TAKE ON THE VALUES 1, 2 OR 3. C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C IERR1 = ERROR FLAG C IERR1 = 0 --> NO ERROR OCCURED DURING CALCULATION. C C NMAX = MAXIMUM NUMBER OF ALLOWED FUNCTION VALUES IN CASE OF C A TABULATED FUNCTION (*** IDAT *** = 2). C C ################################################################## C C IMPLICIT NONE C C .. Parameters .. DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) C .. C .. Scalar Arguments .. DOUBLE PRECISION HT,NU,XI INTEGER IERR1,IVERBO C .. C .. Scalars in Common .. DOUBLE PRECISION SLTN,XISLTN INTEGER IDAT,IFUNC,IREAD,LUERR,LUIN,LUOUT CHARACTER*80 FILEER,FILEIN,FILEOU C .. C .. Local Scalars .. DOUBLE PRECISION DPRCSN,NUMAX,NUMIN,XIMAX,XIMIN INTEGER IERR7,IERRH1,IHPKNS,IMOD CHARACTER*40 DATE C .. C .. External Functions .. DOUBLE PRECISION D1MACH EXTERNAL D1MACH C .. C .. External Subroutines .. EXTERNAL FUNCHE,HNKLNM C .. C .. Common blocks .. COMMON /INOUTE/LUIN,LUOUT,LUERR COMMON /LOGIK/IDAT,IREAD COMMON /NAMEN/FILEIN,FILEOU,FILEER COMMON /SOLUTI/XISLTN,SLTN,IFUNC integer i1mach, nout C .. C .. Save statement .. SAVE IHPKNS C .. nout = i1mach(2) DATE = ' FEBRUARY 1999' C C ------------------------------------------------------------------- C C SET SOME PARAMETERS: C IERR1 = 0 C CODE MODIFIED BY DR. T.R. HOPKINS, 28.04.1998: IHPKNS = 0 C C XISLTN IS NECCESSARY FOR DEBUGGING PURPOSES ONLY. XISLTN = XI C XIMIN = ZERO XIMAX = D1MACH(2) C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C THE CHOICE OF *** NUMIN *** AND *** NUMAX *** C IS SOMEWHAT ARBITRARY. YOU MAY MEET ANOTHER CHOICE. NUMIN = -100.0D0 NUMAX = 100.0D0 C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C IF (IVERBO.NE.0) WRITE(NOUT,FMT=9000) IF (IVERBO.GT.1) WRITE (LUERR,FMT=9000) IF (IVERBO.NE.0) WRITE(NOUT,FMT=9010) DATE IF (IVERBO.GT.1) WRITE (LUERR,FMT=9010) DATE C C ------------------------------------------------------------------- C C CHECK THE INPUT VALUES: C C IS *** XI *** POSITIVE? IF (XI.LE.ZERO) THEN WRITE(NOUT,FMT=9020) WRITE (LUERR,FMT=9020) WRITE(NOUT,FMT=9030) XI WRITE (LUERR,FMT=9030) XI IERR1 = 1 WRITE(NOUT,FMT=9130) C AN ERROR OCCURED, DO NOT CONTINUE, GO TO STOP! GO TO 20 END IF C C CHECK ABSOLUT VALUE OF *** XI *** IF (XI.LT.XIMIN .OR. XI.GT.XIMAX) THEN WRITE(NOUT,FMT=9020) WRITE (LUERR,FMT=9020) WRITE(NOUT,FMT=9040) XI WRITE (LUERR,FMT=9040) XI WRITE (LUERR,FMT=9050) XIMIN,XIMAX IERR1 = 1 WRITE(NOUT,FMT=9130) C AN ERROR OCCURED, DO NOT CONTINUE, GO TO STOP! GO TO 20 END IF C C CHECK ABSOLUT VALUE OF *** NU *** IF (NU.LT.NUMIN .OR. NU.GT.NUMAX) THEN WRITE(NOUT,FMT=9020) WRITE (LUERR,FMT=9020) WRITE(NOUT,FMT=9060) NU WRITE (LUERR,FMT=9060) NU WRITE (LUERR,FMT=9070) NUMIN,NUMAX IERR1 = 1 WRITE(NOUT,FMT=9130) C AN ERROR OCCURED, DO NOT CONTINUE, GO TO STOP! GO TO 20 END IF C IF (IDAT.EQ.1) THEN IF (IVERBO.NE.0) WRITE(NOUT,FMT=9080) IF (IVERBO.GT.1) WRITE (LUERR,FMT=9080) C ELSE IF (IDAT.EQ.2) THEN IF (IVERBO.NE.0) WRITE(NOUT,FMT=9090) IF (IVERBO.GT.1) WRITE (LUERR,FMT=9090) C C THIS IS FOR INTERNAL LOGIC, NEVER CHANGE THE FOLLOWING LINE: IFUNC = 0 C ELSE IF (IDAT.EQ.3) THEN IF (IVERBO.NE.0) WRITE(NOUT,FMT=9100) IF (IVERBO.GT.1) WRITE (LUERR,FMT=9100) C C THIS IS FOR INTERNAL LOGIC, NEVER CHANGE THE FOLLOWING LINE: IFUNC = 0 C ELSE WRITE(NOUT,FMT=9110) IDAT WRITE (LUERR,FMT=9110) IDAT WRITE(NOUT,FMT=9130) C AN ERROR OCCURED, DO NOT CONTINUE, GO TO STOP! GO TO 20 END IF C C ------------------------------------------------------------------- C C CHECK THE FUNCTION *** DHFUNC ***: C CALL FUNCHE(IERR7) C IF (IERR7.NE.0) THEN WRITE(NOUT,FMT=9120) FILEER WRITE (LUERR,FMT=9120) FILEER WRITE(NOUT,FMT=9130) WRITE(NOUT,FMT=9140) FILEER C C CALL DHINFO C STOP END IF C C ------------------------------------------------------------------- C C BEGIN OF CALCULATION! C IF (XI.EQ.ZERO) THEN HT = ZERO GO TO 10 END IF C IMOD = 1 C IF (IMOD.EQ.1) THEN C C !!! BEGIN OF MODE 1: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C MODE 1: HANKEL TRANSFORM BY EXPLICIT CALCULATION OF THE C TRANSFORMATION FORMULA. C I.E. HANKEL TRANSFORM BY INDFINITE INTEGRATION. C IF (IVERBO.NE.0) WRITE(NOUT,FMT=9150) IF (IVERBO.GT.1) WRITE (LUERR,FMT=9160) C CALL HNKLNM(XI,NU,IVERBO,IERRH1,HT) C IF (IVERBO.GT.1) WRITE (LUERR,FMT=9170) XI,NU,HT,IERRH1 C IF (IERRH1.NE.0) THEN WRITE(NOUT,FMT=9180) WRITE (LUERR,FMT=9180) IERR1 = 1 GO TO 20 END IF C C !!! END OF MODE 1: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C ELSE WRITE (LUERR,FMT=9190) IMOD IERR1 = -1 GO TO 20 C C END IF IF-CLAUSE IMOD END IF C C END OF CALCULATION! C C ------------------------------------------------------------------ C C WRITE OUT RESULTS: C 10 IF (IHPKNS.NE.1) THEN C C *** AMENDMENT SUGGESTED BY DR. T.R. HOPKINS, 15.04.1998 *** C WRITE OUT HEAD LINE ONLY ONCE: IHPKNS = 1 C IF (IFUNC.EQ.0 .OR. IDAT.EQ.2 .OR. IDAT.EQ.3) THEN IF (IVERBO.GT.0) WRITE (LUOUT,FMT=9200) ELSE IF (IFUNC.NE.0) THEN IF (IVERBO.GT.0) WRITE (LUOUT,FMT=9210) END IF C C END IF IF_CLAUSE *** IHPKNS *** END IF C IF (IFUNC.EQ.0 .OR. IDAT.EQ.2 .OR. IDAT.EQ.3) THEN C C IFUNC = 0 --> A USER SUPPLIED FUNCTION HAS BEEN TRANSFORMED: IF (IVERBO.GT.0) WRITE(NOUT,FMT=9220) XI,NU,HT IF (IVERBO.GT.0) WRITE (LUERR,FMT=9220) XI,NU,HT IF (IVERBO.GT.0) WRITE (LUOUT,FMT=9230) XI,HT C ELSE IF (IFUNC.NE.0) THEN C C IFUNC NOT 0 --> A TEST FUNCTION FROM WITHIN FUNCTION *** DHFUNC ** C HAS BEEN TRANSFORMED: WRITE(NOUT,FMT=9240) XISLTN,SLTN,HT WRITE (LUOUT,FMT=9250) XISLTN,SLTN,HT IF (IVERBO.GT.0) THEN C CALCULATE THE MACHINE PRECISION *** DPRSCN *** DPRCSN = D1mach(4) IF (SLTN.NE.0.0D0) THEN WRITE (LUERR,FMT=9260) XISLTN,SLTN,HT,HT - SLTN, + (HT-SLTN)/SLTN, (HT-SLTN)/DPRCSN, + (HT-SLTN)/SLTN/DPRCSN ELSE WRITE (LUERR,FMT=9270) XISLTN,SLTN,HT END IF END IF C C END IF-CLAUSE *** IFUNC *** END IF C 20 IF (IVERBO.NE.0) WRITE(NOUT,FMT=9280) FILEER IF (IVERBO.NE.0) WRITE(NOUT,FMT=9290) FILEOU C IF (IVERBO.NE.0) WRITE(NOUT,FMT=9300) IERR1 IF (IVERBO.GT.1) WRITE (LUERR,FMT=9300) IERR1 GO TO 30 C C ------------------------------------------------------------------ C C ERROR HANDLING: C 30 IF (IERR1.EQ.0) THEN IF (IVERBO.NE.0) WRITE(NOUT,FMT=9310) IF (IVERBO.GT.1) WRITE (LUERR,FMT=9310) C ELSE IF (IERR1.NE.0 .AND. IERR1.NE.-1) THEN WRITE(NOUT,FMT=9320) WRITE (LUERR,FMT=9320) C ELSE IF (IERR1.EQ.-1) THEN C C CALL DHINFO C WRITE(NOUT,FMT=9330) CLOSE (LUOUT) CLOSE (LUERR) STOP C C ENDE DER ABFRAGE AUF *** IERR1 ***. END IF C IF (IVERBO.GT.1) WRITE(NOUT,FMT=9340) IF (IVERBO.GT.1) WRITE (LUERR,FMT=9340) C C ------------------------------------------------------------------- C C *** SUBROUTINE DHINFO *** PROVIDES FURTHER INFORMATIONS: C C IF (IVERBO.GT.0) CALL DHINFO C C ------------------------------------------------------------------- C RETURN 9000 FORMAT (/,/,' ####################################',/, + ' MESSAGE HANKEL: BEGIN OF CALCULATION',/, + ' ####################################') 9010 FORMAT (/,' *** PROGRAM HANKEL *** ',/, + ' WRITTEN BY THOMAS WIEDER',/, + ' E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE',/, + ' DATE OF FIRST VERSION: A.D. 1998',/, + ' DATE OF LAST CHANGE:',A40) 9020 FORMAT (/,' MESSAGE HANKEL: WRONG INPUT VALUE!') 9030 FORMAT (/,' MESSAGE HANKEL:',/, + ' PARAMETER *** XI *** IS EQUAL TO ZERO OR NEGATIVE!',/, + ' *** XI *** =',D17.9) 9040 FORMAT (/,' MESSAGE HANKEL:',/, + ' PARAMETER *** XI *** OUTSIDE ALLOWED RANGE!',/, + ' *** XI *** =',D17.9) 9050 FORMAT (/,' ALLOWED INTERVAL [XI_MIN,XI_MAX]:',/,'[',D17.9,',', + D17.9,']') 9060 FORMAT (/,' MESSAGE HANKEL:',/, + ' PARAMETER *** NU *** OUTSIDE ALLOWED RANGE!',/, + ' *** NU *** =',D17.9) 9070 FORMAT (/,' ALLOWED INTERVAL [NU_MIN,NU_MAX]:',/,' [',D17.9,',', + D17.9,']') 9080 FORMAT (/,' MESSAGE HANKEL:',/, + ' WILL TRY TO TRANSFORM *** FUNCTION FUNC ***!') 9090 FORMAT (/,' MESSAGE HANKEL:',/, + ' WILL TRY TO TRANSFORM TABULATED DATA FROM INPUT FILE!') 9100 FORMAT (/,' MESSAGE HANKEL:',/, +' WILL TRY TO TRANSFORM TABULATED DATA TRANSMITTED BY CALL TO HANK +EL!') 9110 FORMAT (/,' MESSAGE HANKEL',/, + ' WRONG VALUE FOR CALL PARAMETER *** IDAT *** !',/, + ' IDAT = ',I6) 9120 FORMAT (/,' MESSAGE HANKEL:',/, +' A TEST OF YOUR FUNCTION *** DHFUNC *** GAVE AN ERRONEOUS R +ESULT!',/,' PLEASE HAVE A LOOK INTO FILE',/,A80) 9130 FORMAT (/,' MESSAGE HANKEL: NO TRANSFORMATION CALCULATED! STOP!') 9140 FORMAT (/,' MESSAGE HANKEL: HAVE A LOOK INTO FILE:',/,A80) 9150 FORMAT (/, +' MESSAGE HANKEL: WILL TRY NUMERICAL HANKEL TRANSFORM...PLEASE WAI +T!') 9160 FORMAT (/, + ' MESSAGE HANKEL: WILL CALL SUBROUTINE *** HNKLNM ** * !' + ) 9170 FORMAT (/, + ' MESSAGE HANKEL: XI NU HT ERROR FLAG:' + ,/,D17.9,3X,D17.9,3X,D17.9,3X,I6) 9180 FORMAT (/,' MESSAGE HANKEL: ERROR DURING TRANSFORMATION!',/, + ' ERROR IN *** SUBROUTINE HNKLNM *** !') 9190 FORMAT (/,' MESSAGE HANKEL:',/,' WRONG CALCULATION MODUS CHOOSEN!' + ,/,' IMOD = ',I6) 9200 FORMAT ('# XI HANKEL_TRANSFORM(XI)') 9210 FORMAT ( +'# XI EXACT SOLUTION NUMERICAL SOL +UTION') 9220 FORMAT (/,' MESSAGE HANKEL:',/, +' ARGUMENT XI ORDER NU NUMERICAL HANKE +L TRANSFORM HT:',/,1X,D17.9,8X,D17.9,8X,D17.9) 9230 FORMAT (E17.9,5X,E17.9) 9240 FORMAT (/, + ' ARGUMENT XI, ANALYTICAL TRANSFORM, NUMERICAL SOLUTION:' + ,/,D17.9,3X,D17.9,3X,D17.9) 9250 FORMAT (E17.9,5X,E17.9,5X,E19.7) 9260 FORMAT (1X,D17.9,1X,D17.9,1X,D17.9,1X,D17.9,1X,D17.9,1X,D17.9,1X, + D17.9) 9270 FORMAT (3X,D17.9,3X,D17.9,3X,D17.9) 9280 FORMAT (/,' MESSAGE HANKEL: FOR ERROR MESSAGES LOOK INTO FILE:',/, + A80) 9290 FORMAT (' MESSAGE HANKEL: FOR RESULTS LOOK INTO FILE:',/,A80) 9300 FORMAT (' IERR1 =',I6) 9310 FORMAT (/,' MESSAGE HANKEL: SUCCESSFUL END OF CALCULATION!',/, + ' NO ERRORS HAVE BEEN DETECTED DURING CALCULATION! GOOD BYE!' + ) 9320 FORMAT (/,' MESSAGE HANKEL:',/,' UNSUCCESSFUL END OF CALCULATION!' + ,/,' ERRORS HAVE BEEN DETECTED DURING CALCULATION!',/, + ' GOOD BYE!') 9330 FORMAT (/,/,' MESSAGE HANKEL:',/, + ' END OF CALCULATION BECAUSE OF UNRCOVERABLE ERROR!',/, + ' GOOD BYE!') 9340 FORMAT (/,' ####################################',/, + ' MESSAGE HANKEL: END OF CALCULATION',/, + ' ####################################',/,/) END C END OF SUBROUTINE *** HANKEL *** C C BEGIN OF SUBROUTINE *** HNKLNM *** SUBROUTINE HNKLNM(XIIN,NUIN,IVERBO,IERR2,HT) C C ################################################################## C C FIRST VERSION: 21.01.1998 C PREVIOUS VERSION: 10.12.1998 C LATEST VERSION: 05.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 PURPOSE: C C SUBROUTINE *** HNKLNM *** CALCULATES THE HANKEL TRANSFROM C OF A REAL-VALUED FUNCTION *** DHFUNC *** C BY EXPLICIT NUMERICAL EVALUATION OF THE TRANSFORMATION FORMULA. C C LITERATURE: C C A. ERDELYI, W. MAGNUS, F. OBERHETTINGER, F.G. TRICOMI: C TABLES OF INTEGRAL TRANSFORMS, VOL. II, C McGRAW-HILL BOOK COMPANY, NEW YORK, 1954. C CHAPTER VIII: HANKEL TRANSFORMS. C C MOST IMPORTANT VARIABLES: C C XI = PARAMETER OF THE HANKEL TRANSFORM. C THE HANKEL TRANSFORM OF *** FUN *** WILL BE CALCULATED C FOR ARGUMENT *** XI ***. C *** FUN *** DEPENDS ON VARIABLE *** X ***. C HT = HANKEL TRANSFORM OF FUNCTION *** FUN ***. C C CALLING SEQUENCE: C C CALLED BY SUBROUTINE *** HANKEL ***. C CALLS FUNCTION *** FUNINT ***. C CALLS SUBROUTINE *** INTHP ***. C CALLS SUBROUTINE *** OSCINT ***. C C ################################################################## C C IMPLICIT NONE C C INTEGER VARIABLES FOR SUBROUTINE *** INTHP *** C C INTEGER VARIABLES FOR SUBROUTINE *** OSCINT *** C PARAMETERS FOR SUBROUTINE *** OSCINT *** C *** 19.02.1998 *** SUCCESSFUL VALUE FOR *** NDIM1 *** = 1000 C *** 19.02.1998 *** SUCCESSFUL VALUE FOR *** NDIM2 *** = 28 C C C C C REAL VARIABLES FOR SUBROUTINE *** INTHP *** C C REAL VARAIBLES FOR SUBROUTINE *** OSCINT *** C C VARIABLES FOR TESTING SUBROUTINE *** INTHP *** C 1 ,TFUNC C C C COMMON BLOCK FOR SUBROUTINE *** INTHP *** C C C EXTERNAL FUNCTION FOR SUBROUTINE *** INTHP *** C EXTERNAL TFUNC C C EXTERNAL FUNCTION FOR SUBROUTINE *** OSCINT *** C C ------------------------------------------------------------------- C C SAVE INPUT C C .. Parameters .. INTEGER NDIM1,NDIM2 PARAMETER (NDIM1=1000,NDIM2=28) INTEGER NMAX PARAMETER (NMAX=1000) C .. C .. Scalar Arguments .. DOUBLE PRECISION HT,NUIN,XIIN INTEGER IERR2,IVERBO C .. C .. Scalars in Common .. DOUBLE PRECISION EXACT,NU,XDATMA,XDATMI,XI,XLIMIT,XSPLIT,YDATMA, + YDATMI INTEGER IDAT,IERR3,IINT,IREAD,LUERR,LUIN,LUOUT,NEND C .. C .. Arrays in Common .. DOUBLE PRECISION XDAT(NMAX),YDAT(NMAX) C .. C .. Local Scalars .. DOUBLE PRECISION ALIMIT,BLIMIT,DCLASS,DUMMY,EPS,PCLASS,PERIOD,PI, + RESULT,RFIRST,XASYMP INTEGER IERR4,IERRME,INF,INTPAR,INTSUB,MCALL,N1,N2,NQUAD C .. C .. Local Arrays .. DOUBLE PRECISION ABSCIS(36),QLIST(0:NDIM1),SAVPER(0:256), + WEIGHT(36),WORK(0:NDIM1,0:NDIM2) INTEGER ISTATE(6) C .. C .. External Functions .. DOUBLE PRECISION FUNINT,GPER,HFUN EXTERNAL FUNINT,GPER,HFUN C .. C .. External Subroutines .. EXTERNAL ERRMSS,FUNK1,FUNK2,FUNK3,G5AND9,INTHP,OSCINT C .. C .. Intrinsic Functions .. INTRINSIC DACOS C .. C .. Common blocks .. COMMON /ARGUME/XI,NU,XLIMIT COMMON /DATEN/XDAT,YDAT,XDATMI,XDATMA,YDATMI,YDATMA,NEND COMMON /FHNKL1/IERR3 COMMON /INOUTE/LUIN,LUOUT,LUERR COMMON /LOGIK/IDAT,IREAD COMMON /LOSUNG/XSPLIT,EXACT,IINT integer i1mach, nout nout = i1mach(2) C .. XI = XIIN NU = NUIN C C ------------------------------------------------------------------- C C SET SOME PARAMETERS: C IERR2 = 0 IERR3 = 0 IERR4 = 0 HT = 0.0D0 RESULT = 0.0D0 PI = 2.0D0*DACOS(0.0D0) DUMMY = 0.0D0 C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C *** 24.02.1998 *** A GOOD VALUE FOR *** XASYMP *** = 100.0 XASYMP = 100.0D0 C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C CHOOSE A TEST FUNCTION WITHIN SUBROUTINE *** TINTHP *** C BY SETTING *** IINT ***: IINT = 9 C C ------------------------------------------------------------------- C C BEGIN THE CALCULATION: C C WE SPLIT THE INFINITE INTEGRAL INTO TWO PARTS. C THE FIRST PART IS A FINITE INTEGRAL FROM 0 TO *** XSPLIT ***. C THE SECOND PART IS AN INFINITE INTEGRAL FROM *** XSPLIT *** C TO INFINITY. C C WE CHOOSE *** XSPLIT *** IN A WAY THAT THE OSCILLATING PART C IN THE INFINITE INTEGRAL HAS A CONSTANT PERIOD OF OSCILLATION!!! C XSPLIT = XASYMP C C IF THE FUNCTION *** FUN *** IS EQUAL TO ZERO FOR *** X *** C GREATER THAN *** XLIMIT *** C OR THE FUNCTION *** FUN *** MAY BETTER BE INTEGRATED FROM C *** XLIMIT *** TO INFINITY THAN FROM *** XASYMP *** C TO INFINITY, C THEN WE SET *** XSPLIT = XLIMIT ***. C C FIRST, CALL SUBROUTINE *** FUNK ** TO GET THE VALUE OF C *** XLIMIT ***. C *** XLIMIT *** APPEARS IN COMMON BLOCK *** /ARGUMENTE/ ***. C XLIMIT = -1.0D0 C IF (IDAT.EQ.1) THEN CALL FUNK1(DUMMY,DUMMY,IERR4) ELSE IF (IDAT.EQ.2) THEN CALL FUNK2(DUMMY,DUMMY,IERR4) ELSE IF (IDAT.EQ.3) THEN CALL FUNK3(DUMMY,DUMMY,XDAT,YDAT,NEND,IERR4) END IF C IF (IERR4.NE.0) THEN WRITE (LUERR,FMT=9000) END IF C IF (XLIMIT.LT.0.0D0) THEN IF (IVERBO.GT.1) THEN WRITE(NOUT,FMT=9010) XSPLIT WRITE (LUERR,FMT=9010) XSPLIT END IF END IF C C NOW WE KNOW *** XLIMIT ***. C IF (XLIMIT.GT.0.0D0) THEN IF (XLIMIT.GT.XASYMP) THEN IF (IVERBO.GT.1) WRITE (LUERR,FMT=9030) XLIMIT IF (IVERBO.GT.1) WRITE (LUERR,FMT=9020) XLIMIT,XASYMP, + XASYMP ELSE XSPLIT = XLIMIT IF (IVERBO.GT.1) WRITE (LUERR,FMT=9030) XLIMIT C END OF IF-CLAUSE *** XASYMP *** END IF C END OF IF-CLAUSE *** XLIMIT *** END IF C IF (IVERBO.GT.1) WRITE (LUERR,FMT=9040) XSPLIT C C *** INPART *** = 1 --> DO THE FINITE INTEGRAL! INTPAR = 1 C 10 IF (INTPAR.LE.2) THEN IF (INTPAR.EQ.1) THEN C DO THE FINITE INTEGRAL: C CHOOSE INTEGRATION SUBROUTINE BY SETTING *** INTSUB ***: INTSUB = 1 ELSE IF (INTPAR.EQ.2) THEN C DO THE INFINITE INTEGRAL: INTSUB = 2 ELSE WRITE (LUERR,FMT=9050) INTPAR STOP C END OF IF-CLAUSE *** INTPAR ***. END IF ELSE WRITE (LUERR,FMT=9050) INTPAR C END OF IF-CLAUSE *** INTPAR *** .LT. 2 END IF C IF (INTSUB.EQ.1) THEN C *** SBROUTINE INTHP *** CALCULATES A FINITE OR C AN INFINITE INTEGRAL. ALIMIT = 0.0D0 BLIMIT = XSPLIT DCLASS = 10.0D0 MCALL = 200000 PCLASS = 0.0D0 C 01.10.1998 EPS=1.0D-06 --> 1.0D-09 EPS = 1.0D-09 C *** INF *** = 3 --> INFINITE OSZILLATING TAIL C *** INF *** = 4 --> FINITE INTEGRAL INF = 4 C CALL INTHP(ALIMIT,BLIMIT,DCLASS,FUNINT,MCALL,PCLASS,EPS,INF, + RESULT) C IF (IVERBO.GT.1) WRITE (LUERR,FMT=9060) RESULT C IF (INF.GT.3) THEN WRITE (LUERR,FMT=9070) INF IERR2 = 1 END IF C ELSE IF (INTSUB.EQ.2) THEN C SUBROUTINE *** OSCINT *** CALCULATES AN INFINITE, OSCILLATING C INTEGRAL. C C THE FOLLOWING FEW LINES ARE TAKEN FROM SOURCE CODE OF ACM 639: C C SET THE INPUT PARAMETERS FOR OSCINT. C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C THE PERIOD OF THE ASYMPTOTIC EXPANSION OF THE BESSEL FUNCTION: C THIS EXPANSION IS THE COSINUS: PERIOD = 2.0D0*PI/XI C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C ALIMIT (AZERO) IS THE LOWER LIMIT OF INTEGRATION ALIMIT = XSPLIT C C RFIRST IS THE RIGHT-HAND ENDPOINT OF THE FIRST INTERVAL RFIRST = 0.0D0 C C NQUAD REPRESENTS THE TYPE OF QUADRATURE RULE USED C NQUAD = -5 5 PT GAUSS RULE C NQUAD = -9 9 PT GAUSS RULE C NQUAD = N,N>0 (N ODD) N POINT TRAPAZOIDAL RULE(N-1 PANELS) NQUAD = -9 C C SET EPSILON, THE REQUESTED TOLERANCE C *** T. WIEDER, 06.03.98 *** A GOOD VALUE FOR *** EPS *** = 1.0D-03 EPS = 1.0D-03 C C NDIM1 AND NDIM2 ARE THE DIMENSIONS OF THE WORK TABLE. THEY C MUST BE LESS THAN OR EQUAL TO THE DIMENSIONS GIVEN TO WORK C IN THE DIMENSION STATEMENT C DO 30 N1 = 0,NDIM1 QLIST(N1) = 0.0D0 DO 20 N2 = 0,NDIM2 WORK(N1,N2) = 0.0D0 20 CONTINUE 30 CONTINUE DO 40 N1 = 1,36 SAVPER(N1) = 0.0D0 WEIGHT(N1) = 0.0D0 ABSCIS(N1) = 0.0D0 40 CONTINUE C CALL OSCINT(ALIMIT,PERIOD,RFIRST,EPS,NQUAD,NDIM1,NDIM2,G5AND9, + HFUN,GPER,WORK,SAVPER,WEIGHT,ABSCIS,QLIST,RESULT, + ISTATE) C IF (IVERBO.GT.1) WRITE (LUERR,FMT=9080) RESULT C IF (ISTATE(1).LT.0) THEN C AN ERROR HAS OCCURRED IN SUBROUTINE *** OSCINT ***. WRITE (LUERR,FMT=9090) WRITE (LUERR,FMT=9100) ISTATE(1),ISTATE(2),ISTATE(3), + ISTATE(4),ISTATE(5),ISTATE(6) IERR2 = 1 C RESET THE UNRELIABLE RESULT: RESULT = 0.0D0 END IF C ELSE IF (INTSUB.EQ.3) THEN C SUBROUTINE *** DQAINF *** CALCULATES AN INFINITE INTEGRAL. WRITE(NOUT,FMT=9110) WRITE (LUERR,FMT=*) 55 C C CALL DTESTI (IERR3,RESULTS) C IERR3 = 1 IF (IERR3.GT.0) THEN WRITE (LUERR,FMT=9120) IERR3 IERR2 = 1 END IF C C END OF IF-CLAUSE *** INTSUB *** END IF C HT = HT + RESULT INTPAR = INTPAR + 1 C CHECK WHETHER THE INFINITE INTEGRAL STILL HAS TO BE EVALUATED: IF (INTPAR.LE.2) GO TO 10 C IF (IERR3.NE.0) THEN C AN ERROR IN FUNCTION *** FUNINT *** HAS OCCURED. WRITE(NOUT,FMT=9130) WRITE (LUERR,FMT=9130) WRITE (LUERR,FMT=9140) XI WRITE (LUERR,FMT=9150) NU IERRME = 10 C CALL ERRMSS(IERRME) C END IF C IF (IERR2.NE.0) THEN C AN ERROR DURING INTEGRATION HAS OCCURED! WRITE(NOUT,FMT=9160) WRITE (LUERR,FMT=9160) WRITE (LUERR,FMT=9140) XI WRITE (LUERR,FMT=9150) NU IERRME = 11 C CALL ERRMSS(IERRME) C END IF C 9000 FORMAT (/,' MESSAGE HNKLNM:',/, + ' ERROR MESSAGE FROM SUBROUTINE * ** FUNK *** ',/, + ' WHILE TRYING TO GET THE VALUE FOR *** XLIMIT * ** !',/, + ' HOWEVER, THE CALCULATION MAY STILL TURN OUT WELL...') 9010 FORMAT (/,' MESSAGE HNKLNM:',/, + ' *** XLIMIT *** LESS OR EQUAL TO ZERO!',/, +' *** XLIMIT HAS BEEN REJECTED AS UPPER LIMIT *** XSPLIT *** FOR +THE FINITE INTEGRAL',/,' *** XSPLIT *** =',D17.9) 9020 FORMAT (/,' MESSAGE HNKLNM: *** XLIMIT.GT.XASYMP *** !',/, + ' *** XLIMIT *** =',D17.9,3X,' XASYMP =',D17.9,/, + ' WILL NOT ACCEPT REQUESTED VALUE FOR *** XLIMIT ***!',/, + ' WILL USE AS UPPER LIMIT *** XSPLIT *** FOR THE FINITE INTEGRAL' + ,/,' THE VALUE OF *** XASYMP *** =',D17.9) 9030 FORMAT (/,' MESSAGE HNKLNM:',/, +' SUBROUTINE *** FUNK (OR DHFUNC) *** REQUESTS TO TAKE THE FINIT +E',/,' INTEGRAL FROM 0 TO *** XLI MIT *** =',D17.9) 9040 FORMAT (/,' MESSAGE HNKLNM: ',/, +' WILL CALCULTATE THE FINITE INTEGRAL FROM 0 TO *** XSPLIT * +** =',D17.9) 9050 FORMAT (/,' MESSAGE HNKLNM: ERROR IN INTERNAL LOGIC!',/, + ' *** INTPAR *** =',I6) 9060 FORMAT (/, + ' MESSAGE HNKLNM: *** INTHP *** RETURNS *** RESULT * ** =' + ,D17.9) 9070 FORMAT (/,' MESSAGE HNKLNM:',/, + ' ERROR ON RETURN FROM INTEGRAT ION SUBROUTINE *** INTHP ***!' + ,/,' *** INF *** =',I6) 9080 FORMAT (/, + ' MESSAGE HNKLNM: *** OSCINT *** RETURNS *** RESULT *** =' + ,D17.9) 9090 FORMAT (/,' MESSAGE HNKLNM:',/, +' ERROR ON RETURN FROM INTEGRAT ION SUBROUTINE *** OSCINT ***!' + ,/,' *** ISTATE *** =',I6) 9100 FORMAT (/,' ISTATE(1)=',I6,2X,' ISTATE(2)=',I6,2X,' ISTATE(3)=', + I6,2X,' ISTATE(4)=',I6,2X,' ISTATE(5)=',I6,2X, + ' ISTATE(6)=',I6) 9110 FORMAT (/,' MESSAGE HNKLNM:',/, + ' SUBROUTINE *** DQAINF *** IS NOT IMPLEMENTED!') 9120 FORMAT (/,' MESSAGE HNKLNM:',/, +' ERROR ON RETURN FROM INTEGRAT ION SUBROUTINE *** DQAINF ***!' + ,/,' *** IFAIL *** =',I6) 9130 FORMAT (/,' MESSAGE HNKLNM:',/, + ' AN ERROR OCCURED DURING EVALU ATION OF THE INTEGRAND!' + ) 9140 FORMAT (/,' MESSAGE HNKLNM: PARAMETER XI = ',D19.7) 9150 FORMAT (/,' MESSAGE HNKLNM: PARAMETER NU = ',D19.7) 9160 FORMAT (/, + ' MESSAGE HNKLNM: AN ERROR OCCURED DURING INTEGRATIO N!' + ) END C END OF SUBROUTINE *** HNKLNM *** C C BEGIN OF FUNCTION *** FUNINT *** DOUBLE PRECISION FUNCTION FUNINT(XIN) C C FIRST VERSION: 21.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 FUNCTION *** FUNINT *** CALCULATES THE INTEGRAND C (XI * X)**(1/2) JNUX(XI * X) FUN(X) C OF THE HANKEL TRANSFORMATION INTEGRAL. C C CALLING SEQUENCE: C CALLED BY SUBROUTINE *** HNKLNM ***. C CALLS SUBROUTINE *** FUNK *** C C MOST IMPORTANT VARIABLES: C C X = INDEPENDENT REAL VARIABLE OF *** FUN ***. C Y = DEPENDENT REAL VARIABLE OF *** FUN ***. C Y = FUN(X) C C XI = PARAMETER OF THE HANKEL TRANSFORM. C THE HANKEL TRANSFORM OF *** FUN *** WILL BE CALCULATED C FOR PARAMETER *** XI ***. C C NU = ORDER OF BESSEL FUNCTION *** DHJNUX ***. C C FUNINT = OUTPUT VALUE = INTEGRAND. C C IERR3 = ERROR FLAG C IERR3 = 0 --> NO ERROR OCCURED DURING EVALUATION OF C FUNCTION *** FUNINT ***. C C ################################################################## C C IMPLICIT NONE C C SAVE INPUT C C .. Parameters .. INTEGER NMAX PARAMETER (NMAX=1000) C .. C .. Scalar Arguments .. DOUBLE PRECISION XIN C .. C .. Scalars in Common .. DOUBLE PRECISION NU,XDATMA,XDATMI,XI,XLIMIT,YDATMA,YDATMI INTEGER IDAT,IERR3,IREAD,LUERR,LUIN,LUOUT,NEND C .. C .. Arrays in Common .. DOUBLE PRECISION XDAT(NMAX),YDAT(NMAX) C .. C .. Local Scalars .. DOUBLE PRECISION JNUXV,NUOUT,X,XIX,Y INTEGER IERR4,IERRJ C .. C .. External Subroutines .. EXTERNAL DHJNUX,FUNK1,FUNK2,FUNK3 C .. C .. Intrinsic Functions .. INTRINSIC DSQRT C .. C .. Common blocks .. COMMON /ARGUME/XI,NU,XLIMIT COMMON /DATEN/XDAT,YDAT,XDATMI,XDATMA,YDATMI,YDATMA,NEND COMMON /FHNKL1/IERR3 COMMON /INOUTE/LUIN,LUOUT,LUERR COMMON /LOGIK/IDAT,IREAD integer i1mach, nout nout = i1mach(2) C .. X = XIN NUOUT = NU C XIX = XI*X C C CALCULATE THE BESSEL FUNCTION *** JNUX *** C OF ORDER *** NU *** FOR ARGUMENT *** XIX ***. IERRJ = 0 C CALL DHJNUX(XIX,NUOUT,IERRJ,JNUXV) C IF (IERRJ.NE.0) THEN WRITE (LUERR,FMT=9000) WRITE (LUERR,FMT=9010) XIX,XI,X,NU IERR3 = 1 WRITE(NOUT,FMT=9020) WRITE (LUERR,FMT=9020) STOP END IF C C CALCULATE *** FUN *** FOR ARGUMENT *** X ***. C IF (IDAT.EQ.1) THEN CALL FUNK1(X,Y,IERR4) ELSE IF (IDAT.EQ.2) THEN CALL FUNK2(X,Y,IERR4) ELSE IF (IDAT.EQ.3) THEN CALL FUNK3(X,Y,XDAT,YDAT,NEND,IERR4) END IF C IF (IERR4.NE.0) THEN WRITE (LUERR,FMT=9030) WRITE (LUERR,FMT=9040) X WRITE (LUERR,FMT=9050) IERR4 IERR3 = 1 FUNINT = 0.0D0 WRITE(NOUT,FMT=9020) WRITE (LUERR,FMT=9020) STOP END IF C C FINALLY CALCULATE THE INTEGRAND: FUNINT = DSQRT(XIX)*JNUXV*Y C RETURN 9000 FORMAT (/,' MESSAGE FUNINT:',/, + ' ERROR DURING BESSEL FUNCTION CALC ULATION!') 9010 FORMAT (' MESSAGE FUNINT:',/,' ARGUMENTS TO BESSEL FUNCTION:',/, + ' XI * X =',D17.9,/,' XI = ',D17.9,' X =',D17.9,/, + ' NU =',D17.9) 9020 FORMAT (/,' MESSAGE FUNINT:',/, + ' AN ERROR DURING FUNCTION EVALUATION HAS OCCURED!',/, + ' THE HANKEL TRANSFORM HAS FAILED! SORRY! STOP!') 9030 FORMAT (/,' MESSAGE FUNINT:',/, + ' ERROR DURING CALCULATION OF FUNCTION *** FUN *** !') 9040 FORMAT (/,' MESSAGE FUNINT:',/, + ' ARGUMENT TO FUNCTION *** FUN ** * :',/,' X =',D17.9) 9050 FORMAT (1X,' *** IERR4 *** =',I6) END C END OF FUNCTION *** FUNINT *** C C BEGIN OF SUBROUTINE *** DHRW *** SUBROUTINE DHRW(IERR) C C ################################################################## C C BEGONNEN AM: 09.12.96 C VORLETZTE AENDERUNG: 01.02.99 C LETZE AENDERUNG AM: 08.02.99 C C PROGRAMMED BY: C THOMAS WIEDER C E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE C C PURPOSE: C LESEN VON X-Y-DATEN AUS EINER EINGABEDATEI. C SORTIEREN VON X-Y-DATEN. C SCHREIBEN VON X-Y-DATEN IN EINE AUSGABEDATEI. C C EINGABEDATEN IM FORMAT: C X_1 Y_1 C X_2 Y_2 C ... C X_N Y_N C C LUIN = LOGICAL NUMBER OF INPUT FILE C C IERR = 0 --> NO ERROR OCCURED. C IERR = 1 --> ERROR WHILE TRYING TO OPEN THE INPUT FILE. C IERR = 2 --> ERROR WHILE READING FROM INPUT FILE. C C ################################################################## C C IMPLICIT NONE C C .. Parameters .. INTEGER NMAX PARAMETER (NMAX=1000) C .. C .. Scalar Arguments .. INTEGER IERR C .. C .. Scalars in Common .. DOUBLE PRECISION XDATMA,XDATMI,YDATMA,YDATMI INTEGER LUERR,LUIN,LUOUT,NEND CHARACTER*80 FILEER,FILEIN,FILEOU C .. C .. Arrays in Common .. DOUBLE PRECISION XDAT(NMAX),YDAT(NMAX) C .. C .. Local Scalars .. INTEGER N C .. C .. External Functions .. DOUBLE PRECISION D1MACH EXTERNAL D1MACH C .. C .. External Subroutines .. EXTERNAL DHSORT C .. C .. Intrinsic Functions .. INTRINSIC MAX,MIN C .. C .. Common blocks .. COMMON /DATEN/XDAT,YDAT,XDATMI,XDATMA,YDATMI,YDATMA,NEND COMMON /INOUTE/LUIN,LUOUT,LUERR COMMON /NAMEN/FILEIN,FILEOU,FILEER integer i1mach, nout nout = i1mach(2) C .. IERR = 0 XDATMI = D1MACH(2) XDATMA = -D1MACH(2) YDATMI = D1MACH(2) YDATMA = -D1MACH(2) C WRITE(NOUT,FMT=9000) C C EINLESEN: C C20 WRITE (*,30) C30 FORMAT (/,' MESSAGE RW: GIVE NAME OF INPUT FILE:',/) C READ (*,40,ERR=20) NAME C40 FORMAT (A40) C FILEIN=NAME C DATEI OEFFNEN: OPEN (LUIN,FILE=FILEIN,STATUS='OLD',ACCESS='SEQUENTIAL',ERR=40) C NEND = 0 DO 10 N = 1,NMAX READ (LUIN,FMT=*,ERR=50,END=20) XDAT(N),YDAT(N) NEND = NEND + 1 IF (NEND.EQ.NMAX) THEN WRITE(NOUT,FMT=9010) NEND GO TO 20 END IF XDATMI = MIN(XDATMI,XDAT(N)) XDATMA = MAX(XDATMA,XDAT(N)) YDATMI = MIN(YDATMI,YDAT(N)) YDATMA = MAX(YDATMA,YDAT(N)) 10 CONTINUE C 20 WRITE(NOUT,FMT=9020) NEND C IF (XDATMI.LT.-D1MACH(2)) THEN WRITE(NOUT,FMT=9030) WRITE (LUERR,FMT=9030) IERR = 1 END IF IF (XDATMI.GT.D1MACH(2)) THEN WRITE(NOUT,FMT=9040) WRITE (LUERR,FMT=9040) IERR = 1 END IF IF (YDATMI.LT.-D1MACH(2)) THEN WRITE(NOUT,FMT=9050) WRITE (LUERR,FMT=9050) IERR = 1 END IF IF (YDATMI.GT.D1MACH(2)) THEN WRITE(NOUT,FMT=9060) WRITE (LUERR,FMT=9060) IERR = 1 END IF C C DATEI SCHLIESSEN: CLOSE (LUIN) C C NACH AUFSTEIGENDEN X-WERTEN SORTIEREN: C CALL DHSORT(XDAT,YDAT,NEND) C DO 30 N = 1,NEND WRITE (LUERR,FMT=9070) N,XDAT(N),YDAT(N) 30 CONTINUE C WRITE(NOUT,FMT=9080) C GO TO 60 40 WRITE(NOUT,FMT=*) 'ERROR ON OPENING FILE!' WRITE (LUERR,FMT=*) 'ERROR ON OPENING FILE! FILE DOES NOT EXIST?' IERR = 1 RETURN 50 WRITE(NOUT,FMT=*) 'ERROR ON READING DATA!' IERR = 2 RETURN C 60 RETURN 9000 FORMAT (/,' MESSAGE RW: BEGIN OF DATA INPUT FROM FILE!') 9010 FORMAT (/,' MELDUNG RW:',/,' WARNING! MAXIMAL NUMBER =',I6, + ' OF DATA PAIRS HAVE BEEN READ!') 9020 FORMAT (/,' MESSAGE RW:',I6,' DATA PAIRS HAVE BEEN READ!') 9030 FORMAT (/,' MESSAGE RW: TOO SMALL VALUE FOR ABCISSA *** X ***!',/, + ' PLEASE CHECK YOUR DATA!') 9040 FORMAT (/,' MESSAGE RW: TOO BIG VALUE FOR ABCISSA *** X ***!',/, + ' PLEASE CHECK YOUR DATA!') 9050 FORMAT (/,' MESSAGE RW: TOO SMALL VALUE FOR ORDINATE *** Y ***!', + /,' PLEASE CHECK YOUR DATA!') 9060 FORMAT (/,' MESSAGE RW: TOO BIG VALUE FOR ORDINATE *** X ***!',/, + ' PLEASE CHECK YOUR DATA!') 9070 FORMAT (' MESSAGE RW: N, X(N), Y(N):',I6,3X,D17.9,3X,D17.9) 9080 FORMAT (/,' MESSAGE RW: END OF READING DATA') END C END OF SUBROUTINE *** DHRW *** C C BEGIN OF SUBROUTINE *** DHSORT *** SUBROUTINE DHSORT(X,Y,N) C C ################################################################## C C BEGONNEN AM: 09.12.96 C LETZTE AENDERUNG: 04.01.98 C C PROGRAMMED BY: C THOMAS WIEDER C E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE C C PURPOSE: C SORTIEREN VON X-Y-DATEN NACH AUFSTEIGENDEN X-WERTEN. C C ################################################################## C C IMPLICIT NONE C C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. DOUBLE PRECISION X(N),Y(N) C .. C .. Local Scalars .. DOUBLE PRECISION XHILF,YHILF INTEGER I,J C .. DO 20 I = 1,N DO 10 J = I + 1,N IF (X(J).LT.X(I)) THEN C DIE BEIDEN DATENPAARE AUSTAUSCHEN: XHILF = X(I) YHILF = Y(I) X(I) = X(J) Y(I) = Y(J) X(J) = XHILF Y(J) = YHILF END IF 10 CONTINUE 20 CONTINUE C RETURN END C END OF SUBROUTINE *** DHSORT *** C C BEGIN OF SUBROUTINE *** DIVDIF *** SUBROUTINE DIVDIF(X,Y,N,M,IDEG,TABLE,XARG,VALUE,TRUBLE) C SUBROUTINE DIVDIF (X,Y,N,M,IDEG,XARG,VALUE,TRUBLE) C C FIRST VERSION: 24.03.1988 C PREVIOUS VERSION: 13.02.1996 C LATEST VERSION: 15.04.1998 C C SOURCE: ??? C I AM TRYING TO FIND THE SOURCE OF SUBROUTINE **** DIVDIF ***. C IN THE MEANTIME, I ASK THE AUTHOR(S) FOR EXCUSE. C I THINK, **** SUBROUTINE DIVDIF *** WAS GIVEN IN A BOOK C ON FORTRAN PROGRAMMING WHICH I READ DURING MY TIME AS PHD STUDENT. C C PURPOSE: C INTERPOLATION C C METHOD: C NEWTON'S DIVIDED DIFFERENCES C C IMPLICIT NONE C .. Scalar Arguments .. DOUBLE PRECISION VALUE,XARG INTEGER IDEG,M,N,TRUBLE C .. C .. Array Arguments .. DOUBLE PRECISION TABLE(N,N),X(N),Y(N) C .. C .. Local Scalars .. DOUBLE PRECISION YEST INTEGER I,IDEGM1,ISUB,ISUB1,ISUB2,J,MAX,NM1 C .. TRUBLE = 0 C CHECK FOR ARGUMENT CONSISTENCY IF (M.LT.N) GO TO 10 TRUBLE = 1 RETURN C CALCULATE FIRST-ORDER DIFFERENCES 10 NM1 = N - 1 DO 20 I = 1,NM1 TABLE(I,1) = (Y(I+1)-Y(I))/ (X(I+1)-X(I)) 20 CONTINUE IF (M.LE.1) GO TO 50 C CALCULATE HIGHER ORDER DIFFERENCES DO 40 J = 2,M DO 30 I = J,NM1 ISUB = I + 1 - J TABLE(I,J) = (TABLE(I,J-1)-TABLE(I-1,J-1))/ + (X(I+1)-X(ISUB)) 30 CONTINUE 40 CONTINUE 50 TRUBLE = 0 C CHECK FOR ARGUMENT CONSISTENCY IF (IDEG.LE.M) GO TO 60 TRUBLE = 2 VALUE = 0.0D0 RETURN C SEARCH X VECTOR FOR ELEMENT .GE. XARG 60 DO 70 I = 1,N IF (I.EQ.N .OR. XARG.LE.X(I)) GO TO 80 70 CONTINUE 80 MAX = I + IDEG/2 C INSURE THAT ALL REQUIRED DIFFERENCES ARE IN TABLE IF (MAX.LE.IDEG) MAX = IDEG + 1 IF (MAX.GT.N) MAX = N C COMPUTE INTERPOLANT VALUE YEST = TABLE(MAX-1,IDEG) IF (IDEG.LE.1) GO TO 100 IDEGM1 = IDEG - 1 DO 90 I = 1,IDEGM1 ISUB1 = MAX - I ISUB2 = IDEG - I YEST = YEST* (XARG-X(ISUB1)) + TABLE(ISUB1-1,ISUB2) 90 CONTINUE 100 ISUB1 = MAX - IDEG TRUBLE = 0 VALUE = YEST* (XARG-X(ISUB1)) + Y(ISUB1) RETURN END C END OF SUBROUTINE *** DIVDIF *** C C BEGIN OF SUBROUTINE *** FUNCHE *** SUBROUTINE FUNCHE(IERR7) C C ################################################################## C C BEGONNEN AM: 28.01.98 C VORLETZTE AENDERUNG: 13.10.98 C LATEST VERSION: 03.02.99 C C PROGRAMMED BY: C THOMAS WIEDER C E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE C C PURPOSE: C SUBROUTINE *** FUNCHE *** CALCULATES FUNCTION *** FUN *** C FOR VARIUS VALUES OF ARGUMENT *** X *** C IN ORDER TO LEARN ABOUT *** FUN ***. C C CALLING SEQUENCE: C CALLED BY SUBROUTINE *** HANKEL ***. C CALLS SUBROUTINE *** FUNK ***. C C MOST IMPORTANT VARIABLES: C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C X IS THE INPUT VARIABLE TO SUBROUTINE *** FUNK ***. C X IS A REAL-VALUED VARIBALE. C X IS THE INDEPENDENT VARIBALE OF FUNCTION *** FUN ***. C X >= 0.0D0 C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Y IS THE OUTPUT VARIABLE FROM SUBROUTINE *** FUNK ***. C Y IS THE DEPENDENT VARIBALE OF FUNCTION *** FUN ***. C Y IS A REAL-VALUED VARIBALE. C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C XMAX = LARGEST POSITVE MAGNITUDE OF DOUBLE PRECISION NUMBERS C PROVDIED BY THE COMPUTER C XTINY = SMALLEST POSITIVE MAGNITUDE OF DOUBLE PRECISION NUMBERS C PROVDIED BY THE COMPUTER C C ################################################################## C C IMPLICIT NONE C .. Parameters .. DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) C .. C .. Scalar Arguments .. INTEGER IERR7 C .. C .. Scalars in Common .. INTEGER LUERR,LUIN,LUOUT C .. C .. Local Scalars .. DOUBLE PRECISION NU,X,XLIMIT,XMAX,Y INTEGER IERR5,N C .. C .. Local Arrays .. DOUBLE PRECISION NUDAT(4) C .. C .. External Functions .. DOUBLE PRECISION D1MACH EXTERNAL D1MACH C .. C .. External Subroutines .. EXTERNAL DHFUNC C .. C .. Common blocks .. COMMON /INOUTE/LUIN,LUOUT,LUERR C .. C .. Data statements .. DATA NUDAT/-0.5D0,0.0D0,3.5D0,10.0D0/ C .. C C SET SOME PARAMETERS: C IERR7 = 0 C XTINY=D1MACH(1) XMAX = D1MACH(2) C C DO THE CALCLULATION: C C WE SELECT SOME VALUES FOR *** NU *** DO 10 N = 1,4 NU = NUDAT(N) C C CHECK *** FUN *** FOR *** X *** = 0. C X = ZERO Y = ZERO C IERR5 = 0 CALL DHFUNC(NU,X,Y,XLIMIT,IERR5) C IF (IERR5.NE.0) THEN WRITE (LUERR,FMT=9000) IERR5 IERR7 = 1 END IF IF (Y.LT.-XMAX) THEN WRITE (LUERR,FMT=9010) X,Y WRITE (LUERR,FMT=9020) END IF IF (Y.GT.XMAX) THEN WRITE (LUERR,FMT=9030) X,Y WRITE (LUERR,FMT=9020) END IF C C CHECK *** FUN *** FOR *** X *** = *** XMAX ***. C X = XMAX Y = ZERO C CALL DHFUNC(NU,X,Y,XLIMIT,IERR5) C IF (IERR5.NE.0) THEN WRITE (LUERR,FMT=9000) IERR5 IERR7 = 1 END IF IF (Y.LT.-XMAX) THEN WRITE (LUERR,FMT=9040) X,Y WRITE (LUERR,FMT=9020) END IF IF (Y.GT.XMAX) THEN WRITE (LUERR,FMT=9030) X,Y WRITE (LUERR,FMT=9020) END IF C IF (IERR7.NE.0) THEN C AN ERROR DURING CALL OF FUNCTION *** DHFUNC *** HAS OCCURRED. WRITE (LUERR,FMT=9050) NU,X,Y END IF C 10 CONTINUE C RETURN 9000 FORMAT (/,' MESSAGE FUNCHE:',/, + ' AN ERROR OCCURED DURING CALL TO *** DHFUNC ***!',/, + ' ERROR FLAG FROM *** DHFUNC *** : *** IERR *** =,',I6) 9010 FORMAT (/,' MESSAGE FUNCHE:',/,' Y = FUN(X=',D19.7,') = ',D19.7,/, + ' Y TOO SMALL, OVERFLOW!') 9020 FORMAT (/,' MESSAGE FUNCHE:',/, + ' THE HANKEL TRANSFORM POSSIBLY WI LL FAIL!',/, + ' PLEASE CHECK YOUR FUNCTION *** FUN *** !') 9030 FORMAT (/,' MESSAGE FUNCHE:',/,' Y = FUN(X=',D19.7,') = ',D19.7,/, + ' Y TOO BIG, OVERFLOW!') 9040 FORMAT (/,' MESSAGE FUNCHE:',/,' Y = FUN(X=',D19.7,') = ',D19.7,/, + ' Y TOO SMALL, OVERFLOW!') 9050 FORMAT (/,' MESSAGE FUNCHE:',/, + ' AN ERROR DURING CALL OF SUBROUTINE *** FUNC ***',/, + ' HAS BEEN REPORTED:',/,' NU, X, Y =',I6,5X,D17.9,5X,D17.9) END C END OF SUBROUTINE *** FUNCHE *** C C BEGIN OF SUBROUTINE *** ERRMSS *** SUBROUTINE ERRMSS(IERRME) C C ################################################################## C C FIRST VERSION: 14.01.1998 C PREVIOUS VERSION: 14.02.1998 C LATEST VERSION: 13.10.1998 C C PROGRAMMED BY: C THOMAS WIEDER C E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE C C PURPOSE: C OUTPUTS ERROR MESSAGES. C C MOST IMPORTANT VARIABLES: C C IERRME = FLAG FOR ERROR MESSAGES C C ################################################################## C C IMPLICIT NONE C C .. Scalar Arguments .. INTEGER IERRME C .. C .. Scalars in Common .. INTEGER LUERR,LUIN,LUOUT CHARACTER*80 FILEER,FILEIN,FILEOU C .. C .. Common blocks .. COMMON /INOUTE/LUIN,LUOUT,LUERR COMMON /NAMEN/FILEIN,FILEOU,FILEER integer i1mach, nout nout = i1mach(2) C .. WRITE(NOUT,FMT=9000) FILEER WRITE (LUERR,FMT=9000) FILEER WRITE (LUERR,FMT=9010) IERRME C IF (IERRME.EQ.-10) THEN WRITE (LUERR,FMT=9020) ELSE IF (IERRME.EQ.10) THEN WRITE (LUERR,FMT=9030) WRITE (LUERR,FMT=9040) ELSE IF (IERRME.EQ.11) THEN WRITE (LUERR,FMT=9050) WRITE (LUERR,FMT=9060) ELSE WRITE(NOUT,FMT=9070) IERRME WRITE (LUERR,FMT=9070) IERRME END IF C RETURN 9000 FORMAT (/,' MESSAGE ERRMSS:',/, + ' AN ERROR OCCURED DURING HANKEL TRANSFORMATION!!!',/, + ' RESULTS WILL BE UNRELIABLE OR WRONG!!!',/, + ' PLEASE LOOK INTO FILE:',/,A80) 9010 FORMAT (/,' MESSAGE ERRMS: ERROR FLAG *** IERRME *** =',I6) 9020 FORMAT (/,' MESSAGE ERRMS: PLEASE CHECK FILE NAMES, PATH ETC.!') 9030 FORMAT (/, + ' MESSAGE ERRMS: THE INTEGRAND COULD NOT BE CALCULATED!!!') 9040 FORMAT (' PLEASE CHECK YOUR FUNCTION *** FUN ***!',/, + 'ARE THERE ANY SINGULARITIES?',/, + 'DOES THE FUNCTION GO TO INFINITY (FOR LARGE ARGUMENTS)?') 9050 FORMAT (/,' MESSAGE ERRMS: THE INTEGRATION COULD NOT BE DONE!!!') 9060 FORMAT (' THE FUNCTION *** FUN *** MUST GO TO ZERO ',/, + ' FOR LA RGE ARGUMENTS!') 9070 FORMAT (/,' MESSAGE ERRMSS:',/, + ' UNKNOWN ERROR IDENTIFIER *** IERRMESS *** =',I6) END C END OF SUBROUTINE *** ERRMSS *** C C BEGIN OF SUBROUTINE *** DHJNUX *** SUBROUTINE DHJNUX(XIN,NUIN,IERR,JNUXV) C ################################################################## C C FIRST VERSION: 22.01.98 C PREVIOUS VERSION: 04.02.99 C LATEST VERSION: 04.02.99 C C PROGRAMMED BY: C THOMAS WIEDER C E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE C C PURPOSE: C C SUBROUTINE *** DHJNUX *** CALCULATES THE BESSEL FUNCTION C OF ORDER *** NU *** FOR ARGUMENT *** X ***. C C CALLING SEQUENCE: C C CALLED BY SUBROUTINE *** FUNINT ***. C CALLS SUBROUTINE *** RJBESL ***. C CALLS SUBROUTINE *** BSSJYA ***. C C MOST IMPORTANT VARIABLES: C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C X IS AN INPUT VARIABLE TO FUNCTION *** DHJNUX ***. C X IS A REAL-VALUED VARIBALE! C X >= 0.0D0 C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C NU IS AN INPUT VARIABLE TO FUNCTION *** JNUX ***. C NU IS A REAL-VALUED VARIBALE! C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C JNUXV IS AN OUTPUT VARIABLE OF FUNCTION *** JNUX ***. C JNUXV IS A REAL-VALUED VARIBALE! C JNUXV IS THE RESULT: JNUXV=JNUX(X) C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C ################################################################## C C IMPLICIT NONE C C SET SOME PARAMETERS: C C .. Parameters .. DOUBLE PRECISION PI,WINZIG,ZERO PARAMETER (PI=3.141592653D0,WINZIG=1.0D-09,ZERO=0.0D0) INTEGER NUMIMA PARAMETER (NUMIMA=101) C .. C .. Scalar Arguments .. DOUBLE PRECISION JNUXV,NUIN,XIN INTEGER IERR C .. C .. Scalars in Common .. INTEGER LUERR,LUIN,LUOUT C .. C .. Local Scalars .. DOUBLE PRECISION ALPHA,ARG,NU,NUMAX,NUMIN,X,XASYMP,XMAX,XMIN,YNUXV INTEGER NB,NCALC C .. C .. Local Arrays .. DOUBLE PRECISION RESULT(NUMIMA) C .. C .. External Functions .. DOUBLE PRECISION D1MACH EXTERNAL D1MACH C .. C .. External Subroutines .. EXTERNAL BSSJYA,RJBESL,RYBESL C .. C .. Intrinsic Functions .. INTRINSIC DBLE,DCOS,DSIN,IDINT C .. C .. Common blocks .. COMMON /INOUTE/LUIN,LUOUT,LUERR integer i1mach, nout nout = i1mach(2) C .. IERR = 0 C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C *** 24.02.1998 *** A GOOD VALUE FOR *** XASYMP *** = 100.0 XASYMP = 100.0D0 C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! NUMIN = -DBLE(NUMIMA) NUMAX = DBLE(NUMIMA) C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! XMIN = ZERO XMAX = D1MACH(2) C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C CHECK INPUT: C X = XIN NU = NUIN C C IS *** X *** POSITIVE? IF (X.LT.ZERO) THEN WRITE(NOUT,FMT=9000) WRITE (LUERR,FMT=9000) WRITE(NOUT,FMT=9010) X WRITE (LUERR,FMT=9010) X IERR = 1 GO TO 10 END IF C C CHECK ABSOLUT VALUE OF *** X *** IF (X.LT.XMIN .OR. X.GT.XMAX) THEN WRITE(NOUT,FMT=9000) WRITE (LUERR,FMT=9000) WRITE(NOUT,FMT=9020) X WRITE (LUERR,FMT=9020) X WRITE (LUERR,FMT=9030) XMIN,XMAX IERR = 1 GO TO 10 END IF C C CHECK ABSOLUT VALUE OF *** NU *** IF (NU.LT.NUMIN .OR. NU.GT.NUMAX) THEN WRITE(NOUT,FMT=9000) WRITE (LUERR,FMT=9000) WRITE(NOUT,FMT=9040) NU WRITE (LUERR,FMT=9040) NU WRITE (LUERR,FMT=9050) NUMIN,NUMAX IERR = 1 GO TO 10 END IF C C DO THE CALCULATION: C C ############################################################### C X = 0.0 C ############################################################### C C THE BESSEL FUNCTION IS DEFINED FOR X=0.0 ALSO: IF (X.EQ.ZERO) THEN IF (NU.EQ.ZERO) JNUXV = 1.0D0 IF (NU.GT.ZERO) JNUXV = 0.0D0 IF (NU.LT.ZERO) JNUXV = 0.0D0 RETURN END IF C C ############################################################### C X > 0.0 C ############################################################### C IF (NU.GE.0) THEN C IF (X.LE.XASYMP) THEN C C PREVENT NUMERICAL INSTABILITIES FOR VERY SMALL ARGUMENTS: IF (X.LT.WINZIG) THEN IF (NU.EQ.ZERO) JNUXV = 1.0D0 IF (NU.GT.ZERO) JNUXV = 0.0D0 RETURN END IF C ALPHA = NU - IDINT(NU) NB = IDINT(NU) + 1 C CALL RJBESL(X,ALPHA,NB,RESULT,NCALC) C IF (NCALC.NE.NB) THEN WRITE (LUERR,FMT=9060) WRITE (LUERR,FMT=9070) NCALC WRITE (LUERR,FMT=9080) X,NB,RESULT(NB) C IERR=1 C WRITE (LUERR,*) 'IERR =,IERR END IF JNUXV = RESULT(NB) C ELSE IF (X.GT.XASYMP) THEN C USE AN ASYMPTOTIC FORMULA FOR THE BESSEL FUNCTION: C CALL BSSJYA(X,NU,JNUXV) C C END OF IF-CLAUSE *** X *** END IF C C ################################################################# C NU < 0 C ################################################################# C ELSE IF (NU.LT.0.0D0) THEN C IF (X.LE.XASYMP) THEN C C PREVENT NUMERICAL INSTABILITIES FOR VERY SMALL ARGUMENTS: IF (X.LT.WINZIG) THEN JNUXV = 0.0D0 RETURN END IF C NU = -NU ALPHA = NU - IDINT(NU) NB = IDINT(NU) + 1 C CALL RJBESL(X,ALPHA,NB,RESULT,NCALC) C IF (NCALC.NE.NB) THEN WRITE (LUERR,FMT=9060) WRITE (LUERR,FMT=9070) NCALC WRITE (LUERR,FMT=9080) X,NB,RESULT(NB) IERR = 2 WRITE (LUERR,FMT=*) 'IERR =',IERR END IF NU = -NU JNUXV = RESULT(NB) C ELSE IF (X.GT.XASYMP) THEN C USE AN ASYMPTOTIC FORMULA FOR THE BESSEL FUNCTION: C CALL BSSJYA(X,-NU,JNUXV) C C END OF IF-CLAUSE *** X *** END IF C C CALCULATE THE BESSEL FUNCTION OF THE SECOND KIND Y(NU,X): C NU = -NU ALPHA = NU - IDINT(NU) NB = IDINT(NU) + 1 C C TO CIRCUMVENT PROBLEMS WITH SMALL ARGUMENTS *** X *** C IN CASE OF NEGATIVE *** NU ***: C IF (X.LT.1.0D-03) THEN C YNUXV=-1.0D0+X YNUXV = ZERO C ELSE C CALL RYBESL(X,ALPHA,NB,RESULT,NCALC) C IF (NCALC.LT.0) THEN WRITE (LUERR,FMT=9090) WRITE (LUERR,FMT=9100) NCALC WRITE (LUERR,FMT=9080) X,NB,RESULT(NB) IERR = 3 WRITE (LUERR,FMT=*) 'IERR =',IERR END IF NU = -NU YNUXV = RESULT(NB) C C END OF IF-CLAUSE *** X *** END IF C C FOR NEGATIVE *** NU *** SEE: C G. KORN AND T. KORN, MC GRAW HILL, 1968, CHAPTER 21.8-4. C W.H. PRESS, S.A. TEUKOLSKY, W.T. VETTERLING, B.P. FLANNERY: C NUMERICAL RECIPES, SECOND EDITION 1994, CAMBRDIGE UNIVERSITY PRESS C EQUATION (6.7.19): C C ARG=DMOD(-NU*PI,2.0D0*PI) ARG = -NU*PI JNUXV = DCOS(ARG)*JNUXV - DSIN(ARG)*YNUXV C C END OF IF-CLAUSE *** NU.LT.0.0D0 ***. END IF C 10 RETURN 9000 FORMAT (/,' MESSAGE JNUX: WRONG INPUT VALUE!') 9010 FORMAT (/,' MESSAGE JNUX:',/,' PARAMETER *** X *** IS NEGATIVE!', + /,' *** X *** =',D17.9) 9020 FORMAT (/,' MESSAGE JNUX:',/, + ' PARAMETER *** X *** OUTSIDE ALLOWED RANGE!',/, + ' *** X *** =',D17.9) 9030 FORMAT (/,' ALLOWED INTERVAL [X_MIN,X_MAX]:',/,'[',D17.9,',', + D17.9,']') 9040 FORMAT (/,' MESSAGE HANKEL:',/, + ' PARAMETER *** NU *** OUTSIDE ALLOWED RANGE!',/, + ' *** NU *** =',D17.9) 9050 FORMAT (/,' ALLOWED INTERVAL [NU_MIN,NU_MAX]:',/,'[',D17.9,',', + D17.9,']') 9060 FORMAT (/,'MESSAGE JNUX: NCALC NOT EQUAL NB! WARNING!',/, + ' THE ACCURACY OF THE RESULTS MAY BE LOW!') 9070 FORMAT (/,' MESSAGE JNUX: ',/, + ' NCALC ON RETURN FORM SUBROUTINE *** RJBESL *** =',I6) 9080 FORMAT (1X,'X, NB, RESULT(NB):',D17.9,3X,I6,3X,D17.9) 9090 FORMAT (/,'MESSAGE JNUX: NCALC NOT EQUAL NB! ERROR!') 9100 FORMAT (/,' MESSAGE JNUX: ',/, + ' NCALC ON RETURN FORM SUBROUTINE *** RYBESL *** =',I6) END C END OF SUBROUTINE *** DHJNUX *** C C BEGIN OF SUBROUTINE *** BSSJYA *** SUBROUTINE BSSJYA(X,NU,RESULT) C C ################################################################## C C FIRST VERSION: 12.11.1997 C PREVIOUS VERSION: 29.01.1998 C LATEST VERSION: 03.02.1998 C C PROGRAMMED BY: C THOMAS WIEDER C E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE C C THIS PROGRAM CALCULATES THE ASYMPTOTIC EXPANSION OF C THE FIRST BESSEL FUNCTION J(NU,X) FOR LARGE ARGUMENT X. C C LITERATURE: C C 1.) E.T. WHITTAKER AND G.N. WATSON: C ACOURSE OF MODERN ANALYSIS, C CAMBRIDGE UNIVERSITY PRESS, LONDON, C 4TH EDITION, REPRINT 1978, C P. 368. C 2.) R. COURANT AND D. HILBERT: C METHODEN DER MATHEMATISCHEN PHYSIK, C SPRINGER VERLAG, BERLIN, 3TE AUFLAGE, 1968, C KAPITEL VII. C C ################################################################## C C IMPLICIT NONE C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C *** 24.02.1998 *** A GOOD VALUE FOR *** XASYMP *** = 100.0 C .. Parameters .. DOUBLE PRECISION PI,ZERO PARAMETER (PI=3.141592653D0,ZERO=0.0D0) C .. C .. Scalar Arguments .. DOUBLE PRECISION NU,RESULT,X C .. C .. Local Scalars .. DOUBLE PRECISION ARG,XASYMP C .. C .. Intrinsic Functions .. INTRINSIC DCOS,DSQRT integer i1mach, nout nout = i1mach(2) C .. XASYMP = 100.0D0 C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C IF (X.LT.XASYMP) THEN WRITE (68,FMT=9000) X WRITE(NOUT,FMT=9000) X STOP END IF C IF (NU.LT.ZERO) THEN WRITE (68,FMT=9010) NU WRITE(NOUT,FMT=9010) NU STOP END IF C ARG = X - (PI*NU/2.0D0) - PI/4.0D0 RESULT = DSQRT(2.0D0/ (PI*X))*DCOS(ARG) C RETURN 9000 FORMAT (/,'MESSAGE BSSJYA: ARGUMENT *** X *** TOO LOW:',D17.9) 9010 FORMAT (/,'MESSAGE BSSJYA: ARGUMENT *** NU *** NEGATIVE:',D17.9) END C END OF SUBROUTINE *** BSSJYA *** C C BEGIN OF FUNCTION *** TFUNC *** DOUBLE PRECISION FUNCTION TFUNC(XIN) C C ################################################################# C C FIRST VERSION: 05.02.1998 C PREVIOUS VERSION: 05.02.1998 C LATEST VERSION: 06.02.1998 C C PROGRAMMED BY: C THOMAS WIEDER C E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE C C PURPOSE: C TEST SUBROUTINE *** INTHP ***. C C ################################################################# C C IMPLICIT NONE C C .. Scalar Arguments .. DOUBLE PRECISION XIN C .. C .. Scalars in Common .. DOUBLE PRECISION EXACT,XSPLIT INTEGER IINT C .. C .. Local Scalars .. DOUBLE PRECISION JNUXV,NU,X INTEGER IERRJ C .. C .. External Subroutines .. EXTERNAL DHJNUX C .. C .. Intrinsic Functions .. INTRINSIC DSIN,EXP C .. C .. Common blocks .. COMMON /LOSUNG/XSPLIT,EXACT,IINT integer i1mach, nout nout = i1mach(2) C .. IERRJ = 0 X = XIN C IF (IINT.EQ.1) THEN TFUNC = EXP(-X) C SOLUTION: 1.0 EXACT = 1.0D0 C ELSE IF (IINT.EQ.2) THEN TFUNC = 1.0D0/ (X**2+1.0D0) C SOLUTION: PI/2 EXACT = 1.570796327D0 C ELSE IF (IINT.EQ.3) THEN TFUNC = 1.0D0/ (X**1.1D0+1.0D0) C SOLUTION: 10.13724986 EXACT = 10.13724986D0 C ELSE IF (IINT.EQ.4) THEN TFUNC = DSIN(X)/ (X+1.0D0) C SOLUTION:0.6214496243 EXACT = 0.6214496243D0 C ELSE IF (IINT.EQ.5) THEN NU = 1.0D0 C CALL DHJNUX(X,NU,IERRJ,JNUXV) C TFUNC = JNUXV*EXP(-X) IF (IERRJ.NE.0) THEN WRITE(NOUT,FMT=9000) IERRJ END IF C SOLUTION: 1/(2+sqrt(2)) = 0.292893218 EXACT = 0.292893218D0 C ELSE IF (IINT.EQ.6) THEN NU = 1.0D0 C CALL DHJNUX(X,NU,IERRJ,JNUXV) C TFUNC = JNUXV IF (IERRJ.NE.0) THEN WRITE(NOUT,FMT=9000) IERRJ END IF C SOLUTION: 1.0 EXACT = 1.0D0 C ELSE IF (IINT.EQ.7) THEN NU = 1.0D0 C CALL DHJNUX(X,NU,IERRJ,JNUXV) C TFUNC = (1.0D0/ (X+1.0D0))*JNUXV IF (IERRJ.NE.0) THEN WRITE(NOUT,FMT=9000) IERRJ END IF C SOLUTION: 0.3980927698 EXACT = 0.3980927698D0 C ELSE IF (IINT.EQ.8) THEN NU = 3.0D0 C CALL DHJNUX(X,NU,IERRJ,JNUXV) C TFUNC = (1.0D0/ (X+1.0D0))*JNUXV IF (IERRJ.NE.0) THEN WRITE(NOUT,FMT=9000) IERRJ END IF C SOLUTION: 0.2464041090 EXACT = 0.2464041090D0 C ELSE IF (IINT.EQ.9) THEN NU = 5.0D0 C CALL DHJNUX(X,NU,IERRJ,JNUXV) C TFUNC = (1.0D0/ (X+1.0D0))*JNUXV IF (IERRJ.NE.0) THEN WRITE(NOUT,FMT=9000) IERRJ END IF C SOLUTION: 0.1659094729 EXACT = 0.1659094729D0 C ELSE WRITE(NOUT,FMT=*) 'UNKNOWN INTEGRAND IDENTIFICATION NUMBER!' STOP END IF C RETURN 9000 FORMAT (/,' AN ERROR OCCURED DURING CALL OF *** DHJNUX ***',I6) END C END OF FUNCTION *** TFUNC *** C C BEGIN OF FUNCTION *** HFUN *** DOUBLE PRECISION FUNCTION HFUN(X) C C ################################################################# C C FIRST VERSION: 05.02.1998 C PREVIOUS VERSION: 15.04.1998 C LATEST VERSION: 08.10.1998 C C PROGRAMMED BY: C THOMAS WIEDER C E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE C C PURPOSE: C CALCULATE SUBROUTINE *** DHFUNC ***. C C CALLING SEQUENCE: C CALLED BY SUBROUTINE *** OSCINT ***. C C ################################################################# C C IMPLICIT NONE C C .. Parameters .. INTEGER NMAX PARAMETER (NMAX=1000) C .. C .. Scalar Arguments .. DOUBLE PRECISION X C .. C .. Scalars in Common .. INTEGER IDAT,IREAD C .. C .. Local Scalars .. DOUBLE PRECISION Y INTEGER IDUMMY,IERR4 C .. C .. Local Arrays .. DOUBLE PRECISION DUMMYA(NMAX) C .. C .. External Subroutines .. EXTERNAL FUNK1,FUNK2,FUNK3 C .. C .. Common blocks .. COMMON /LOGIK/IDAT,IREAD C .. IDUMMY = 12 DUMMYA(1) = 0.0D0 C IF (IDAT.EQ.1) THEN CALL FUNK1(X,Y,IERR4) ELSE IF (IDAT.EQ.2) THEN CALL FUNK2(X,Y,IERR4) ELSE IF (IDAT.EQ.3) THEN CALL FUNK3(X,Y,DUMMYA,DUMMYA,IDUMMY,IERR4) END IF C IF (IERR4.NE.0) THEN WRITE (68,FMT=9000) END IF C HFUN = Y C RETURN 9000 FORMAT (/,' MESSAGE HFUN:',/, + ' ERROR MESSAGE FROM SUBROUTINE *** DHFUNK *** ',/, + ' *** IERR4 *** = ',I6) END C END OF FUNCTION *** HFUN *** C C BEGIN OF FUNCTION *** GPER *** DOUBLE PRECISION FUNCTION GPER(X) C C ################################################################# C C FIRST VERSION: 05.02.1998 C PREVIOUS VERSION: 19.02.1998 C LATEST VERSION: 24.02.1998 C C PROGRAMMED BY: C THOMAS WIEDER C E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE C C PURPOSE: C CALCULATE THE OSCILLATING PART OF THE INTEGRAND. C C CALLING SEQUENCE: C CALLED BY SUBROUTINE *** OSCINT ***. C C ################################################################# C C IMPLICIT NONE C C .. Scalar Arguments .. DOUBLE PRECISION X C .. C .. Scalars in Common .. DOUBLE PRECISION NU,XI,XLIMIT C .. C .. Local Scalars .. DOUBLE PRECISION JNUXV,XIX INTEGER IERRJ C .. C .. External Subroutines .. EXTERNAL DHJNUX C .. C .. Intrinsic Functions .. INTRINSIC DSQRT C .. C .. Common blocks .. COMMON /ARGUME/XI,NU,XLIMIT integer i1mach, nout nout = i1mach(2) C .. IERRJ = 0 C XIX = XI*X C CALL DHJNUX(XIX,NU,IERRJ,JNUXV) C GPER = DSQRT(XIX)*JNUXV C IF (IERRJ.NE.0) THEN WRITE(NOUT,FMT=9000) IERRJ WRITE(NOUT,FMT=9010) X END IF C RETURN 9000 FORMAT (/,' MESSAGE GPER:',/, + ' AN ERROR OCCURED DURING CALL OF *** DHJNUX ***',/, + ' *** IERRJ *** =',I6) 9010 FORMAT (' ARGUMENT *** X *** =',D17.9) END C END OF FUNCTION *** GPER *** C C BEGIN OF SUBROUTINE *** DGAMMA *** DOUBLE PRECISION FUNCTION DGAMMA(XIN) C C ################################################################# C C FIRST VERSION: 05.02.1998 C PREVIOUS VERSION: 05.02.1998 C LATEST VERSION: 06.02.1998 C C PROGRAMMED BY: C THOMAS WIEDER C E-MAIL: WIEDER@HRZ.UNI-KASSEL.DE C C PURPOSE: C C SUPPORT TO FUNCTION *** DLGAMA ***. C C ################################################################# C C IMPLICIT NONE C .. Parameters .. DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) C .. C .. Scalar Arguments .. DOUBLE PRECISION XIN C .. C .. Local Scalars .. DOUBLE PRECISION X C .. C .. External Functions .. DOUBLE PRECISION DLGAMA EXTERNAL DLGAMA C .. C .. Intrinsic Functions .. INTRINSIC EXP integer i1mach, nout nout = i1mach(2) C .. X = XIN IF (X.GT.ZERO) THEN DGAMMA = DLGAMA(X) DGAMMA = EXP(DGAMMA) ELSE IF (X.EQ.ZERO) THEN DGAMMA = ZERO ELSE WRITE(NOUT,FMT=9000) X STOP END IF RETURN 9000 FORMAT (/,' MESSAGE DGAMMA: WRONG INPUT VALUE X =',D17.9) END SHAR_EOF fi # end of overwriting check cd .. cd .. cd .. # End of shell archive exit 0