c=========================================================== c bsidpa_drv: Driver routine for bsidpa. It parses the c inline arguments and sets all the parametes necessary c for bsidpa. After bsidpa, it output the solution c (a, alpha, phi, pp, pi1, pi2, m, z) to .sdf file. Its c executable is called 'bsidpa'. c c bsidpa: Uses LSODA to integrate ODEs which define c the initial data for a spherically symmetric (1d), c complex, scalar field with self-interaction potential c in polar-areal coordinates. c c usage: bsidpa [ c ] c c or c c bsidpa c [ ] c c Tail types encoded in this function: c 1 ----> A/r + B/r^2 + C/r^3 c 2 ----> A exp(-Br) + C/r c 3 ----> A exp(-Br) c c c := Central amplitude for the scalar field c := Sets the level of solution resolution: c (#of grid points) n = 2^level + 1 c := Maximum range of integration domain. c := LSODA tolerance parameter. c := Tolerance for estimate the eigenvalue c w, coming from the ansatz exp(iwt) for c static solutions. Therefore, for any , c there is a single w0 +- w_tol which results c in a solution of the system satisfing the c fields and metric components boundary c conditions. c c := Output step. We use this as a c parameter rather than an output c level for ease in extending c integrations to larger as c eigenvalue E is determined more c precisely as well as for avoiding LSODA c crashes in the bracketing phase of the c code. c := Sets an initial guess for the bracket [whi,wlo] c Note that usually whi < wlo. whi stands for the c values of w that has its tail to blow up. c Analogously, the tail blows down for wlo. c := The other limit of the bracket. c c Note that the hard bracket [ ] encloses the parameters that are c optional. Therefore, the only really mandatory in line parameters c are , and . c c=========================================================== program bsidpa_drv implicit none integer iargc integer i4arg real*8 r8arg character*2 itoc real*8 r8_never parameter ( r8_never = -1.0d-60 ) c----------------------------------------------------------- c Command-line arguments: phi0, level, rmax, lsoda_tol, c w_tol, dr, whi, wlo c----------------------------------------------------------- real*8 phi0 integer level real*8 rmax real*8 lsoda_tol real*8 w_tol real*8 dr real*8 whi real*8 wlo real*8 w c----------------------------------------------------------- c Flag to let the function know if and what parameters c default should be used. dft = 0 assigns lsoda_tol, w_tol c dr, whi and wlo their the default values. dft = 1 assigns c only dr, whi and wlo their default values. dft = 2 assigns c no default value, i.e. the user must specify all the c relevant parameters values. c----------------------------------------------------------- integer dft c----------------------------------------------------------- c Tail type determines the type of tail to be matched c as the boson star "atmosphere". When phi(r) reaches c a cutoff point or when the concavity of the boson c star tail changes its sign (remember: the shooting c process to find the boson star solution ends up blowing c up or down. This code makes sure the last trial always c blows down) a tail is matched at the cutoff point until c the end of the prescribed domain. Several different tails c are encoded. c c 1) A/r + B/r^2 + C/r^3 c 2) A exp(-Br) + C/r c 3) A exp(-Br) c c The default tail is 1. However I have detected for poor c resolutions a mismatch of the smoother tails. A warning c message is output and other tail tried. The less smooth c one, 4, should work for all cases. To insist in one c particular tail, you ought to increase the resolution c of the boson star solution sought. c----------------------------------------------------------- integer tail_type c----------------------------------------------------------- c Locals. c----------------------------------------------------------- integer n, i c----------------------------------------------------------- c Define the potential type and the vector of adimensional c parameters. c----------------------------------------------------------- integer pn integer pt real*8 p(100) c----------------------------------------------------------- c Fields and geometric functions used in spherical symmetric c case. m is the derived ADM mass c----------------------------------------------------------- integer nmax parameter ( nmax = 2**16 + 1) real*8 phi(nmax), pp(nmax), alpha(nmax), a(nmax), & r(nmax), m(nmax), z(nmax) real*8 pi1(nmax), pi2(nmax) c----------------------------------------------------------- c bbhutil output related variables. c 'gft' output routine:'gft_out_bbox': bbox = bounding c box which specifies the limits of the grid function's c domain. c c bbox := ( xmin, xmax ) c----------------------------------------------------------- real*8 time real*8 bbox(2) integer shape(1) integer rank c----------------------------------------------------------- c Debug auxiliary variables c----------------------------------------------------------- character*(10) cdnm parameter ( cdnm = 'bsidpa_drv' ) logical ltracet parameter ( ltracet = .true. ) logical ltracef parameter ( ltracef = .false. ) c----------------------------------------------------------- c Parse command line arguments. c----------------------------------------------------------- if( iargc() .lt. 3 ) go to 900 if( iargc() .eq. 3 ) dft = 0 if( iargc() .eq. 4 ) go to 900 if( iargc() .eq. 5 ) dft = 1 if( iargc() .eq. 6 .or. iargc() .eq. 7 .or. iargc() .eq. 8) & go to 900 if( iargc() .eq. 9 ) dft = 2 if( iargc() .ge. 10 ) go to 900 phi0 = r8arg(1,r8_never) level = i4arg(2,r8_never) rmax = r8arg(3,r8_never) lsoda_tol = r8arg(4,r8_never) w_tol = r8arg(5,r8_never) dr = r8arg(6,r8_never) whi = r8arg(7,r8_never) wlo = r8arg(8,r8_never) tail_type = i4arg(9,r8_never) if( phi0 .eq. r8_never .or. level .eq. r8_never .or. & rmax .eq. r8_never ) go to 900 if ( ltracef ) then write (0,*) cdnm, ': command line arguments parsed' write (0,*) cdnm, ': phi0 = ', phi0 write (0,*) cdnm, ': level = ', level write (0,*) cdnm, ': rmax = ', rmax write (0,*) cdnm, ': lsoda_tol = ', lsoda_tol write (0,*) cdnm, ': w_tol = ', w_tol write (0,*) cdnm, ': dr = ', dr write (0,*) cdnm, ': whi = ', whi write (0,*) cdnm, ': wlo = ', wlo write (0,*) cdnm, ': dft = ', dft write (0,*) cdnm, ': tail_type = ', tail_type write (0,*) end if c----------------------------------------------------------- c Set up the parameter vector for the self-interaction c potential and the potential type. pt = 0 is massless c einstein-klein-gordon system. pt = 1 is a polynomial c self-interacting potential for the ekg system (default) c----------------------------------------------------------- pt = 1 pn = 4 p(1)=0.0d0 p(2)=1.0d0 p(3)=0.0d0 p(4)=0.0d0 c----------------------------------------------------------- c Set up the resolution of the grid. c----------------------------------------------------------- n = 2**level + 1 c----------------------------------------------------------- c Shoot the value of w in order to find a 1d boson star c solution. Returns the solution for a specified level. c----------------------------------------------------------- call bsidpa(a,alpha,phi,pp,m,r,p,phi0,rmax, & lsoda_tol,w_tol,dr,whi,wlo,w,n,pn,pt,dft,tail_type) write(0,*)cdnm,': Current value for w: ', w do i = 1 , n pi1(i) = 0.0d0 pi2(i) = - w * phi0 * a(i) / alpha(i) z(i) = 2.0d0 * m(i) / r(i) enddo c----------------------------------------------------------- c Output the solution. c----------------------------------------------------------- time = 0.0d0 bbox(1) = 0.0d0 bbox(2) = rmax shape(1) = n rank = 1 call gft_out_bbox('a_'//itoc(level),time,shape,rank,bbox,a) call gft_out_bbox('alpha_'//itoc(level),time,shape,rank,bbox, & alpha) call gft_out_bbox('phi_'//itoc(level),time,shape,rank,bbox,phi) call gft_out_bbox('pp_'//itoc(level),time,shape,rank,bbox,pp) call gft_out_bbox('pi1_'//itoc(level),time,shape,rank,bbox,pi1) call gft_out_bbox('pi2_'//itoc(level),time,shape,rank,bbox,pi2) call gft_out_bbox('m_'//itoc(level),time,shape,rank,bbox,m) call gft_out_bbox('z_'//itoc(level),time,shape,rank,bbox,z) stop 900 continue write(0,*) 'usage: bsidpa '// & ' [ '// & ' ] ' write(0,*) '' write(0,*) ' or' write(0,*) '' write(0,*) ' bsidpa '// & ' '// & ' [ ] ' write(0,*) '' write(0,*) 'Tail types encoded in this function:' write(0,*) ' 1 ----> A/r + B/r^2 + C/r^3' write(0,*) ' 2 ----> A exp(-Br) + C/r' write(0,*) ' 3 ----> A exp(-Br)' write(0,*) '' stop end