hankel变换fortran程序.docx
《hankel变换fortran程序.docx》由会员分享，可在线阅读，更多相关《hankel变换fortran程序.docx（54页珍藏版）》请在冰豆网上搜索。
hankel变换fortran程序
0，1阶hankel变换的程序,滤波算法
2010年05月17日星期一12:
51P.M.
c=======================================================================
C说在前面：
c 关于hankel变换的程序在国内外都有许多程序，但是可能是出于某某种原因或目的，
c 如：
没时间公开，没地方公开，保护自己的知识产权等。
。
。
c 你很难从某人，某单位直接获得可用的，解释详细的源代码，这种代码保护行为并形成恶性循环。
c严重了影响了国内科技的发展，使用多数初学者做了大量重复性的工作，不能站在巨人的肩膀上进行创新。
c浪费广大青少年的青春和国内资源。
这认为这是国内科技落后于国外的主要原因。
c在这里这鄙视一下一部分人，鄙视“你”是大家的合法的权利。
c这里我整理出了几个hankel变换，并有详细的使用说明，供广大科研工作者，特别是初级研究者使用。
由于baidu空间不法贴太多的代码，如果你有需要，请与我邮件联系。
c （1如涉及了你的知识产权问题，请与我联系。
c （2如果你已经获得了源代码，请不要保护起来自己用，请分发给同样需要此代码的人。
让更多的人中收受益。
谢谢。
c=======================================================================
c一个用于测试0，1阶hankel变换的程序。
c 采用两种变换:
滤波算法（迟滞褶积，相关褶积和高斯积分。
c 其中使用了六个子函数，
c1. SUBROUTINEHY2FHT（滤波算法与高斯积分，
c2. SUBROUTINEHANKELS（滤波算法
c3. COMPLEX HANKELY（滤波算法120,140
c4. COMPLEX HANKELK（滤波算法61,121,241
c5. COMPLEX*16FUNCTIONHANKELD（滤波算法801，
c6. COMPLEX*16FUNCTIONHANKELZ（滤波算法283。
c
c说明：
（1提高计算速度;计算精度；
c （2一次计算多个变换参数或多个相关核函数的变换值,
c （3减少对核函数求值计算次数和调用Bessel函数的次数。
c （4数值卷积方法，无需计算Bessel函数值，则比直接数值积分速度快。
c （5高阶变换可通过关系式jn-1（x+jn+1（x=2*n*jn（x/x.拆分实现。
c=======================================================================
下面代码是用来测试6个hankel变换的驱动函数。
平时使用时并不需要，你可以从中学习如果调用这6个hankel变换函数的语法。
下面的代码并不全，只是说明部分，如有需要，请与我联系，特别说明是GP免费。
如果你已获得，请不要独享，希望是全国需要的人都有一份。
把这个科技阻碍因素消除掉！
！
！
c=======================================================================
c-----------------------------------------------------------------------
c程序整理:
yiling（HuangLingEmail:
yilingr0@
c Jilinuniversitychina.吉林大学地球探测科学与技术学院（2010
C=====================================================================
programtest_hankelTransform
c=======================================================================
c语法：
c1. SUBROUTINEHY2FHT（B,NB,NREL,TOL,NTOL,RERR,AERR,NORD,FUN1,REL1,
c *ZWORK,FWORK,BWORK,NPCS,IPATH,MAXDIM,MAXREL,ZANS,NFUN1,IERR
c 2. SUBROUTINEHANKELS（BMAX,NB,NREL,TOL,NTOL,NORD,FUN1,IJREL,ZWORK,
c *ZANS,ARG,NOFUN1,IERR
c 3. COMPLEXFUNCTIONHANKELY（NORD,B,FUN,NF
c 4. COMPLEXFUNCTIONHANKELK（NORD,B,FUN,FK,NF
c 5. COMPLEX*16FUNCTIONHANKELD（N,B,FUN,TOL,NF,NEW
c 6. COMPLEX*16FUNCTIONHANKELZ（N,B,FUN,TOL,NF,NEW
c=======================================================================
c=======================================================================
c注意：
变换精度是与核函数有关的。
c 由不同的核函数计算的滤波系数，则不能适用于所有核函数（精度不足。
c=======================================================================
c=======================================================================
c%变量类型声名。
implicitnone
c---------------------------------------------------------------公共参数声名
c变换结果
complex*16Z,Zr1,Zr2
c变换阶次NORD_H：
0or1,调用核函数次NF,运行结果错误信息IERR
integerNORD_H,NF,IERR
c测试的核函数编号J,本程序中循环过程用的变量K,Kb
integerJ,K,Kb
c变换参数BI,希望的误差限（公差tolerancefactor:
HY2FHT:
0.1d-14,HANKELD,HANKELZ:
0.1d-5
real*8BI,TOL_H
c解析变换结果值EXACT，滤波计算结果值FILT，计算结果的绝对误差ABSERR，计算结果的相对误差RELERR
real*8EXACT,FILT,ABSERR,RELERR
c相关褶积,第一次非相关new=1,作相关褶积new=0
integerNEW
c核函数声名
complex*16C1,C2,C3,C4,C5
externalC1,C2,C3,C4,C5
c临时
character（len=6:
:
AAA
c--------------------
c用于循环不同的B
integer,parameter:
:
NB_LIST=5
real*8:
:
B_LIST（NB_LIST
DATAB_LIST/0.001,0.01,0.5,1.0,10.0/
c---------------------------------------------------------------HY2FHT调用参数声名
c%两个常数。
如：
PARAMETER（MAXDIM=100,MAXREL=3
integer,parameter:
:
maxdim=100,maxrel=3
c%对应各个相关核函数的hankel变换阶次。
integernord（maxrel,nfun1（2
c%固定形式。
real*8b（maxdim,bwork（maxdim,9,TOL
c希望的绝对误差ABSERR，相对误差RELERR
real*8RERR,AERR
c%固定形式。
complex*16zwork（801,maxrel,zans（maxdim,maxrel,fwork（maxrel
c%相似相关核函数声名。
COMPLEX*16REL11,REL31
externalREL11,REL31
c连续达到公差限的次数。
对于平滑非震荡核函数，ntol=1即可。
其实情况可设ntol=2,3.防止在自适就卷积过程中过早地截断。
integerNTOL
c相似相关核函数个数NREL,运行路径设置IPATH,用于besax2。
cNPCS依赖于ipath参数，表示部分积分的子块数。
通常设置npcs=1。
对于比较难的核函数，npcs>1。
integerNB,NREL,IPATH,NPCS
c---------------------------------------------------------------HANKELD调用参数声名
COMPLEX*16HANKELD
c用于相关核函数计算时的保存数据空间，核函数值FSAVE_D（801
COMPLEX*16FSAVE_D（801
c核函数自变量值GSAVE_D（801
real*8GSAVE_D（801
c核函数值个数NSAVE_D
integerNSAVE_D
c核函数值及自变量公共空间
COMMON/SAVE_D/FSAVE_D,GSAVE_D,NSAVE_D
c---------------------------------------------------------------HANKELZ调用参数声名
COMPLEX*16HANKELZ
c用于相关核函数计算时的保存数据空间，核函数值FSAVE_D（283，核函数自变量值GSAVE_D（283
COMPLEX*16FSAVE_Z（283
c核函数自变量值GSAVE_Z（283
real*8GSAVE_Z（283
c核函数值个数NSAVE_D
integerNSAVE_Z
c核函数值及自变量公共空间
COMMON/SAVE_Z/FSAVE_Z,GSAVE_Z,NSAVE_Z
c---------------------------------------------------------------HANKELS调用参数声名
c变换参数个数NB_S,相关函数个数NREL_S,达到容差的次数ntol_S,函数最终调用核函数求值的次数nofun1
integerNB_S,NREL_S,ntol_S,nofun1
c变换参数最大值BMAX_S,容差限tol_Sm
realBMAX_S,tol_S
c变换结果ZANS_S（:
:
用于工作的数据空间ZWORK_S
COMPLEX,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=g**ijrel（1,k*fun1（g**ijrel（2,k。
integer,allocatable:
:
nord_S（:
ijrel（:
:
c核函数C3_S
complexC5_S,C4_S,C3_S,C1_S
externalC5_S,C4_S,C3_S,C1_S
c---------------------------------------------------------------HANELY调用参数声名
realB_Y
complexHANKELY
externalHANKELY
c---------------------------------------------------------------HANELK调用参数声名
c电导率,磁导率,频率,gamma.
realsigma,u0,f,pi
complexk0
realRR
c
realB_K
complexHANKELK,TST
externalHANKELK,TST
common/TST_K/pi,sigma,f,u0
c---------------------------------------HANELK参数设置
sigma=3.2
pi=3.1415926
u0=4*pi*1.d-7
f=1.
k0=Csqrt（CMPLX（0.0,2*pi*u0*f*sigma;
c=======================================================================
hanke变换（滤波算法-HANKELS变换函数2010年05月17日星期一12:
21P.M.c=======================================================================
C===========================HankelTransform====================
C=======================================================================
c程序整理:
yiling（HuangLingEmail:
yilingr0@
cJilinuniversitychina.吉林大学地球探测科学与技术学院（2010
C=====================================================================
SUBROUTINEHANKELS（BMAX,NB,NREL,TOL,NTOL,NORD,FUN1,IJREL,ZWORK,
*ZANS,ARG,NOFUN1,IERR
C=======================================================================
INTEGERNB,NREL,NTOL,NORD（NREL,IJREL（2,NREL,NOFUN1,IERR
REALBMAX,TOL,ARG（NB
COMPLEXZWORK（283,NREL,ZANS（NB,NREL
C=======================================================================
COMPLEXC,CMAX,ZSUM,ZERO,FUN1
INTEGERKEY（283
DOUBLEPRECISIONE,ER,Y1,Y,ABSCIS
DIMENSIONT（2,TMAX（2
DIMENSIONWT0（283,WT1（283
C-----WEDEFINECOMPLEXC,CMAXTOBEEQUIVALENTTOREALARRAYST（2,
CTMAX（2,RESPECTIVELY,FORUSEINTHETRUNCATIONCRITERIONTESTS,
CWHERECISANYCONVOLUTIONPRODUCTANDCMAXISTHEMAXIMUM
CCONVOLVEDPRODUCTINTHEFIXEDABSCISSARANGE（SEEPARAMETERTOL.
EQUIVALENCE（C,T（1,（CMAX,TMAX（1
DATAZERO/（0.0,0.0/
C-----ABSCIS=BASECONSTANTFORFILTERABSCISSAGENERATION
DATAABSCIS/0.7358852661479794460D0/
C-----E=DEXP（.2D0,ER=1.0D0/E（ALSOUSEDINABSCISSAGENERATION
DATAE,ER/1.221402758160169834D0,0.818730753077981859D0/
C-----WT0（I=J0HANKELTRANSFORMFILTERWEIGHTSFORI=1,283
DATA
*WT0（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/
*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.8693427E-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（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.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-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,WT0（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-
升级会员 