chave算法代码.docVIP

cc 以下为chave算法 subroutine hankel(BESR,BESI) implicit double precision(a-h,o-z) parameter (pi=3.14159265358979323846D0,nterm=50) PARAMETER(E0=8.85D-12,V0=4*PI*1.0D-7) double precision karg,kern dimension karg(255,nterm),kern(510,nterm),nk(nterm) COMPLEX*16 K0_2,K1_2 COMMON /PKC/DQ,XX,YY,WU,K0_2,K1_2 COMMON /PKD/II common /test/ng,nf,ni common /besint/nk,np,nps,karg,kern external f1 data rerr,AERR/1.e-11,1.0E-12/ R=dabs(XX) new=1 call besaut(besr,besi,0,1,7,r,f1,rerr,aerr,1,new,ierr) end c******************************************************************** subroutine besaut(besr,besi,norder,nl,nu,r,funct,rerr,aerr, * npcs,new,ierr) implicit double precision (a-h,o-z) dimension xsum(1) parameter (nsum=0) common /test/ng,nf,ni external funct nf=0 if(nl.gt.nu)then besr=0. besi=0. ierr=1 return endif nw=max(new,1) call bestrn(besr,besi,norder,nl,r,funct,.1*rerr,.1*aerr, * npcs,xsum,nsum,nw,ierr) if((ierr.ne.0).and.(nl.eq.7))then ng=nl return else oldr=besr oldi=besi do 10 n=nl+1,nu call bestrn(besr,besi,norder,n,r,funct,.1*rerr,.1*aerr, * npcs,xsum,nsum,2,ierr) if((ierr.ne.0).and.(n.eq.7))then besr=oldr besi=oldi ng=n return elseif((abs(besr-oldr).le.rerr*abs(besr)+aerr).and. * (abs(besi-oldi).le.rerr*abs(besi)+aerr))then ng=n return else oldr=besr oldi=besi endif 10 continue endif ng=7 ierr=1 return end C******************************************************************** subroutine bestrn(besr,besi,norder,ng,r,funct,rerr,aerr, * npcs,xsum,nsum,new,ierr) parameter (nterm=50, nstop=100) implicit double precision (a-h,o-z) double precision karg,kern,lastr,lasti dimension karg(255,nterm),kern(510,nterm),sr(nstop),si(nstop), * nk(nterm),xsum(1) common /besint/ nk,np,nps,karg,ke