"""Scalar Schwarzschild circular modes from official pybhpt homogeneous solutions. G=c=M=q_s=1, Box Phi=-4 pi rho, Phi=sum R_lm(r)Y_lm exp(-i m Omega t). The delta source gives [Delta R_lm']=-4 pi conjugate(Y_lm(pi/2,0))/u^t. Fr=q partial_r Phi is covariant Schwarzschild radial force. Both one-sided limits and the individual complex m>=0 contributions are retained. """ import argparse import importlib.metadata import json import math import platform import time from pathlib import Path import numpy as np from scipy.special import sph_harm_y from pybhpt.radial import RadialTeukolsky def encoded(z): return [float(z.real), float(z.imag)] def flux_check(radius,multipoles,lmax): """Independent inhomogeneous-source amplitudes and scalar stress fluxes.""" from pybhpt.geo import KerrGeodesic from pybhpt.teuk import TeukolskyMode from pybhpt.flux import FluxMode geo=KerrGeodesic(0.0,radius,0.0,1.0,nsamples=256) modes=[] for ell in range(1,lmax+1): for emm in range(1 if ell%2 else 2,ell+1,2): mode=TeukolskyMode(0,ell,emm,0,0,geo) mode.solve(geo,method='AUTO') flux=FluxMode(geo,mode).energy modes.append(dict(ell=ell,m=emm,weight=2, infinity=2*flux['I'],horizon=2*flux['H'])) infinity=math.fsum(m['infinity'] for m in modes) horizon=math.fsum(m['horizon'] for m in modes) local=math.sqrt(1-3/radius)*math.fsum( m['force_time'] for m in multipoles[:lmax+1]) difference=local-infinity-horizon if abs(difference)>1e-10*(infinity+horizon): raise ArithmeticError('Independent scalar energy-flux balance failed') return dict(lmax=lmax,infinity=infinity,horizon=horizon, local_Ft_over_ut=local,difference=difference, relative_difference=difference/(infinity+horizon),modes=modes) def main(): p=argparse.ArgumentParser(description=__doc__) p.add_argument('--radius',type=float,required=True) p.add_argument('--lmax',type=int,default=80) p.add_argument('--method',default='AUTO') p.add_argument('--rtol',type=float,default=None) p.add_argument('--flux-lmax',type=int,default=0) p.add_argument('--output',type=Path,required=True) args=p.parse_args() r=args.radius if r not in (6.0,10.0) or args.lmax < 0 or args.lmax>100: raise ValueError('This bounded benchmark requires r=6 or10 and 0<=lmax<=100') if not 0<=args.flux_lmax<=args.lmax: raise ValueError('Flux cutoff must lie between zero and the force cutoff') if args.output.exists(): raise FileExistsError(args.output) ut=1/math.sqrt(1-3/r) delta=r*(r-2) start=time.monotonic() modes=[] multipoles=[] for ell in range(args.lmax+1): fp=[]; fm=[]; ft=[] for emm in range(ell % 2,ell+1,2): check_radii=np.array([r-.05,r,r+.05]) radial=RadialTeukolsky(0,ell,emm,0.0,emm/r**1.5,check_radii) radial.solve(method=args.method,rtol=args.rtol) ri_all=radial.radialsolutions('In') di_all=radial.radialderivatives('In') ro_all=radial.radialsolutions('Up') do_all=radial.radialderivatives('Up') ri,di,ro,do=(v[1] for v in (ri_all,di_all,ro_all,do_all)) wr=delta*(ri*do-ro*di) wr_checks=check_radii*(check_radii-2)*(ri_all*do_all-ro_all*di_all) wr_drift=max(abs(wr_checks/wr-1)) if wr_drift>1e-10: raise ArithmeticError((ell,emm,'radial Wronskian is not constant')) # At the equator this factorial expression is exact before the # final float conversion. It avoids the angular recursion's # growing relative error in a cancellation-sensitive mode sum. angular_weight=((2*ell+1)*math.comb(ell+emm,(ell+emm)//2) *math.comb(ell-emm,(ell-emm)//2))/2**(2*ell) yy=angular_weight/(4*math.pi) angular_check=abs(sph_harm_y(ell,emm,math.pi/2,0.0))**2 if abs(yy-angular_check)>1e-12*yy: raise ArithmeticError('Equatorial harmonic normalization mismatch') pref=-angular_weight/ut/wr plus=pref*ri*do minus=pref*ro*di field=pref*ri*ro time_force=-1j*(emm/r**1.5)*field jump=-angular_weight/(ut*delta) if not all(np.isfinite(v) for v in (plus,minus,wr)): raise ArithmeticError((ell,emm,'nonfinite radial output')) if abs(plus-minus-jump)>1e-12*max(abs(jump),1e-6): raise ArithmeticError((ell,emm,'delta-source jump mismatch')) weight=1 if emm==0 else 2 fp.append(weight*plus.real) fm.append(weight*minus.real) ft.append(weight*time_force.real) modes.append(dict(ell=ell,m=emm,weight=weight, radial_in=encoded(ri),radial_in_derivative=encoded(di), radial_up=encoded(ro),radial_up_derivative=encoded(do), wronskian=encoded(wr),angular_norm=yy, wronskian_relative_drift=float(wr_drift), force_plus=encoded(plus),force_minus=encoded(minus), field=encoded(field),force_time=encoded(time_force), source_jump=jump, jump_residual=abs(plus-minus-jump))) plus=math.fsum(fp) minus=math.fsum(fm) expected=-(2*ell+1)/(ut*delta) if ell==0 and max(abs(plus+1/(ut*delta)),abs(minus))>1e-14: raise ArithmeticError('Static monopole analytic check failed') multipoles.append(dict(ell=ell,force_plus=plus,force_minus=minus, force_average=(plus+minus)/2, force_time=math.fsum(ft), jump_expected=expected, jump_residual=plus-minus-expected)) if ell%10==0: print(f'r={r:g} l={ell:2d} F+={plus:.15g} F-={minus:.15g} ' f'jump_error={plus-minus-expected:.3g} ' f'seconds={time.monotonic()-start:.2f}',flush=True) result=dict(radius=r,lmax=args.lmax,method=args.method,rtol=args.rtol, conventions='G=c=M=q_s=1; covariant Fr; scalar Box Phi=-4pi rho; ' 'unit-normalized spherical harmonics; m>0 pairs weighted2', python=platform.python_version(), dependencies={x:importlib.metadata.version(x) for x in ('pybhpt','numpy','scipy')}, elapsed_seconds=time.monotonic()-start,modes=modes,multipoles=multipoles) if args.flux_lmax: result['flux_check']=flux_check(r,multipoles,args.flux_lmax) print('Flux balance relative difference',result['flux_check']['relative_difference'],flush=True) result['elapsed_seconds']=time.monotonic()-start args.output.write_text(json.dumps(result,indent=2)+'\n') print('Saved',args.output,flush=True) if __name__=='__main__': main()