************************************************************************ *** WQfunc, Wmfunc, Vmfunc *** WGfunc, Wmufunc, Vmufunc *** WGHFfunc, WGHFDfunc : WGfunc with spetial cut for hydro-foil *** WGNKfunc : WGfunc with spetial cut for semi-submerged body *** WGdfunc, VGdfunc : A pair of wave vortex *** WGintNK, VGintNK, LogNK, VLog *** = WGint, = VGint *** WGintHF, VGintHF, LogHF *** Expint *** StrmLnPlt, TOKOSN *** fixed on 2015/12/01 ************************************************************************ complex function WQfunc(z,zeta,K,ThetaCut,sWQ) * * ThetaCut : Angle of Cut of Log(z-zeta) must be POSITIVE * * ThetaCut - 2 * Pi < Argument of Log(z-zeta) =< ThetaCut * .lt. .le. * - 3/2 Pi < Argument of Log(z-zetaBar) & E1 =< Pi/2 * * 0 < K < 50 : Usual Wave Making Condition * K >= 50 : Limit Case of Low Speed ( Mirror Image Flow) * K = 0 : Limit Case of High Speed ( Inverse Mirror Image Flow) * K < 0 : Infinite Region Flow * real K complex i/(0,1)/ complex z,zeta,zetaBar,miKZ,E1 complex LogZ,LogZBar,EmiKZ,Skappa,sWQ pi=acos(0.0)*2 pi2=2*pi pih=pi/2 * LogZ=log(z-zeta) Alog=real(LogZ) ArgZ=imag(LogZ) if(ArgZ.le.ThetaCut-pi2) then ArgZ=ArgZ+pi2 elseif(ArgZ.gt.ThetaCut) then ArgZ=ArgZ-pi2 endif LogZ=cmplx(ALog,ArgZ) * zetaBar=conjg(zeta) LogZBar=log(z-zetaBar) Alog=real(LogZBar) ArgZ=imag(LogZBar) if(ArgZ.gt.pih) ArgZ=ArgZ-pi2 LogZBar=cmplx(ALog,ArgZ) * if(K.lt.0) then sWQ=0 WQfunc=LogZ elseif(K.eq.0) then sWQ = -LogZBar WQfunc=LogZ-LogZBar elseif(K.lt.50) then miKZ=-i*K*(z-zetaBar) EmiKZ=exp(miKZ) call expint(miKZ,E1) xsign=real(z-zeta) if(xsign.gt.0) then Skappa=EmiKZ*E1 else Skappa=EmiKZ*E1+pi2*i*EmiKZ endif sWQ = LogZBar+2*Skappa WQfunc=LogZ+LogZBar+2*Skappa else sWQ = LogZBar WQfunc=LogZ+LogZBar endif end ************************************************************************ complex function Wmfunc(z,zeta,K,sWm) * * 0 < K < 50 : Usual Wave Making Condition * K >= 50 : Limit Case of Low Speed ( Mirror Image Flow) * K = 0 : Limit Case of High Speed ( Inverse Mirror Image Flow) * K < 0 : Infinite Region Flow * real K complex i/(0,1)/ complex z,zeta,zetaBar,miKZ,E1 complex InvZ,InvZBar,EmiKZ,Skappa,sWm pi=acos(0.0)*2 pi2=2*pi * InvZ=1/(z-zeta) zetaBar=conjg(zeta) InvZBar=1/(z-zetaBar) * if(K.lt.0) then sWm=0 Wmfunc=-InvZ elseif(K.eq.0) then sWm = InvZBar Wmfunc=-InvZ+InvZBar elseif(K.lt.50) then miKZ=-i*K*(z-zetaBar) EmiKZ=exp(miKZ) call expint(miKZ,E1) xsign=real(z-zeta) if(xsign.ge.0) then Skappa=EmiKZ*E1 else Skappa=EmiKZ*E1+pi2*i*EmiKZ endif sWm = InvZBar+2*i*K*Skappa Wmfunc=-InvZ+InvZBar+2*i*K*Skappa else sWm = -InvZBar Wmfunc=-InvZ-InvZBar endif end ************************************************************************ complex function Vmfunc(z,zeta,K,sVm) * * 0 < K < 50 : Usual Wave Making Condition * K >= 50 : Limit Case of Low Speed ( Mirror Image Flow) * K = 0 : Limit Case of High Speed ( Inverse Mirror Image Flow) * K < 0 : Infinite Region Flow * real K complex i/(0,1)/ complex z,zeta,zetaBar,miKZ,E1 complex InvZ2,InvZBar,InvZBar2,EmiKZ,Skappa,sVm pi=acos(0.0)*2 pi2=2*pi * InvZ2= 1/(z-zeta)**2 zetaBar=conjg(zeta) InvZBar= 1/(z-zetaBar) InvZBar2=1/(z-zetaBar)**2 * miKZ=-i*K*(z-zetaBar) EmiKZ=exp(miKZ) call expint(miKZ,E1) xsign=real(z-zeta) if(xsign.ge.0) then Skappa=EmiKZ*E1 else Skappa=EmiKZ*E1+pi2*i*EmiKZ endif sVm= -InvZBar2-2*K*i*InvZBar+2*K**2*Skappa Vmfunc=InvZ2-InvZBar2-2*K*i*InvZBar+2*K**2*Skappa end ************************************************************************ complex function WGfunc(z,zeta,K,ThetaCut,sWG) * * 0 < K < 50 : Usual Wave Making Condition * K >= 50 : Limit Case of Low Speed ( Mirror Image Flow) * K = 0 : Limit Case of High Speed ( Inverse Mirror Image Flow) * K < 0 : Infinite Region Flow * * ThetaCut : Angle of Cut of Log(z-zeta) must be POSITIVE * * ThetaCut - 2 * Pi < Argument of Log(z-zeta) =< ThetaCut * .lt. .le. * - 3/2 Pi < Argument of Log(z-zetaBar) & E1 =< Pi/2 * * real K complex i/(0,1)/ complex z,zeta,zetaBar,miKZ,E1 complex LogZ,LogZBar,EmiKZ,Skappa,sWG pi=acos(0.0)*2 pi2=2*pi pih=pi/2 * LogZ=log(z-zeta) Alog=real(LogZ) ArgZ=imag(LogZ) if(ArgZ.le.ThetaCut-pi2) then ArgZ=ArgZ+pi2 elseif(ArgZ.gt.ThetaCut) then ArgZ=ArgZ-pi2 endif LogZ=cmplx(ALog,ArgZ) * zetaBar=conjg(zeta) LogZBar=log(z-zetaBar) Alog=real(LogZBar) ArgZ=imag(LogZBar) if(ArgZ.gt.pih) ArgZ=ArgZ-pi2 LogZBar=cmplx(ALog,ArgZ) * if(K.lt.0) then WGfunc=-i*LogZ sWG=0 elseif(K.eq.0) then sWG = -i*LogZBar WGfunc=-i*LogZ-i*LogZBar elseif(K.lt.50) then miKZ=-i*K*(z-zetaBar) EmiKZ=exp(miKZ) call expint(miKZ,E1) xsign=real(z-zeta) if(xsign.gt.0) then Skappa=EmiKZ*E1 else Skappa=EmiKZ*E1+pi2*i*EmiKZ endif sWG = i*LogZBar+2*i*Skappa WGfunc=-i*LogZ+i*LogZBar+2*i*Skappa else sWG = i*LogZBar WGfunc=-i*LogZ+i*LogZBar endif end ************************************************************************ complex function Wmufunc(z,zeta,K,sWmu) * * 0 < K < 50 : Usual Wave Making Condition * K >= 50 : Limit Case of Low Speed ( Mirror Image Flow) * K = 0 : Limit Case of High Speed ( Inverse Mirror Image Flow) * K < 0 : Infinite Region Flow * real K complex i/(0,1)/ complex z,zeta,zetaBar,miKZ,E1 complex InvZ,InvZBar,EmiKZ,Skappa,sWmu pi=acos(0.0)*2 pi2=2*pi * InvZ=1/(z-zeta) zetaBar=conjg(zeta) InvZBar=1/(z-zetaBar) * if(K.lt.0) then sWmu=0 Wmufunc=i*InvZ elseif(K.eq.0) then Wmufunc=i*InvZ+i*InvZBar sWmu= i*InvZBar elseif(K.lt.50) then miKZ=-i*K*(z-zetaBar) EmiKZ=exp(miKZ) call expint(miKZ,E1) xsign=real(z-zeta) if(xsign.gt.0) then Skappa=EmiKZ*E1 else Skappa=EmiKZ*E1+pi2*i*EmiKZ endif Wmufunc=i*InvZ+i*InvZBar-2*K*Skappa sWmu = i*InvZBar-2*K*Skappa else Wmufunc=i*InvZ-i*InvZBar sWmu = -i*InvZBar endif end ************************************************************************** complex function Vmufunc(z,zeta,K,sVmu) * * 0 < K < 50 : Usual Wave Making Condition * K >= 50 : Limit Case of Low Speed ( Mirror Image Flow) * K = 0 : Limit Case of High Speed ( Inverse Mirror Image Flow) * K < 0 : Infinite Region Flow * real K complex i/(0,1)/ complex z,zeta,zetaBar,miKZ,E1 complex InvZ2,InvZBar,InvZBar2,EmiKZ,Skappa,sVmu pi=acos(0.0)*2 pi2=2*pi * InvZ2= 1/(z-zeta)**2 zetaBar=conjg(zeta) InvZBar= 1/(z-zetaBar) InvZBar2=1/(z-zetaBar)**2 * miKZ=-i*K*(z-zetaBar) EmiKZ=exp(miKZ) call expint(miKZ,E1) xsign=real(z-zeta) if(xsign.ge.0) then Skappa=EmiKZ*E1 else Skappa=EmiKZ*E1+pi2*i*EmiKZ endif sVmu = -i*InvZBar2+2*K*InvZBar+2*i*K**2*Skappa Vmufunc=-i*InvZ2-i*InvZBar2+2*K*InvZBar+2*i*K**2*Skappa end ************************************************************************ complex function WGHFfunc(z,zeta,K,zetaE,sWG) * * Cut of Log(z-zeta) is set as * (-Inf,etaE) - zetaE - zeta * * 0 < K < 50 : Usual Wave Making Condition * K >= 50 : Limit Case of Low Speed ( Mirror Image Flow) * K = 0 : Limit Case of High Speed ( Inverse Mirror Image Flow) * K < 0 : Infinite Region Flow * real K complex i/(0,1)/ complex z,zeta,zetaE,zetaBar,miKZ,E1 complex LogZ,LogZBar,EmiKZ,Skappa,sWG pi=acos(0.0)*2 pi2=2*pi pih=pi/2 * LogZ=log(z-zeta) Alog=real(LogZ) ArgZ=imag(LogZ) x=real(z) y=imag(z) xsi=real(zeta) eta=imag(zeta) xsiE=real(zetaE) etaE=imag(zetaE) if(x.lt.xsiE) then if(y.ge.etaE) then if(ArgZ.lt.0) then ArgZ=ArgZ+pi2 endif else if(ArgZ.gt.0) then ArgZ=ArgZ-pi2 endif endif else if(x.lt.xsi) then ul=(y-etaE)*(xsi-xsiE)-(eta-etaE)*(x-xsiE) if(ul.ge.0) then if(ArgZ.lt.0) then ArgZ=ArgZ+pi2 endif else if(ArgZ.gt.0) then ArgZ=ArgZ-pi2 endif endif endif endif LogZ=cmplx(ALog,ArgZ) * zetaBar=conjg(zeta) LogZBar=log(z-zetaBar) Alog=real(LogZBar) ArgZ=imag(LogZBar) if(ArgZ.gt.pih) ArgZ=ArgZ-pi2 LogZBar=cmplx(ALog,ArgZ) if(K.lt.0) then WGHFfunc=-i*LogZ sWG=0 elseif(K.eq.0) then WGHFfunc=-i*LogZ-i*LogZBar sWG = -i*LogZBar elseif(K.lt.50) then miKZ=-i*K*(z-zetaBar) EmiKZ=exp(miKZ) call expint(miKZ,E1) xsign=real(z-zeta) if(xsign.gt.0) then Skappa=EmiKZ*E1 else Skappa=EmiKZ*E1+pi2*i*EmiKZ endif sWG = i*LogZBar+2*i*Skappa WGHFfunc=-i*LogZ+i*LogZBar+2*i*Skappa else sWG = i*LogZBar WGHFfunc=-i*LogZ+i*LogZBar endif end ************************************************************************ complex function WGHFDfunc(z,zeta,K,zetaM,zetaE,swG) * * Cut of Log(z-zeta) is set as * (-Inf,etaE) - zetaE - zetaM - zeta * * 0 < K < 50 : Usual Wave Making Condition * K >= 50 : Limit Case of Low Speed ( Mirror Image Flow) * K = 0 : Limit Case of High Speed ( Inverse Mirror Image Flow) * K < 0 : Infinite Region Flow * real K complex i/(0,1)/ complex z,zeta,zetaE,zetaM,zetaBar,miKZ,E1 complex LogHF,LogZ,LogZBar,EmiKZ,Skappa,swG pi=acos(0.0)*2 pi2=2*pi pih=pi/2 * LogZ=LogHF(z,zeta,zetaM,zetaE) * zetaBar=conjg(zeta) LogZBar=log(z-zetaBar) Alog=real(LogZBar) ArgZ=imag(LogZBar) if(ArgZ.gt.pih) ArgZ=ArgZ-pi2 LogZBar=cmplx(ALog,ArgZ) if(K.lt.0) then swG=0 WGHFDfunc=-i*LogZ elseif(K.eq.0) then swG = -i*LogZBar WGHFDfunc=-i*LogZ-i*LogZBar elseif(K.lt.50) then miKZ=-i*K*(z-zetaBar) EmiKZ=exp(miKZ) call expint(miKZ,E1) xsign=real(z-zeta) if(xsign.gt.0) then Skappa=EmiKZ*E1 else Skappa=EmiKZ*E1+pi2*i*EmiKZ endif swG = i*LogZBar+2*i*Skappa WGHFDfunc=-i*LogZ+i*LogZBar+2*i*Skappa else swG = i*LogZBar WGHFDfunc=-i*LogZ+i*LogZBar endif end ************************************************************************ complex function WGNKfunc(z,zeta,K,zetaM,sWG) * * Wave Vortex with cut of LogNK * * 0 < K < 50 : Usual Wave Making Condition * K >= 50 : Limit Case of Low Speed ( Mirror Image Flow) * K = 0 : Limit Case of High Speed ( Inverse Mirror Image Flow) * K < 0 : Infinite Region Flow * * Cut of log(z-zeta) is set as * (xsiM,inf) - zeta - zeta * real K complex i/(0,1)/ complex z,zeta,zetaM,zetaBar,miKZ,E1,LogNK complex LogZ,LogZBar,EmiKZ,SKappa,swG pi=acos(0.0)*2 pi2=2*pi pih=pi/2 * LogZ=LogNK(z,zeta,zetaM) * zetaBar=conjg(zeta) LogZBar=log(z-zetaBar) Alog=real(LogZBar) ArgZ=imag(LogZBar) if(ArgZ.gt.pih) ArgZ=ArgZ-pi2 LogZBar=cmplx(ALog,ArgZ) * if(K.lt.0) then swG = 0 WGNKfunc=-i*LogZ elseif(K.eq.0) then swG = -i*LogZBar WGNKfunc=-i*LogZ-i*LogZBar elseif(K.lt.50) then miKZ=-i*K*(z-zetaBar) EmiKZ=exp(miKZ) call expint(miKZ,E1) xsign=real(z-zeta) if(xsign.gt.0) then Skappa=EmiKZ*E1 else Skappa=EmiKZ*E1+pi2*i*EmiKZ endif swG = i*LogZBar+2*i*Skappa WGNKfunc=-i*LogZ+i*LogZBar+2*i*Skappa else swG=i*LogZBar WGNKfunc=-i*LogZ+swG endif end *********************************************************************** complex function WGdfunc(z,zeta2,zeta1,K) * * a pair of wave VORTEX WG(a,zeta2) - WG(z,zeta1) * * Cut : zeta1 -- zeta2 * real K complex WG1,WG2,swG1,swG2 complex WGNKfunc complex z,zeta2,zeta1 complex dz /(0,0.000001)/ * * for zeta1 WG1=WGNKfunc(z,zeta1,K,zeta2,swG1) * for zeta2 ** WG2=WGNKfunc(z,zeta2,K,zeta2,swG2) WG2=WGNKfunc(z,zeta2,K,zeta2+dz,swG2) * Difference WGdfunc=WG2-WG1 end *********************************************************************** complex function VGdfunc(z,zeta2,zeta1,K) * * Complex Velocity due to * a couple of wave VORTEX WG(a,zeta2) - WG(z,zeta1) * real K complex VG1,VG2 complex Wmufunc,Wmuf complex z,zeta2,zeta1 * * for zeta1 VG1=-Wmufunc(z,zeta1,K,Wmuf) * for zeta2 VG2=-Wmufunc(z,zeta2,K,Wmuf) * Difference VGdfunc=VG2-VG1 end ************************************************************************ complex function WGintNK(z,zeta2,zeta1,K,zetaM) * WGintNK * * The Revised Version * * For a Flow aroud a Water-Surface-Piercing-Body * Integrated Wave Vortex Segment * with zeta1, zeta2 and zetaM located clockwisely * Cut of log(z-zeta) is set as * zeta1,2 -- zetaM -- (xsiM,etaE) -- (-inf,etaE) * etaE > 10 * Cut of log(z-zetaBar) as * zeatBar -- (xsiBar,inf) * * Cut of E1 as * zeatBar -- (xsiBar,inf) * * 0 < K < 50 : Usual Wave Making Condition * K >= 50 : Limit Case of Low Speed ( Mirror Image Flow) * K = 0 : Limit Case of High Speed ( Inverse Mirror Image Flow) * K < 0 : Infinite Region Flow * real K complex i/(0,1)/ complex WGint1,WGint2,LogNK,VLog complex z,zeta1,zeta2,zetaM complex dzeta,zetaBar,zdif,zdifBar complex LogZ,LogZBar complex EiG,miKz,EmiKz,E1,SKappa pi=acos(0.0)*2 pi2=pi*2 * x=real(z) xsi1=real(zeta1) xsi2=real(zeta2) dzeta=zeta2-zeta1 * * For zeta2 * LogZ=LogNK(z,zeta2,zetaM) zdif=z-zeta2 * zetaBar=conjg(zeta2) zdifBar=z-zetaBar LogZBar=VLog(zdifBar) EiG=dzeta/abs(dzeta) if(K.lt.0) then WGint2=i*conjg(EiG)*zdif*(LogZ-1) elseif(K.eq.0) then WGint2=i*conjg(EiG)*zdif*(LogZ-1)+i*EiG*zdifBar*(LogZBar-1) elseif(K.lt.50) then miKz=-i*K*zdifBar EmiKz=exp(miKz) call expint(miKZ,E1) if(x-xsi2.gt.0) then SKappa=EmiKz*E1 else SKappa=EmiKz*E1+pi2*i*EmiKz endif WGint2=i*conjg(EiG)*zdif*(LogZ-1)-i*EiG*zdifBar*(LogZBar-1) * +2/K*EiG*LogZBar+2/K*EiG*SKappa else WGint2=i*conjg(EiG)*zdif*(LogZ-1)-i*EiG*zdifBar*(LogZBar-1) endif * * For zeta1 * LogZ=LogNK(z,zeta1,zetaM) zdif=z-zeta1 * zetaBar=conjg(zeta1) zdifBar=z-zetaBar LogZBar=VLog(zdifBar) ** EiG=dzeta/abs(dzeta) if(K.lt.0) then WGint1=i*conjg(EiG)*zdif*(LogZ-1) elseif(K.eq.0) then WGint1=i*conjg(EiG)*zdif*(LogZ-1)+i*EiG*zdifBar*(LogZBar-1) elseif(K.lt.50) then miKz=-i*K*zdifBar EmiKz=exp(miKz) call expint(miKZ,E1) if(x-xsi1.gt.0) then SKappa=EmiKz*E1 else SKappa=EmiKz*E1+pi2*i*EmiKz endif WGint1=i*conjg(EiG)*zdif*(LogZ-1)-i*EiG*zdifBar*(LogZBar-1) * +2/K*EiG*LogZBar+2/K*EiG*SKappa else WGint1=i*conjg(EiG)*zdif*(LogZ-1)-i*EiG*zdifBar*(LogZBar-1) endif WGintNK=WGint2-WGint1 end ************************************************************************ complex function VGintNK(z,zeta2,zeta1,K,zetaM) * * The Revised Version * * VGintNK * Points zeta2 and zeta1 are located clockwisely * Cut of log(z-zeta) is set as * zeta -- zetaM -- (xsiM,inf) * Cut of log(z-zetaBar) & E1 as * zetaBar -- (xsiBar,inf) * * 0 < K < 50 : Usual Wave Making Condition * K >= 50 : Limit Case of Low Speed ( Mirror Image Flow) * K = 0 : Limit Case of High Speed ( Inverse Mirror Image Flow) * K < 0 : Infinite Region Flow * real K complex i/(0,1)/ complex VGint1,VGint2,LogNK,VLog complex z,zeta1,zeta2,zetaM,zetaE complex dzeta,zetaBar,zdifBar complex LogZ,LogZBar complex EiG,miKz,EmiKz,E1,SKappa real etaE/10/ pi=acos(0.0)*2 pi2=pi*2 * x=real(z) xsi1=real(zeta1) xsi2=real(zeta2) dzeta=zeta2-zeta1 * For zetaE zetaE=conjg(zetaM)+i*etaE * * For zeta2 * LogZ=LogNK(z,zeta2,zetaM,zetaE) * zetaBar=conjg(zeta2) zdifBar=z-zetaBar LogZBar=VLog(zdifBar) EiG=dzeta/abs(dzeta) if(K.lt.0) then VGint2=i*conjg(EiG)*LogZ elseif(K.eq.0) then VGint2=i*conjg(EiG)*LogZ+i*EiG*LogZBar elseif(K.lt.50) then miKz=-i*K*zdifBar EmiKz=exp(miKz) call expint(miKZ,E1) if(x-xsi2.gt.0) then SKappa=EmiKz*E1 else SKappa=EmiKz*E1+pi2*i*EmiKz endif VGint2=i*conjg(EiG)*LogZ-i*EiG*LogZBar * -2*i*EiG*SKappa else VGint2=i*conjg(EiG)*LogZ-i*EiG*LogZBar endif * * For zeta1 * LogZ=LogNK(z,zeta1,zetaM,zetaE) * zetaBar=conjg(zeta1) zdifBar=z-zetaBar LogZBar=VLog(zdifBar) ** EiG=dzeta/abs(dzeta) if(K.lt.0) then VGint1=i*conjg(EiG)*LogZ elseif(K.eq.0) then VGint1=i*conjg(EiG)*LogZ+i*EiG*LogZBar elseif(K.lt.50) then miKz=-i*K*zdifBar EmiKz=exp(miKz) call expint(miKZ,E1) if(x-xsi1.gt.0) then SKappa=EmiKz*E1 else SKappa=EmiKz*E1+pi2*i*EmiKz endif VGint1=i*conjg(EiG)*LogZ-i*EiG*LogZBar * -2*i*EiG*SKappa else VGint1=i*conjg(EiG)*LogZ-i*EiG*LogZBar endif VGintNK=VGint2-VGint1 end *************************************************************************** complex function VLog(z) * * Complex Logarithmic Function for Iamg(z) > 0 * complex z,pi2i * pi=acos(0.0)*2 pi2i=cmplx(0,pi*2) pih=pi/2 * VLog=log(z) if(imag(VLog).gt.pih) VLog=VLog-pi2i end *************************************************************************** complex function LogNK(z,zeta,zetaM) * * Complex Logarithmic Function for Neumann-Kelvin Problem * ( See Section 3.3.4 ) * * Cut of Log ; zeta - zetaM - (xiM,inf) * complex z,zeta,zetaM,LogZ * pi=acos(0.0)*2 pi2=2*pi * LogZ=log(z-zeta) Alog=real(LogZ) ArgZ=imag(LogZ) * x=real(z) y=imag(z) xsi=real(zeta) eta=imag(zeta) xsiM=real(zetaM) etaM=imag(zetaM) * Case I if(xsi.ge.xsiM .and. eta.ge.etaM) then if(x.lt.xsiM .and. y.ge.eta) then ArgZ=ArgZ-pi2 elseif(x.ge.xsiM .and. y.lt.eta .and. * (y-etaM)*(xsi-xsiM).ge.(eta-etaM)*(x-xsiM)) then ArgZ=ArgZ+pi2 endif * Case II elseif(xsi.lt.xsiM .and. eta.ge.etaM) then if(x.lt.xsiM .and. * (y.ge.eta .or.(y-eta)*(xsiM-xsi).ge.(etaM-eta)*(x-xsi))) then ArgZ=ArgZ-pi2 endif * Case III elseif(xsi.lt.xsiM .and. eta.lt.etaM) then if(x.lt.xsiM .and. y.ge.eta .and. * (y-eta)*(xsiM-xsi).ge.(etaM-eta)*(x-xsi)) then ArgZ=ArgZ-pi2 endif * Case IV else if(y.ge.eta .and. * (x.lt.xsiM .or.(y-etaM)*(xsi-xsiM).lt.(eta-etaM)*(x-xsiM))) then ArgZ=ArgZ-pi2 endif endif LogNK=cmplx(ALog,ArgZ) end ************************************************************************ complex function WGintHF(z,zeta2,zeta1,K,zetaM,zetaE) * WGintHF * For a Flow aroud a Submerged Hydro-Foil * Integrated Wave Vortex Segment * with zeta1, zeta2 and zetaM located clockwisely * Cut of log(z-zeta) is set as * zeta1,2 -- zetaM -- zetaE(T.E.) -- (-inf,etaE) * Cut of log(z-zetaBar) as * zetaBar -- (-inf,etaBar) * Note : Use vlog if you want to change the Cut to * zetaBar -- (xsiBar,inf) * Cut of E1 as * zeteBar --(xsiBar,inf) * * 0 < K < 50 : Usual Wave Making Condition * K >= 50 : Limit Case of Low Speed ( Mirror Image Flow) * K = 0 : Limit Case of High Speed ( Inverse Mirror Image Flow) * K < 0 : Infinite Region Flow * real K complex i/(0,1)/ complex WGint1,WGint2 complex z,zeta1,zeta2,zetaM,zetaE complex dzeta,zetaBar,zdif,zdifBar complex LogZ,LogZBar,LogHF complex EiG,miKz,EmiKz,E1,SKappa pi=acos(0.0)*2 pi2=pi*2 * x=real(z) xsi1=real(zeta1) xsi2=real(zeta2) dzeta=zeta2-zeta1 * * For zeta2 * LogZ=LogHF(z,zeta2,zetaM,zetaE) zdif=z-zeta2 * zetaBar=conjg(zeta2) zdifBar=z-zetaBar LogZBar=log(zdifBar) EiG=dzeta/abs(dzeta) if(K.lt.0) then WGint2=i*conjg(EiG)*zdif*(LogZ-1) elseif(K.eq.0) then WGint2=i*conjg(EiG)*zdif*(LogZ-1)+i*EiG*zdifBar*(LogZBar-1) elseif(K.lt.50) then miKz=-i*K*zdifBar EmiKz=exp(miKz) call expint(miKZ,E1) if(x-xsi2.ge.0) then SKappa=EmiKz*E1 else SKappa=EmiKz*E1+pi2*i*EmiKz endif WGint2=i*conjg(EiG)*zdif*(LogZ-1)-i*EiG*zdifBar*(LogZBar-1) * +2/K*EiG*LogZBar+2/K*EiG*SKappa else WGint2=i*conjg(EiG)*zdif*(LogZ-1)-i*EiG*zdifBar*(LogZBar-1) endif * * For zeta1 * LogZ=LogHF(z,zeta1,zetaM,zetaE) zdif=z-zeta1 * zetaBar=conjg(zeta1) zdifBar=z-zetaBar LogZBar=log(zdifBar) ** EiG=dzeta/abs(dzeta) if(K.lt.0) then WGint1=i*conjg(EiG)*zdif*(LogZ-1) elseif(K.eq.0) then WGint1=i*conjg(EiG)*zdif*(LogZ-1)+i*EiG*zdifBar*(LogZBar-1) elseif(K.lt.50) then miKz=-i*K*zdifBar EmiKz=exp(miKz) call expint(miKZ,E1) if(x-xsi1.ge.0) then SKappa=EmiKz*E1 else SKappa=EmiKz*E1+pi2*i*EmiKz endif WGint1=i*conjg(EiG)*zdif*(LogZ-1)-i*EiG*zdifBar*(LogZBar-1) * +2/K*EiG*LogZBar+2/K*EiG*SKappa else WGint1=i*conjg(EiG)*zdif*(LogZ-1)-i*EiG*zdifBar*(LogZBar-1) endif WGintHF=WGint2-WGint1 end ************************************************************************ complex function VGintHF(z,zeta2,zeta1,K,zetaM,zetaE) * VGintHF * For a Flow aroud a Submerged Hydro-Foil * Integrated Wave Vortex Segment * with zeta1, zeta2 and zetaM located clockwisely * Cut of log(z-zeta) is set as * zeta1,2 -- zetaM -- zetaE(T.E.) -- (-inf,etaE) * Cut of log(z-zetaBar) as * zetaBar -- (-inf,etaBar) * Note : Use vlog if you want to change the Cut to * zetaBar -- (xsiBar,inf) * Cut of E1 as * zeatBar --(xsiBar,inf) * * 0 < K < 50 : Usual Wave Making Condition * K >= 50 : Limit Case of Low Speed ( Mirror Image Flow) * K = 0 : Limit Case of High Speed ( Inverse Mirror Image Flow) * K < 0 : Infinite Region Flow * real K complex i/(0,1)/ complex VGint1,VGint2 complex z,zeta1,zeta2,zetaM,zetaE complex dzeta,zetaBar,zdifBar complex LogZ,LogZBar,LogHF complex EiG,miKz,EmiKz,E1,SKappa pi=acos(0.0)*2 pi2=pi*2 * x=real(z) xsi1=real(zeta1) xsi2=real(zeta2) dzeta=zeta2-zeta1 * * For zeta2 * LogZ=LogHF(z,zeta2,zetaM,zetaE) * zetaBar=conjg(zeta2) zdifBar=z-zetaBar LogZBar=Log(zdifBar) EiG=dzeta/abs(dzeta) if(K.lt.0) then VGint2=i*conjg(EiG)*LogZ elseif(K.eq.0) then VGint2=i*conjg(EiG)*LogZ+i*EiG*LogZBar elseif(K.lt.50) then miKz=-i*K*zdifBar EmiKz=exp(miKz) call expint(miKZ,E1) if(x-xsi2.ge.0) then SKappa=EmiKz*E1 else SKappa=EmiKz*E1+pi2*i*EmiKz endif VGint2=i*conjg(EiG)*LogZ-i*EiG*LogZBar * -2*i*EiG*SKappa else VGint2=i*conjg(EiG)*LogZ-i*EiG*LogZBar endif * * For zeta1 * LogZ=LogHF(z,zeta1,zetaM,zetaE) * zetaBar=conjg(zeta1) zdifBar=z-zetaBar LogZBar=Log(zdifBar) ** EiG=dzeta/abs(dzeta) if(K.lt.0) then VGint1=i*conjg(EiG)*LogZ elseif(K.eq.0) then VGint1=i*conjg(EiG)*LogZ+i*EiG*LogZBar elseif(K.lt.50) then miKz=-i*K*zdifBar EmiKz=exp(miKz) call expint(miKZ,E1) if(x-xsi1.ge.0) then SKappa=EmiKz*E1 else SKappa=EmiKz*E1+pi2*i*EmiKz endif VGint1=i*conjg(EiG)*LogZ-i*EiG*LogZBar * -2*i*EiG*SKappa else VGint1=i*conjg(EiG)*LogZ-i*EiG*LogZBar endif VGintHF=VGint2-VGint1 end ************************************************************************ complex function LogHF(z,zeta,zetaM,zetaE) * * Complex Logarithmic Function for Hydro-Foil ( See Section 3.3.4 ) * * Cut of Log(z-zeta): * (-Inf,etaE) - zetaE - zetaM -zeta * * complex z,zeta,zetaE,zetaM,LogZ * pi=acos(0.0)*2 pi2=2*pi * LogZ=log(z-zeta) Alog=real(LogZ) ArgZ=imag(LogZ) * y=imag(z) eta=imag(zeta) etaE=imag(zetaE) etaM=imag(zetaM) * AlpE=imag(log(zetaE-zeta)) AlpM=imag(log(zetaM-zeta)) ArgEM=imag(log(zetaE-zetaM)) ArgME=imag(log(zetaM-zetaE)) ArgZE=imag(log(z-zetaE)) ArgZM=imag(log(z-zetaM)) if(AlpE.lt.0) then * Case I if(AlpM.ge.AlpE+pi) then * Case I-1 if(y.lt.eta) then if(y.ge.etaE .and. ArgZM.lt.ArgEM) then ArgZ=ArgZ+pi2 endif else if(ArgZM.ge.ArgEM .and. ArgZM.lt.AlpM-pi) then ArgZ=ArgZ-pi2 endif endif elseif(AlpM.ge.0) then * Case I-2 if(y.ge.etaE) then if(y.lt.eta.and.ArgZ.lt.AlpE .or. * ArgZ.ge.AlpE.and.ArgZ.lt.AlpM.and.ArgZM.lt.ArgEM) then ArgZ=ArgZ+pi2 endif endif elseif(AlpM.ge.AlpE) then * Case I-3 if(y.ge.etaE.and.y.lt.eta.and.ArgZ.lt.AlpE .or. * ArgZ.lt.AlpM.and.ArgZ.ge.AlpE.and.ArgZE.ge.ArgME) then ArgZ=ArgZ+pi2 endif elseif(ArgME.ge.AlpE+pi) then * Case I-4 if(y.lt.eta.and.y.ge.etaM.and.ArgZ.lt.AlpM .or. * y.lt.etaM.and.y.ge.etaE.and.ArgZE.ge.ArgME) then ArgZ=ArgZ+pi2 endif else * Case I-5 Out of Rule if(y.ge.etaE.and.y.lt.eta.and.ArgZ.lt.AlpE) then ArgZ=ArgZ+pi2 endif endif else * Case II if(AlpM.lt.AlpE-pi) then * Case II-1 if(y.ge.eta) then if(y.lt.etaE .and. ArgZM.ge.ArgEM) then ArgZ=ArgZ-pi2 endif else if(ArgZM.lt.ArgEM .and. ArgZM.ge.AlpM+pi) then ArgZ=ArgZ+pi2 endif endif elseif(AlpM.lt.0) then * Case II-2 if(y.lt.etaE) then ArgZM=imag(log(z-zetaM)) if(y.ge.eta.and.ArgZ.ge.AlpE .or. * ArgZ.le.AlpE.and.ArgZ.ge.AlpM.and.ArgZM.ge.ArgEM) then ArgZ=ArgZ-pi2 endif endif elseif(AlpM.lt.AlpE) then * Case II-3 if(y.lt.etaE.and.y.ge.eta.and.ArgZ.ge.AlpE .or. * ArgZ.ge.AlpM.and.ArgZ.lt.AlpE.and.ArgZE.lt.ArgME) then ArgZ=ArgZ-pi2 endif elseif(ArgME.lt.AlpE-pi) then * Case II-4 if(y.ge.eta.and.y.lt.etaM.and.ArgZ.ge.AlpM .or. * y.ge.etaM.and.y.lt.etaE.and.ArgZE.lt.ArgME) then ArgZ=ArgZ-pi2 endif else * Case II-5 Out of Rule if(y.lt.etaE.and.y.ge.eta.and.ArgZ.ge.AlpE) then ArgZ=ArgZ-pi2 endif endif endif * LogHF=cmplx(ALog,ArgZ) end ************************************************************************ complex function WGint(z,zeta2,zeta1,K,zetaM) * WGintNK * real K complex z,zeta1,zeta2,zetaM complex WGintNK * WGint=WGintNK(z,zeta2,zeta1,K,zetaM) end ************************************************************************ complex function VGint(z,zeta2,zeta1,K,zetaM) * VGintNK * real K complex z,zeta1,zeta2,zetaM complex VGintNK * VGint=VGintNK(z,zeta2,zeta1,K,zetaM) end ************************************************************************ SUBROUTINE EXPINT(z,e1) * * EXPONENTIAL INTEGRAL OF COMPLEX ARGUMENT * CODED BY KATSUO SUZUKI (BODAI) * IMPLICIT COMPLEX (Z) REAL INF/9.99999999999999D09/ COMPLEX E1,E1GL,EXPE1,LOGZ REAL GAMMA/0.5772156649D00/ REAL PAI/3.14159265359D0/ ** RETURN * ** ENTRY EXPI1 (Z,E1) ** ENTRY EXPI2 (Z,E1GL) ** ENTRY EXPI3 (Z,EXPE1) ** ENTRY EXPI4 (Z,ZEXPE1) ** ENTRY EXPIAL (Z,E1,E1GL,EXPE1,ZEXPE1) * * E1 =COMPLEX EXPONENTIAL INTEGRAL OF COMPLEX ARGUMENT Z * E1GL =E1 + GAMMA + LOG( Z ) * EXPE1 = EXP( Z ) * E1 * ZEXPE1 = Z * EXPE1 * R=ABS(Z) X=real(Z) Y=IMAG(Z) IF(R.EQ.0.0) GO TO 10 LOGZ=LOG(Z) IF(R.GT.25.0) GO TO 40 IF(R.GT.4.0.AND.X.GE.0.0) GO TO 30 IF(ABS(Y).GT.4.0) GO TO 30 GO TO 20 * * R = 0 * 10 CONTINUE E1=DCMPLX(INF,0.0) E1GL=(1.0,0.0) EXPE1=E1 ZEXPE1=DCMPLX(GAMMA,0.0) WRITE(6,100) 100 FORMAT(' *** (EXPINT) Z=(0.0,0.0) ***') RETURN * * SERIES EXPANSION * * E1GL= - S (-Z)**N/(N*FACT(N)) * N=1 * 20 CONTINUE I=1 ZI=Z ZSUM=ZI 21 I=I+1 A=-DFLOAT(I-1)/DFLOAT(I*I) ZI=A*ZI*Z ZSUM=ZSUM+ZI RE=ABS(ZSUM) AD=ABS(ZI)*1.0E8 IF(RE.LE.AD) GO TO 21 E1GL=ZSUM E1=E1GL-GAMMA-LOGZ EXPE1=EXP(Z)*E1 ZEXPE1=Z*EXPE1 RETURN * * CONTINUED FRACTION * * EXPE1 = 1/Z+1/1+1/Z+2/1+2/Z+3/1+3/Z+ * 30 CONTINUE I=11 Z0=(0.0,0.0) ZI=Z0 31 IF(I.EQ.1) GO TO 35 I=I-1 AI=DFLOAT(I) ZD=Z+ZI IF(ZD.NE.0.0) GO TO 32 ZI=Z0 GO TO 33 32 ZI=AI/ZD 33 CONTINUE ZD=1.0+ZI IF(ZD.NE.0.0) GO TO 34 ZI=Z0 GO TO 31 34 CONTINUE ZI=AI/ZD GO TO 31 35 CONTINUE EXPE1=1.0/(Z+ZI) E1=EXP(-Z)*EXPE1 E1GL=E1+GAMMA+LOGZ ZEXPE1=Z*EXPE1 RETURN * * ASYMPTOTIC EXPANSION * * ZEXPE1 = 1-1/Z+FACT(2)/Z**2-FACT(3)/Z**3+ * 40 CONTINUE ZINV=1.0/Z I=0 ZI=(1.0,0.0) ZSUM=ZI 41 I=I+1 AI=DFLOAT(I) ZI=-AI*ZINV*ZI ZSUM=ZSUM+ZI RE=ABS(ZSUM) AD=ABS(ZI)*1.0E8 IF(RE.LE.AD) GO TO 41 IF(Y.EQ.0.0.AND.X.LT.0.0) ZSUM=ZSUM-Z*EXP(Z)*DCMPLX(0.0,PAI) ZEXPE1=ZSUM EXPE1=ZEXPE1*ZINV IF(-X.GE.200.0) GO TO 42 E1=EXP(-Z)*EXPE1 E1GL=E1+GAMMA+LOGZ RETURN 42 E1=INF*EXP(DCMPLX(0.0,-Y))*EXPE1 E1GL=E1+GAMMA+LOGZ RETURN END ************************************************************************ subroutine StrmLnPlt(NxD, Nx,Ny,xpos,ypos, Phi,Psi, L0,z00,ze0,dp) * real xpos(NxD,*),ypos(NxD,*) real Phi(NxD,*),Psi(NxD,*) * Nxp1=Nx+1 Nyp1=Ny+1 dp=(ze0-z00)/L0 * call lcolor(3,4) call minmaxV(NxD,Nxp1,Nyp1,Phi,PhiMax,PhiMin,0) write(6,'(a,2f10.3)') '*** Min,Max of Phi(i,j) =',PhiMax,PhiMin call minmaxV(NxD,Nxp1,Nyp1,Psi,PsiMax,PsiMin,0) write(6,'(a,2f10.3)') '*** Min,Max of Psi(i,j) =',PsiMax,PsiMin * * Phi > 0 z0=0 ze=ze0 L=L0*(ze-z0)/(ze0-z00) call color(2) call TOKOSN(NxD,Nxp1,Nyp1,xpos,ypos,Phi, z0,ze,L) call display * Phi =< 0 z0=z00 ze=0 L=L0*(ze-z0)/(ze0-z00) call color(6) call TOKOSN(NxD,Nxp1,Nyp1,xpos,ypos,Phi, z0,ze,L) call display * Psi > 0 z0=0 ze=ze0 L=l0*(ze-z0)/(ze0-z00) call color(5) call TOKOSN(NxD,Nxp1,Nyp1,xpos,ypos,Psi, z0,ze,L) call display * Psi =< 0 z0=z00 ze=0 L=l0*(ze-z0)/(ze0-z00) call color(4) call TOKOSN(NxD,Nxp1,Nyp1,xpos,ypos,Psi, z0,ze,L) call display end ************************************************************************ SUBROUTINE TOKOSN(nd,n,m,xpos,ypos,zval, ZI,ZE,L) * real xpos(nd,*),ypos(nd,*),zval(nd,*) parameter ( lmax=400, lmaxp1=lmax+1 ) DIMENSION X(3),Y(3),Z(3) COMMON /T1/ZEQ(lmaxp1),zmax,zmin,LL * IF(L.GT.lmax.OR.L.LT.0) then WRITE(6, * '('' *** (TOKOSN) 9TH ARG. .GT.'',i4,'' .OR. .LT.0 ***'')') * lmax STOP else IF(N.LE.0.OR.M.LE.0) then WRITE(16, * '('' *** (TOKOSN) ARGUMENT N OR M .LE. 0 ***'')') STOP else LL=L llp1=ll+1 ZMin=MIN(ZI,ZE) zmax=max(zi,ze) if(ll.eq.0) then dz=0 else dz=(zmax-zmin)/ll endif DO I=1,Llp1 ZEQ(I)=ZMin+DZ*(I-1) enddo DO J=1,M-1 x(2)=xpos(1,j) y(2)=ypos(1,j) z(2)=zval(1,j) x(3)=xpos(1,j+1) y(3)=ypos(1,j+1) z(3)=zval(1,j+1) DO I=1,N-1 ** write(6,*) ' i,j =',i,j X(1)=X(3) Y(1)=Y(3) Z(1)=Z(3) X(3)=X(2) Y(3)=Y(2) Z(3)=Z(2) X(2)=xpos(i+1,j) Y(2)=ypos(i+1,j) Z(2)=zval(i+1,j) CALL TPLOT(X,Y,Z) ** call display ** pause X(3)=xpos(i+1,j+1) Y(3)=ypos(i+1,j+1) Z(3)=zval(i+1,j+1) CALL TPLOT(X,Y,Z) ** call display ** pause enddo enddo endif endif call penup RETURN END ************************************************************************ SUBROUTINE TPLOT(X,Y,T) * * ***SUBPROGRAM-1 OF TOKOSN *** * parameter ( lmax=400, lmaxp1=lmax+1 ) COMMON /T1/ Z(lmaxp1),zmax,zmin,L DIMENSION X(3),Y(3),T(3) dimension XP(lmaxp1),YP(lmaxp1),XQ(lmaxp1),YQ(lmaxp1) DATA XE,YE/1.0E30,1.0E30/ * ZMA=MAX(T(1),T(2),T(3)) ZMI=MIN(T(1),T(2),T(3)) IF(ZMA.ge.Zmin.and.ZMI.le.Zmax) then if(zma.eq.t(1)) then nma=1 if(zmi.eq.t(3)) then nmi=3 nmd=2 else nmi=2 nmd=3 endif elseif(zma.eq.t(2)) then nma=2 if(zmi.eq.t(1)) then nmi=1 nmd=3 else nmi=3 nmd=1 endif else nma=3 if(zmi.eq.t(1)) then nmi=1 nmd=2 else nmi=2 nmd=1 endif endif ZMD=T(NMD) CALL TNODE(ZMA,ZMI,ZMD,IMA,IMI,IMAM,IMIM,IND) ** write(6,*) 'ind,ima-imi =',ind,ima,imim,imam,imi IF(IND.ne.0) then XA=X(NMA) YA=Y(NMA) ZA=T(NMA) XB=X(NMI) YB=Y(NMI) ZB=T(NMI) XC=X(NMD) YC=Y(NMD) ZC=T(NMD) CX=(XA-XB)/(ZA-ZB) CY=(YA-YB)/(ZA-ZB) DO I=IMI,IMA XP(I)=XB+CX*(Z(I)-ZB) YP(I)=YB+CY*(Z(I)-ZB) ** write(6,*) 'P,i,x,y =',i,xp(i),yp(i) enddo if(imi.le.imam.and.zc.ne.zb) then CX=(XC-XB)/(ZC-ZB) CY=(YC-YB)/(ZC-ZB) DO I=IMI,IMAM XQ(I)=XB+CX*(Z(I)-ZB) YQ(I)=YB+CY*(Z(I)-ZB) ** write(6,*) 'Q,i,x,y =',i,xq(i),yq(i) enddo endif if(imim.le.ima.and.za.ne.zc) then CX=(XA-XC)/(ZA-ZC) CY=(YA-YC)/(ZA-ZC) DO I=IMIM,IMA XQ(I)=XC+CX*(Z(I)-ZC) YQ(I)=YC+CY*(Z(I)-ZC) ** write(6,*) 'Q,i,x,y =',i,xq(i),yq(i) enddo endif DO I=IMI,IMA if(mod(i-1,5).eq.0) then call newpen(3) else call newpen(1) endif IF(I.EQ.IMI.AND.XP(I).EQ.XE.AND.YP(I).EQ.YE) then CALL DRAW(XQ(I),YQ(I)) XE=XQ(I) YE=YQ(I) else IF(I.ne.IMI.or.XQ(I).ne.XE.or.YQ(I).ne.YE) then CALL PENUP IF(XP(I).GT.XQ(I)) then CALL DRAW(XQ(I),YQ(I)) CALL DRAW(XP(I),YP(I)) XE=XP(I) YE=YP(I) else CALL DRAW(XP(I),YP(I)) CALL DRAW(XQ(I),YQ(I)) XE=XQ(I) YE=YQ(I) endif else CALL DRAW(XP(I),YP(I)) XE=XP(I) YE=YP(I) endif endif enddo endif endif END ************************************************************************ SUBROUTINE TNODE(ZMA,ZMI,ZMD,IMA,IMI,IMAM,IMIM,IND) * * ***SUBPROGRAM-2 OF TOKOSN *** * parameter ( lmaxp1=400+1 ) COMMON /T1/Z(lmaxp1),zmax,zmin,L * ** write(6,'(a,3f10.3)') 'ZMA,ZMI,ZMD =',ZMA,ZMI,ZMD ** write(6,'(a/3(5(i3,f10.3)/))') 'i,z(i) =',(i,z(i),i=1,L+1) if(zmi.gt.zmax) then ind=0 else Lp1=L+1 do i=1,Lp1 ** write(6,'(a,i5,2f10.3)') 'NODE-1 i,z(i) =',i,z(i),zmi if(z(i).ge.zmi) goto 10 enddo 10 continue imi=i ** write(6,'(a,i5)') '***** imi =',imi if(zma.lt.zmin) then ind=0 else do i=Lp1,imi,-1 ** write(6,'(a,i5,2f10.3)') 'NODE-2 i,z(i) =',i,z(i),zma if(z(i).lt.zma) goto 20 enddo 20 continue ima=i ** write(6,'(a,i5)') '***** ima =',ima do i=imi,ima ** write(6,'(a,i5,2f10.3)') 'NODE-3 i,z(i) =',i,z(i),zmd if(z(i).ge.zmd) goto 30 enddo 30 continue imim=i ** write(6,'(a,i5)') '***** imim =',imim do i=imim,imi,-1 ** write(6,'(a,i5,2f10.3)') 'NODE-4 i,z(i) =',i,z(i),zmd if(z(i).lt.zmd) goto 40 enddo 40 continue imam=i ** write(6,'(a,i5)') '***** imam =',imam ind=1 endif endif END