program buckling_stability_diagram *--------------------------------------------------------------------------- * * This Fortran program calculates the stability diagram of irrotational * and 3D-incompressible buckling modes (lambda>0) * * Delta a_0=(q/r)*a_{0,r} * * with free surface boundary conditions * * -i sigma h + uh0_x=0 (inner,outer) * -i sigma h + wh0_z=0 (upper,lower) * * The equation for a_0 is a leading order solution to a cylindrical * expansion * * Phi = Sigma a_n(r,t,theta) z^n, * * that decouples horizontal from vertical motion. * * The stabilitiy diagram is calculated for arbitrarily wide tori for the * purpose of studying its gravitational wave emissions. Gravitational * wave emission will be long lasting when energized by feedback from * a central rotating black hole. * * We extend to wide tori [1] the wave instabilities for infinitesimally * slender tori of Papaloizou-Pringle [2,3]. Our new results demonstrate * that the lowest-order non-axisymmetric modes become unstable first. * Quadrupole gravitational wave emission is hereby expected to be the * dominant emission channel [4,5]. Around stellar mass black holes, the * emission frequencies fall into the bandwidth of sensitivity of LIGO-Virgo * and KAGRA. * * References: * 1. van Putten, M.H.P.M., 2002, ApJ Lett., 575, L71 * 2. Papaloizou, J. C. B., & Pringle, J. E. 1984, MNRAS, 208, 721 * 3. Goldreich, P., Goodman, J., & Narayan, R. 1986, MNRAS, 221, 339 * 4. van Putten, M.H.P.M., & Levinson, A., 2003, ApJ, 584, 937 * 5. van Putten, M.H.P.M., 2012, Prog. Theor. Phys., 127, 331 * See further * van Putten, M.H.P.M., 2005, Gravitational Radiation, Luminous Black * Holes and Gamma-Ray Burst Supernovae (Cambridge University Press) * van Putten, M.H.P.M., and Levinson, A.,2012, Relativistic Astrophysics * of the Transient Universe (Cambridge University Press) * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * * Copyright 2002,2004,2012,2016 Maurice H.P.M. van Putten * This program is distributed under the GNU General Public Licence * * You should have received a copy of the GNU General Public License * along with this program. If not, see . * *---------------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) real*8 FU(2),dx(2),DF(2,2) real*8 n_n,n_p,NP real derr,err,max_err,err_b parameter (Nb=999) q=2. b=0.3 do 20,m=1,3 do 20,i=1,Nb omega=-2.+4.*real(i)/Nb write(40+m,*)real(omega),real(G(omega,q,m,b)) 20 close(40+m) IUNIT=300 JUNIT=400 i_select=1 if(i_select.eq.2)then omega0=-0.05 q0= 2.5 else q0=1.733 omega0=0. endif do 2,m=1,7 if(m.eq.3.and.i_select.eq.2)then b=0.2 omega0=-0.253033111299451 q0=2.149100 fla=flambda(omega0,m,q0,b) fla=2.5759772071693 fs=CHECK(omega0,fla,q0,m,b) write(*,*)'---- m=3 initial ----',real(fs),real(fla) write(*,*)'etc:',real(q0),real(omega0),real(fla) endif omega=omega0 q=q0 err=1.E-00 max_err=0.E-00 max_ite=0 iwords=0 do 1,i_b=1,Nb q_old=q if(i_select.eq.2)then b=0.2+0.0008*real(i_b) else b=0.001*real(i_b) endif err=1.E-00 iter=0 10 iter=iter+1 call Equations(FU,m,q,b,omega) call Jacobian(DF,m,q,b,omega) call solve(DF,dx,FU) omega=omega-dx(1) q=q-dx(2) err=max(abs(FU(1)),abs(FU(2))) if(err.gt.1.E-06.and.iter.lt.40)goto 10 call Equations(FU,m,q,b,omega) om=omega write(90+m,*)real(b),real(FU(1)),real(FU(2)),real(G(om,q,m,b)) if(q.gt.1.7.and.err.lt.1.D-06)then call zerosH0(q,x,b,it0,f_e) fla=flambda(omega,m,q,b) fs=CHECK(omega,fla,q,m,b) write(iUNIT+m,*)real(b),real(q),real(omega),real(fla),real(fs) s1=omega-m*dom(q,x) s2=omega-m*dom(q,b) call turning_point(omega,q,y,m,i_f,f_e) write(jUNIT+m,*)real(b),real(x),real(y),real(s1),real(s2) iwords=iwords+1 derr=err its=iter if(q_old.lt.2.D-00 .and. q.gt.2.D-00)then call turning_point(omega,q,x,m,i_f,f_e) fla=flambda(omega,m,q,b) write(200,*)real(fla),real(omega),real(x/b),i_f,real(f_e) call zerosH0(q,x,b,i_f,f_e) write(100,*)real(b),real(fla),real(omega),real(x/b),real(m*b) rs1=phi(x,fla,m,q)*m/(1+x)+dom(q,x)*(1+x) rs2=phi(b,fla,m,q)*m/(1+b)+dom(q,b)*(1+b) rs3=phi_x(x,fla,m,q) rs4=phi_x(b,fla,m,q) write(101,*)real(fla),real(rs1),real(rs2),real(rs3),real(rs4) p=sqrt(q*q+4*m*m) n_p=(q+p)/(2*m) n_n=(q-p)/(2*m) siP=omega-m*dom(q,b) S2 =siP**2 NP= -m*H1(q,b)-(2-q)*siP write(105,*)real(fla),real(siP),real(S2),real(NP),real(S2/NP) else q_old=q endif if(i_b.eq.1)then call turning_point(omega,q,x,m,i_f,f_e) fla=flambda(omega,m,q,b) write(102,*)real(fla),real(omega),real(x/b),i_f,real(f_e) p=sqrt(q*q+4*m*m) n_p=(q+p)/(2*m) n_n=(q-p)/(2*m) siP=omega-m*dom(q,b) S2 =siP**2 NP= -m*H1(q,b)-(2-q)*siP write(104,*)real(fla),real(siP),real(S2),real(NP),real(S2/NP) call zerosH0(q,x,b,i_f,f_e) rs1=phi(x,fla,m,q)*m/(1+x) rs2=phi(b,fla,m,q)*m/(1+b) rs3=phi_x(x,fla,m,q) rs4=phi_x(b,fla,m,q) write(103,*)real(fla),real(rs1),real(rs2),real(rs3),real(rs4) endif endif if(i_b.eq.1)then omega0=omega q0=q endif max_err=max(max_err,derr) 1 max_ite=max(max_ite,its) write(*,*)'lower results for m=',m write(*,*)' q=',q write(*,*)' omega=',omega write(*,*)'err,iter=',max_err,max_ite 2 write(*,*)' iwords=',iwords end subroutine solve(DF,dx,FU) *----------------------------------------------------------------------------- * *----------------------------------------------------------------------------- real*8 DF(2,2),dx(2),FU(2),det det=DF(1,1)*DF(2,2)-DF(1,2)*DF(2,1) dx(1)= DF(2,2)*FU(1)-DF(1,2)*FU(2) dx(2)=-DF(2,1)*FU(1)+DF(1,1)*FU(2) dx=dx/det return end subroutine Jacobian(DF,m,q,b,omega) *----------------------------------------------------------------------------- * *----------------------------------------------------------------------------- integer m real*8 DF(2,2),FU(2),fP(2),fM(2),omega real*8 q,b,eps parameter (eps=1.D-07) call Equations(fP,m,q,b,omega+eps) call Equations(fM,m,q,b,omega-eps) DF(:,1)=(fP-fM)/(2*eps) call Equations(fP,m,q+eps,b,omega) call Equations(fM,m,q-eps,b,omega) DF(:,2)=(fP-fM)/(2*eps) return end subroutine Equations(FU,m,q,b,omega) *----------------------------------------------------------------------------- * *----------------------------------------------------------------------------- integer m real*8 q,b,omega,FU(2),G,DG FU(1)= G(omega,q,m,b) FU(2)=DG(omega,q,m,b) return end real*8 function G(om,q,m,b) *----------------------------------------------------------------------------- * *----------------------------------------------------------------------------- implicit real*8 (a-z) integer m,iter p=sqrt(q*q+4*m*m) n_p=(q+p)/(2*m) n_n=(q-p)/(2*m) call zerosH0(q,x,b,iter,err) siP=om-m*dom(q,b) siM=om-m*dom(q,x) SP=siP**2 SM=siM**2 kP=m/(1+b) kM=m/(1+x) NP= -kP*H1(q,b)-(2-q)*siP/(1+b)**q NM= -kM*H1(q,x)-(2-q)*siM/(1+x)**q alpha=dabs(1+x)**p-(1+b)**p beta =dabs(1+x)**p+(1+b)**p beta1=n_n*(1+b)**p-n_p*dabs(1+x)**p beta2=n_p*(1+b)**p-n_n*dabs(1+x)**p G=alpha*(SP*SM+NP*NM*n_p*n_n)+beta1*NM*SP+beta2*NP*SM return end real*8 function DG(om,q,m,b) *----------------------------------------------------------------------------- * *----------------------------------------------------------------------------- implicit real*8 (a-z) integer m parameter (eps=1.D-05) DG=(G(om+eps,q,m,b)-G(om-eps,q,m,b))/(2.D-00*eps) return end real*8 function dom(q,x) *----------------------------------------------------------------------------- * *----------------------------------------------------------------------------- implicit real*8 (a-z) dom=1/(1+x)**q-1 return end subroutine zerosH0(q,x,y,iter,err) *----------------------------------------------------------------------------- * *----------------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) Hy =H(q,y) err =1.D-00 iter=0 xL =-0.9 xR = 0. 1 iter=iter+1 x=(xL+xR)/2 if(sign(1.D-00,H(q,xL)-Hy)*sign(1.D-00,H(q,x)-Hy).lt.0.)then xR=x else xL=x endif err=abs(H(q,x)-Hy) if(err.gt.2.D-14.and.iter.lt.100)goto 1 if(iter.eq.100)write(99,*)real(q),real(x),real(y),real(err),iter return end real*8 function H(q,x) *----------------------------------------------------------------------------- * *----------------------------------------------------------------------------- implicit real*8 (a-z) e=2-2*q H=1/(1+x)+(1+x)**e/e return end real*8 function H1(q,x) implicit real*8 (a-z) e=1-2*q H1=(1+x)**e-1/(1+x)**2 return end subroutine turning_point(om,q,x,m,iter,err) *------------------------------------------------------------------------------- * Buckling modes *------------------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) err =1.D-00 iter=0 xL =-0. xR = 4. 1 iter=iter+1 x=(xL+xR)/2 if(sign(1.D-00,m*dom(q,xL)-om)*sign(1.D-00,m*dom(q,x)-om).lt.0.) & then xR=x else xL=x endif err=abs(m*dom(q,x)-om) if(err.gt.2.D-10.and.iter.lt.100)goto 1 return end real*8 function flambda(om,m,q,b) *------------------------------------------------------------------------------- * *------------------------------------------------------------------------------- implicit real*8 (a-z) integer m p=sqrt(q*q+4*m*m) nh_p=(q+p)/(2*m) nh_n=(q-p)/(2*m) siP=om-m*dom(q,b) S2 =siP**2 kP=m/(1+b) NP= -kP*H1(q,b)-(2-q)*siP/(1+b)**q flambda=(NP*nh_p/S2-1)/(1-NP*nh_n/S2)*(1+b)**p return end real*8 function phi(x,la,m,q) *------------------------------------------------------------------------------- * *------------------------------------------------------------------------------- implicit real*8 (a-z) integer m p =sqrt(q*q+4*m*m) n_p=(q+p)/2.D-00 n_n=(q-p)/2.D-00 phi=dabs(1+x)**n_p+la*dabs(1+x)**n_n return end real*8 function phi_x(x,la,m,q) *------------------------------------------------------------------------------- * *------------------------------------------------------------------------------- implicit real*8 (a-z) integer m parameter (eps=1.D-06) p=sqrt(q*q+4*m*m) n_p=(q+p)/2 n_n=(q-p)/2 phi_x=n_p*(1+x)**(n_p-1)+la*n_n*(1+x)**(n_n-1) return end real*8 function CHECK(om,la,q,m,b) *------------------------------------------------------------------------------- * *------------------------------------------------------------------------------- implicit real*8 (a-z) integer m,iter call zerosH0(q,x,b,iter,err) G1=Gx(om,la,q,m,x) G2=Gx(om,la,q,m,b) G1=Gy(om,la,q,m,x) G2=Gy(om,la,q,m,b) CHECK=max(abs(G1),abs(G2)) return end real*8 function Gy(om,la,q,m,x) *------------------------------------------------------------------------------- * *------------------------------------------------------------------------------- implicit real*8 (a-z) integer m epsilon=(m/(1+x))*phi(x,la,m,q)/phi_x(x,la,m,q) siX=om-m*dom(q,x) SX =siX**2 NX =-(m/(1+x))*H1(q,x)-(2-q)*siX/(1+x)**q Gy=epsilon*SX-Nx return end real*8 function Gx(om,la,q,m,x) *---------------------------------------------------------- * Separate out lambda coefficient *---------------------------------------------------------- implicit real*8 (a-z) integer m siX=om-m*dom(q,x) SX=siX**2 kX=m/(1+x) NX=-kX*H1(q,x)-(2-q)*siX/(1+x)**q p=sqrt(q*q+4*m*m) nh_p=(q+p)/(2*m) nh_n=(q-p)/(2*m) Gx = la-(nh_p*NX/SX-1)/(1-nh_n*NX/SX)*(1+x)**p return end real*8 function G_B(om,q,m,b) *------------------------------------------------------------------------------- * *------------------------------------------------------------------------------- implicit real*8 (a-z) integer m,iter la=0.D-00 call zerosH0(q,x,b,iter,err) p=sqrt(q*q+4*m*m) nh_p=(q+p)/(2*m) nh_n=(q-p)/(2*m) siP=om-m*dom(q,b) siM=om-m*dom(q,x) SP=siP**2 SM=siM**2 kP=m/(1+b) kM=m/(1+x) NP= -kP*H1(q,b)-(2-q)*siP/(1+b)**q NM= -kM*H1(q,x)-(2-q)*siM/(1+x)**q fs=SP*SM*(1-nh_n*NP/SP)*(1-nh_n*NM/SM) G_B=(Gx(om,la,q,m,b)-Gx(om,la,q,m,x))*fs return end