hankel变换fortran程序.docx

1、hankel变换fortran程序0，1阶hankel变换的程序,滤波算法2010年05月17日 星期一 12:51 P.M.c=C 说在前面：c 关于hankel变换的程序在国内外都有许多程序，但是可能是出于某某种原因或目的，c 如：没时间公开，没地方公开，保护自己的知识产权等。c 你很难从某人，某单位直接获得可用的，解释详细的源代码，这种代码保护行为并形成恶性循环。 c 严重了影响了国内科技的发展，使用多数初学者做了大量重复性的工作，不能站在巨人的肩膀上进行创新。c 浪费广大青少年的青春和国内资源。这认为这是国内科技落后于国外的主要原因。c 在这里这鄙视一下一部分人，鄙视“你”是大家的合法

2、的权利。c 这里我整理出了几个hankel变换，并有详细的使用说明，供广大科研工作者，特别是初级研究者使用。由于baidu空间不法贴太多的代码，如果你有需要，请与我邮件联系。c (1如涉及了你的知识产权问题，请与我联系。c (2如果你已经获得了源代码，请不要保护起来自己用，请分发给同样需要此代码的人。让更多的人中收受益。谢谢。c=c 一个用于测试0，1阶hankel变换的程序。c 采用两种变换: 滤波算法(迟滞褶积，相关褶积和高斯积分。c 其中使用了六个子函数，c 1. SUBROUTINE HY2FHT(滤波算法与高斯积分，c 2. SUBROUTINE HANKELS(滤波算法c 3. C

3、OMPLEX HANKELY(滤波算法120,140c 4. COMPLEX HANKELK(滤波算法61,121,241c 5. COMPLEX*16 FUNCTION HANKELD(滤波算法801，c 6. COMPLEX*16 FUNCTION HANKELZ(滤波算法283。cc 说明：(1提高计算速度;计算精度；c (2一次计算多个变换参数或多个相关核函数的变换值,c (3减少对核函数求值计算次数和调用Bessel函数的次数。c (4数值卷积方法，无需计算Bessel函数值，则比直接数值积分速度快。c (5高阶变换可通过关系式jn-1(x+jn+1(x=2*n*jn(x/x.拆分实

4、现。c=下面代码是用来测试6个hankel变换的驱动函数。平时使用时并不需要，你可以从中学习如果调用这6个hankel变换函数的语法。下面的代码并不全，只是说明部分，如有需要，请与我联系，特别说明是GP免费。如果你已获得，请不要独享，希望是全国需要的人都有一份。把这个科技阻碍因素消除掉！c=c-c 程序整理: yiling (Huang Ling Email: yilingr0c Jilin university china. 吉林大学 地球探测科学与技术学院(2010C=program test_hankelTransformc=c 语法：c 1. SUBROUTINE HY2FHT(B,N

5、B,NREL,TOL,NTOL,RERR,AERR,NORD,FUN1,REL1,c * ZWORK,FWORK,BWORK,NPCS, IPATH,MAXDIM,MAXREL, ZANS,NFUN1,IERRc 2. SUBROUTINE HANKELS(BMAX,NB,NREL,TOL,NTOL,NORD,FUN1,IJREL,ZWORK,c * ZANS,ARG,NOFUN1,IERRc 3. COMPLEX FUNCTION HANKELY(NORD,B,FUN,NF c 4. COMPLEX FUNCTION HANKELK(NORD,B,FUN,FK,NF c 5. COMPLEX

6、*16 FUNCTION HANKELD(N,B,FUN,TOL,NF,NEWc 6. COMPLEX*16 FUNCTION HANKELZ(N,B,FUN,TOL,NF,NEW c=c=c 注意：变换精度是与核函数有关的。c 由不同的核函数计算的滤波系数，则不能适用于所有核函数(精度不足。c=c=c % 变量类型声名。implicit nonec- 公共参数 声名c 变换结果complex*16 Z,Zr1,Zr2c 变换阶次NORD_H：0 or1, 调用核函数次NF, 运行结果错误信息IERRinteger NORD_H,NF,IERRc 测试的核函数编号J, 本程序中循环过程用的变量

7、K,Kbinteger J,K,Kbc 变换参数BI, 希望的误差限(公差tolerance factor: HY2FHT: 0.1d-14, HANKELD, HANKELZ: 0.1d-5 real*8 BI,TOL_Hc 解析变换结果值EXACT， 滤波计算结果值FILT，计算结果的绝对误差ABSERR， 计算结果的相对误差RELERR real*8 EXACT,FILT,ABSERR,RELERRc 相关褶积,第一次非相关new=1,作相关褶积new=0integer NEWc 核函数声名complex*16 C1,C2,C3,C4,C5external C1,C2,C3,C4,C5c

8、 临时character(len=6 : AAAc-c 用于循环不同的Binteger,parameter : NB_LIST=5real*8 :B_LIST(NB_LISTDATA B_LIST/0.001, 0.01, 0.5,1.0,10.0/c- HY2FHT 调用参数 声名 c % 两个常数。 如： PARAMETER (MAXDIM= 100,MAXREL= 3 integer,parameter : maxdim= 100,maxrel= 3c % 对应各个相关核函数的hankel变换阶次。 integer nord(maxrel,nfun1(2 c % 固定形式。 real*8

9、 b(maxdim,bwork(maxdim,9,TOLc 希望的绝对误差ABSERR， 相对误差RELERRreal*8 RERR,AERRc % 固定形式。 complex*16 zwork(801,maxrel,zans(maxdim,maxrel,fwork(maxrelc % 相似相关核函数声名。COMPLEX*16 REL11,REL31 external REL11,REL31c 连续达到公差限的次数。对于平滑非震荡核函数，ntol=1即可。其实情况可设ntol=2,3.防止在自适就卷积过程中过早地截断。 integer NTOLc 相似相关核函数个数NREL, 运行路径设置IP

10、ATH, 用于besax2。c NPCS依赖于ipath参数，表示部分积分的子块数。通常设置npcs=1。对于比较难的核函数，npcs1。integer NB,NREL,IPATH,NPCSc- HANKELD 调用参数 声名 COMPLEX*16 HANKELDc 用于相关核函数计算时的保存数据空间，核函数值FSAVE_D(801COMPLEX*16 FSAVE_D(801c 核函数自变量值GSAVE_D(801real*8 GSAVE_D(801c 核函数值个数NSAVE_Dinteger NSAVE_Dc 核函数值及自变量公共空间 COMMON/SAVE_D/FSAVE_D,GSAVE_

11、D,NSAVE_Dc- HANKELZ 调用参数 声名 COMPLEX*16 HANKELZc 用于相关核函数计算时的保存数据空间，核函数值FSAVE_D(283，核函数自变量值GSAVE_D(283COMPLEX*16 FSAVE_Z(283c 核函数自变量值GSAVE_Z(283real*8 GSAVE_Z(283c 核函数值个数NSAVE_Dinteger NSAVE_Z c 核函数值及自变量公共空间COMMON/SAVE_Z/FSAVE_Z,GSAVE_Z,NSAVE_Zc- HANKELS 调用参数 声名c 变换参数个数NB_S, 相关函数个数NREL_S, 达到容差的次数ntol_

12、S, 函数最终调用核函数求值的次数nofun1integer NB_S,NREL_S,ntol_S,nofun1c 变换参数最大值BMAX_S, 容差限tol_Smreal BMAX_S,tol_Sc 变换结果ZANS_S(:,:, 用于工作的数据空间ZWORK_SCOMPLEX,allocatable : ZANS_S(:,:,ZWORK_S(:,:c 变换参数(输出arg，rg(i,i=1,nb，且有arg(i+1/arg(i=exp(.2。 real,allocatable : arg(:c 各核函数及相关核函数的变换阶次nord_S, 各相关核函数对应核函数的关系(i,j。funk(g

13、=g*ijrel(1,k * fun1(g*ijrel(2,k。integer,allocatable: nord_S(:,ijrel(:,:c 核函数C3_Scomplex C5_S,C4_S,C3_S,C1_Sexternal C5_S,C4_S,C3_S,C1_Sc- HANELY 调用参数 声名real B_Ycomplex HANKELYexternal HANKELYc- HANELK 调用参数 声名c 电导率,磁导率,频率,gamma.real sigma,u0,f,picomplex k0real RRcreal B_Kcomplex HANKELK,TSTexternal H

14、ANKELK,TSTcommon/TST_K/pi,sigma,f,u0c- HANELK 参数设置sigma = 3.2pi = 3.1415926u0 = 4*pi*1.d-7f = 1.k0=Csqrt(CMPLX(0.0,2*pi*u0*f*sigma;c=hanke变换(滤波算法-HANKELS变换函数2010年05月17日 星期一 12:21 P.M.c=C= Hankel Transform =C=c 程序整理: yiling (Huang Ling Email: yilingr0c Jilin university china. 吉林大学 地球探测科学与技术学院(2010 C=

15、SUBROUTINE HANKELS(BMAX,NB,NREL,TOL,NTOL,NORD,FUN1,IJREL,ZWORK,* ZANS,ARG,NOFUN1,IERR C= INTEGER NB,NREL,NTOL,NORD(NREL,IJREL(2,NREL,NOFUN1,IERRREAL BMAX,TOL,ARG(NBCOMPLEX ZWORK(283,NREL,ZANS(NB,NREL C= COMPLEX C,CMAX,ZSUM,ZERO,FUN1INTEGER KEY(283DOUBLE PRECISION E,ER,Y1,Y,ABSCISDIMENSION T(2,TMAX(2

16、DIMENSION WT0(283,WT1(283C-WE DEFINE COMPLEX C,CMAX TO BE EQUIVALENT TO REAL ARRAYS T(2,C TMAX(2, RESPECTIVELY, FOR USE IN THE TRUNCATION CRITERION TESTS,C WHERE C IS ANY CONVOLUTION PRODUCT AND CMAX IS THE MAXIMUMC CONVOLVED PRODUCT IN THE FIXED ABSCISSA RANGE (SEE PARAMETER TOL.EQUIVALENCE (C,T(1,

17、(CMAX,TMAX(1DATA ZERO/(0.0,0.0/C-ABSCIS=BASE CONSTANT FOR FILTER ABSCISSA GENERATIONDATA ABSCIS/0.7358852661479794460D0/C-E=DEXP(.2D0, ER=1.0D0/E (ALSO USED IN ABSCISSA GENERATION DATA E,ER/1.221402758160169834D0,0.818730753077981859D0/C-WT0(I=J0 HANKEL TRANSFORM FILTER WEIGHTS FOR I=1,283 DATA* WT0

18、( 1,WT0( 2,WT0( 3,WT0( 4,* WT0( 5,WT0( 6,WT0( 7,WT0( 8, * WT0( 9,WT0( 10,WT0( 11,WT0( 12, * WT0( 13,WT0( 14,WT0( 15,WT0( 16, * WT0( 17,WT0( 18,WT0( 19,WT0( 20, * WT0( 21,WT0( 22,WT0( 23,WT0( 24, * WT0( 25,WT0( 26,WT0( 27,WT0( 28, * WT0( 29,WT0( 30,WT0( 31,WT0( 32, * WT0( 33,WT0( 34,WT0( 35,WT0( 36/

19、* 2.1969101E-11, 4.1201161E-09,-6.1322980E-09, 7.2479291E-09,*-7.9821627E-09, 8.5778983E-09,-9.1157294E-09, 9.6615250E-09,*-1.0207546E-08, 1.0796633E-08,-1.1393033E-08, 1.2049873E-08,*-1.2708789E-08, 1.3446466E-08,-1.4174300E-08, 1.5005577E-08,*-1.5807160E-08, 1.6747136E-08,-1.7625961E-08, 1.8693427

20、E-08,*-1.9650840E-08, 2.0869789E-08,-2.1903555E-08, 2.3305308E-08,*-2.4407377E-08, 2.6033678E-08,-2.7186773E-08, 2.9094334E-08,*-3.0266804E-08, 3.2534013E-08,-3.3672072E-08, 3.6408936E-08,*-3.7425022E-08, 4.0787921E-08,-4.1543242E-08, 4.5756842E-08/DATA* WT0( 37,WT0( 38,WT0( 39,WT0( 40,* WT0( 41,WT0

21、( 42,WT0( 43,WT0( 44, * WT0( 45,WT0( 46,WT0( 47,WT0( 48, * WT0( 49,WT0( 50,WT0( 51,WT0( 52, * WT0( 53,WT0( 54,WT0( 55,WT0( 56, * WT0( 57,WT0( 58,WT0( 59,WT0( 60, * WT0( 61,WT0( 62,WT0( 63,WT0( 64, * WT0( 65,WT0( 66,WT0( 67,WT0( 68, * WT0( 69,WT0( 70,WT0( 71,WT0( 72/ *-4.6035233E-08, 5.1425075E-08,-5

22、.0893896E-08, 5.7934897E-08,*-5.6086570E-08, 6.5475248E-08,-6.1539913E-08, 7.4301996E-08,*-6.7117043E-08, 8.4767837E-08,-7.2583120E-08, 9.7366568E-08,*-7.7553611E-08, 1.1279873E-07,-8.1416723E-08, 1.3206914E-07,*-8.3217217E-08, 1.5663185E-07,-8.1482581E-08, 1.8860593E-07,*-7.3963141E-08, 2.3109673E-

23、07,-5.7243707E-08, 2.8867452E-07,*-2.6163525E-08, 3.6808773E-07, 2.7049871E-08, 4.7932617E-07,* 1.1407365E-07, 6.3720626E-07, 2.5241961E-07, 8.6373487E-07,* 4.6831433E-07, 1.1916346E-06, 8.0099716E-07, 1.6696015E-06/DATA* WT0( 73,WT0( 74,WT0( 75,WT0( 76,* WT0( 77,WT0( 78,WT0( 79,WT0( 80, * WT0( 81,W

24、T0( 82,WT0( 83,WT0( 84,* WT0( 85,WT0( 86,WT0( 87,WT0( 88, * WT0( 89,WT0( 90,WT0( 91,WT0( 92, * WT0( 93,WT0( 94,WT0( 95,WT0( 96, * WT0( 97,WT0( 98,WT0( 99,WT0(100, * WT0(101,WT0(102,WT0(103,WT0(104, * WT0(105,WT0(106,WT0(107,WT0(108/ * 1.3091334E-06, 2.3701475E-06, 2.0803829E-06,3.4012978E-06,* 3.2456774E-06, 4.9240402E-06, 5.0005198E-06,7.1783540E-06,* 7.6367633E-06, 1.0522038E-05, 1.1590021E-