c=========================================================== c bsidpa: Uses LSODA to integrate ODEs which define c c=========================================================== subroutine bsidpa(a,alpha,phi,pp,m,r,p,phi0,rmax, & lsoda_tol,w_tol,dr,whi,wlo,om,n,pn,pt,dft,tail_type) implicit none external flag c----------------------------------------------------------- c Command-line arguments: phi0, rmax, lsoda_tol, c w_tol, dr, whi, wlo c----------------------------------------------------------- real*8 phi0 real*8 rmax real*8 lsoda_tol real*8 w_tol real*8 dr real*8 whi real*8 wlo 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----------------------------------------------------------- real*8 n real*8 pcents parameter ( pcents = 1.0d-2 ) integer i logical success real*8 rtbis c----------------------------------------------------------- c Parameters default values (local) c----------------------------------------------------------- real*8 tol_dft parameter ( tol_dft = 1.0d-10 ) real*8 w_tol_dft parameter ( w_tol_dft = 1.0d-10 ) real*8 drout_dft parameter ( drout_dft = 1.0d-4 ) real*8 whi_dft parameter ( whi_dft = 1.0d0 ) real*8 wlo_dft parameter ( wlo_dft = 1.0d0 + whi_dft*pcents) integer tail_tp_dft parameter ( tail_tp_dft = 1 ) c----------------------------------------------------------- c Stores the initial values of the user optional c parameters (local) c----------------------------------------------------------- real*8 tol_init real*8 w_tol_init real*8 drout_init real*8 whi_init real*8 wlo_init integer tail_tp_init c----------------------------------------------------------- c Define the potential type and the vector of adimensional c parameters. c----------------------------------------------------------- integer pn integer pt real*8 p(pn) c----------------------------------------------------------- c Fields and geometric functions used in spherical symmetric c case. m is the derived ADM mass c----------------------------------------------------------- real*8 phi(n), pp(n), alpha(n), a(n), & r(n), m(n) c----------------------------------------------------------- c Common communication with routine 'fcn' in 'fcn.f'. c----------------------------------------------------------- include 'fcn.inc' c----------------------------------------------------------- c Stores the w value to be communicated outside the subroutine c----------------------------------------------------------- real*8 om c----------------------------------------------------------- c Debug auxiliary variables c----------------------------------------------------------- character*(6) cdnm parameter ( cdnm = 'bsidpa' ) logical ltracet parameter ( ltracet = .true. ) logical ltracef parameter ( ltracef = .false. ) c----------------------------------------------------------- c Assign the values of the external parameters parsed @ c bsidpa_drv to the values of the internal parameters c to be stored in a common block (@ fcn.inc) and used by c a variaty of other subroutines. c----------------------------------------------------------- domain_factor = 3.0d0 ph0 = phi0 rmx = domain_factor * rmax drout = dr tol = lsoda_tol pty = pt pdim = pn tail_tp = tail_type w = om do i = 1, pdim pvec(i) = p(i) enddo pie = 4.0d0 * atan(1.0d0) if ( ltracef ) then write (0,*) cdnm, ': arguments saved in the common block sofar' write (0,*) cdnm, ': ph0 = ', ph0 write (0,*) cdnm, ': rmx = ', rmx write (0,*) cdnm, ': drout = ', drout write (0,*) cdnm, ': tol = ', tol write (0,*) cdnm, ': pty = ', pty write (0,*) cdnm, ': pdim = ', pdim write (0,*) cdnm, ': pvec(i) = ', (pvec(i),i=1,pdim) write (0,*) cdnm, ': pie = ', pie write (0,*) cdnm, ': dft = ', dft write (0,*) cdnm, ': tail_tp = ', tail_tp write (0,*) end if if (dft .eq. 0)then tol_init = tol_dft w_tol_init = w_tol_dft drout_init = drout_dft whi_init = whi_dft wlo_init = wlo_dft tail_tp_init = tail_tp_dft write (0,*) cdnm, ': WARNING! Using default value for lsoda_tol' & //': ', tol_dft write (0,*) cdnm, ': WARNING! Using default value for w_tol' & //': ', w_tol_dft write (0,*) cdnm, ': WARNING! Using default value for dr' & //': ', drout_dft write (0,*) cdnm, ': WARNING! Using default value for whi' & //': ', whi_dft write (0,*) cdnm, ': WARNING! Using default value for wlo' & //': ', wlo_dft write (0,*) cdnm, ': WARNING! Using default value for tail_type' & //': ', tail_tp_dft else if (dft.eq.1)then tol_init = tol w_tol_init = w_tol drout_init = drout_dft whi_init = whi_dft wlo_init = wlo_dft tail_tp_init = tail_tp_dft write (0,*) cdnm, ': WARNING! Using default value for dr' & //': ', drout_dft write (0,*) cdnm, ': WARNING! Using default value for whi' & //': ', whi_dft write (0,*) cdnm, ': WARNING! Using default value for wlo' & //': ', wlo_dft write (0,*) cdnm, ': WARNING! Using default value for tail_type' & //': ', tail_tp_dft else if (dft.eq.2)then tol_init = tol w_tol_init = w_tol drout_init = drout whi_init = whi wlo_init = wlo tail_tp_init = tail_tp else write(0,*)cdnm,'Invalid value for dft!!!' endif tol = tol_init w_tol = w_tol_init drout = drout_init whi = whi_init wlo = wlo_init tail_tp = tail_tp_init if ( ltracef ) then write (0,*) cdnm, ': Default values assigned?' write (0,*) cdnm, ': ph0 = ', ph0 write (0,*) cdnm, ': n = ', n write (0,*) cdnm, ': rmx = ', rmx write (0,*) cdnm, ': tol = ', tol write (0,*) cdnm, ': w_tol = ', w_tol write (0,*) cdnm, ': drout = ', drout write (0,*) cdnm, ': whi = ', whi write (0,*) cdnm, ': wlo = ', wlo write (0,*) cdnm, ': dft = ', dft write (0,*) cdnm, ': tail_tp = ', tail_tp write (0,*) end if success = .false. call first_brac(flag,whi,wlo,success) if (.not.success) then write(0,*)cdnm,': ' write(0,*)'Default values for whi: ', whi_dft write(0,*)'Default values for wlo: ', wlo_dft write(0,*)'Initial values for whi: ', whi_init write(0,*)'Initial values for wlo: ', wlo_init write(0,*)'Current values for whi: ', whi write(0,*)'Current values for wlo: ', wlo stop endif w = rtbis(flag,whi,wlo,w_tol) write(0,*)cdnm,': Current value for w: ', w call tail_polareal(a,alpha,phi,pp,m,r,n) call polareal(a,alpha,phi,pp,m,r,n) call rescale(a,alpha,n) om = w return end c=========================================================== c Subroutine to find the bracket to be used in binary c search c=========================================================== subroutine first_brac(func,x1,x2,success) implicit none integer ntry, func external func real*8 x1, x2 real*8 mfactor, dfactor parameter (mfactor = 2.0d0, dfactor = 0.5d0, ntry = 50) integer f1, f2, j logical success c----------------------------------------------------------- c Debug auxiliary variables c----------------------------------------------------------- character*(10) cdnm parameter ( cdnm = 'first_brac' ) logical ltracet parameter ( ltracet = .true. ) logical ltracef parameter ( ltracef = .false. ) if(x1.eq.x2) then write(0,*)cdnm,': You have to guess an initial '// & 'range for the bracket in first_brac' stop endif f1 = func(x1) f2 = func(x2) if(f1.ne.1) then write(0,*)cdnm,': The first element of the bracket must give '// & 'a positive value ' write(0,*)' for the function. '// & 'Try a new guess for it!!!' stop endif success = .true. do j = 1, ntry if (f1*f2.lt.0) return if(f2.eq.1)then x2 = x2 + mfactor * (x2 - x1) f2 = func(x2) if(ltracet)then write(0,*)cdnm,': f2 = 1' write(0,*)cdnm,': x1 = ',x1, ', f1 = ',f1 write(0,*)cdnm,': x2 = ',x2, ', f2 = ',f2 endif endif if(f2.eq.0)then x2 = x2 - dfactor * (x2 - x1) f2 = func(x2) if(ltracet)then write(0,*)cdnm,': f2 = 0' write(0,*)cdnm,': x1 = ',x1, ', f1 = ',f1 write(0,*)cdnm,': x2 = ',x2, ', f2 = ',f2 endif endif enddo write(0,*)cdnm,': A bracket for binarying search'// & ' the values of w could not be found!' write(0,*)' after ',ntry,'times ' write(0,*)' Please assign better values '// & 'for whi and wlo!' success = .false. return end c=========================================================== c Using bisection, find the root of a function func known c to lie between x1 and x2. The root, returned as rtbis, c will be refined until its accuracy is +- xacc. c This routine was borrowed, and then modified, from c numerical recipes in Fortran. c=========================================================== real*8 function rtbis(func,x1,x2,xacc) implicit none integer jmax real*8 x1, x2, xacc integer func parameter (jmax = 40) !Maximum allowed number of bisections integer j real*8 dx, f, fmid, xmid c----------------------------------------------------------- c Debug auxiliary variables c----------------------------------------------------------- character*(5) cdnm parameter ( cdnm = 'rtbis' ) logical ltracet parameter ( ltracet = .true. ) logical ltracef parameter ( ltracef = .false. ) fmid = func(x2) f = func(x1) if (f*fmid.ge.0.0d0)then write(0,*)cdnm,': Root must be bracketed in rtbis!' endif if(f.lt.0.0d0)then rtbis=x1 dx=x2-x1 else rtbis=x2 dx=x1-x2 endif do j = 1, jmax dx=dx*0.5d0 xmid=rtbis+dx fmid=func(xmid) if(fmid.le.0.0d0)rtbis=xmid if(abs(dx).lt.xacc.or.fmid.eq.0.0d0)return enddo write(0,*)cdnm,': Too many bysections in rtbis!!!' return end c=========================================================== c flag: 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 := Central amplitude for the scalar field c := Estimate eigenvalue, for any , c there is a single which results in c a solution of the system satisfing the c fields and metric components boundary c conditions c c := Maximum range of integration 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. c := LSODA tolerance parameter. c c=========================================================== integer function flag(om) implicit none c----------------------------------------------------------- c LSODA Variables. c----------------------------------------------------------- external fcn, jac integer neq parameter ( neq = 4 ) real*8 y(neq) real*8 r, rout integer itol real*8 rtol, atol integer itask, istate, iopt integer lrw parameter ( lrw = 22 + neq * 16 ) real*8 rwork(lrw) integer liw parameter ( liw = 20 + neq ) integer iwork(liw) integer jt c----------------------------------------------------------- c Common communication with routine 'fcn' in 'fcn.f'. c----------------------------------------------------------- include 'fcn.inc' c----------------------------------------------------------- c Locals. c----------------------------------------------------------- integer last_flag, icount real*8 om c----------------------------------------------------------- c Debug auxiliary variables c----------------------------------------------------------- character*(4) cdnm parameter ( cdnm = 'flag' ) logical ltracet parameter ( ltracet = .true. ) logical ltracef parameter ( ltracef = .false. ) c----------------------------------------------------------- c Set LSODA parameters. Use same value for absolute c and relative tolerance. c----------------------------------------------------------- itol = 1 rtol = tol atol = tol itask = 1 iopt = 0 jt = 2 c----------------------------------------------------------- c Passes the argument of flag (local but external variable) c to the internal variable w, i.e. specified in the fcn.inc c and shared be several procedures. c----------------------------------------------------------- w = om c----------------------------------------------------------- c Initialize the solution, and output it. c----------------------------------------------------------- c r = 1.0d-8 r = 0.0 y(1) = 1.0d0 ! a(0) = 1 - Regularity condition y(2) = 1.0d0 ! alpha(0) = 1 - Rescale along with w later y(3) = ph0 ! phi(0) = ph0 - Central field y(4) = 0.0d0 ! pp(0) = 0 - Regularity condition icount = 0 !Counts the number of turns the data does flag = int(sign(1.0d0,y(3))) c----------------------------------------------------------- c Do the integration. c----------------------------------------------------------- c call xsetun (0) istate = 1 do while( r .lt. rmx ) last_flag = flag rout = r + drout call lsoda(fcn,neq,y,r,rout, & itol,rtol,atol,itask, & istate,iopt,rwork,lrw,iwork,liw,jac,jt) if( istate .lt. 0 ) then write(0,*) 'flag: Error return ', istate, & ' from LSODA ' write(0,*) 'flag: Current interval ', & r, r + drout return end if flag = int(sign(1.0d0,y(3))) if (last_flag.ne.flag)then icount=icount+1 if (icount.eq.3) then flag = 0 return endif endif if(abs(y(3)).ge.1.1d0*ph0) return end do return end c=========================================================== c tail_polareal: 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. Fits a tail to the solution c when a threshold is reached. In this code, 3 diffent types c of tails were implemented: c c 1) A/r + B/r^2 + C/r^3 c 2) A exp(-B*r) + C/r c 3) A exp(-B*r) c c=========================================================== subroutine tail_polareal(a,alpha,phi,pp,m,r,n) implicit none c----------------------------------------------------------- c Fields and geometric functions used in spherical symmetric c case. m is the derived ADM mass c----------------------------------------------------------- integer n real*8 phi(n), pp(n), alpha(n), a(n), & r(n), m(n) c----------------------------------------------------------- c LSODA Variables. c----------------------------------------------------------- external fcn, jac integer neq parameter ( neq = 4 ) real*8 y(neq) integer itol real*8 rtol, atol integer itask, istate, iopt integer lrw parameter ( lrw = 22 + neq * 16 ) real*8 rwork(lrw) integer liw parameter ( liw = 20 + neq ) integer iwork(liw) integer jt c----------------------------------------------------------- c Common communication with routine 'fcn' in 'fcn.f'. c----------------------------------------------------------- include 'fcn.inc' c----------------------------------------------------------- c Locals. c----------------------------------------------------------- real*8 pp_prime(n) real*8 h integer i, ifit real*8 rcut logical set_tail integer pos_concav real*8 AA, B, C, D real*8 fc(n), gc(n), hc(n) c----------------------------------------------------------- c To add an exponential tail: c The criteria is the following: when c phi(r) reaches phi_tol and its second derivative c (its concavity) changes from + to - this code adds an c exponential tail to the solution for phi(r). Note that c this criteria doesn't work for all cases. For example c when pphi(r) changes the sign from - to + in the case c the solution blows up at the tail, we could say that the c stellar radius was reached as well. So in order to rely only c on the first criteria, a more robust one, YOU should make c sure that the w value chosen gives a blown down tail! c----------------------------------------------------------- real*8 phi_tol c----------------------------------------------------------- c Debug auxiliary variables c----------------------------------------------------------- character*(13) cdnm parameter ( cdnm = 'tail_polareal' ) logical ltracet parameter ( ltracet = .true. ) logical ltracef parameter ( ltracef = .false. ) phi_tol = 1.0d-4 * ph0 c----------------------------------------------------------- c Defines the grid for the domain [0,rmx]. Note that c domain_factor brings rmx to its original value. It was c changed in order to shoot and get a more precise value c for w. c----------------------------------------------------------- rmx = rmx / domain_factor if ( n.eq.1 ) then write(0,*)'The number of grid points must be ' // & ' different from "1"' stop else h = rmx / ( n - 1 ) end if do i = 1, n r(i) = (i-1)*h enddo r(n) = rmx c----------------------------------------------------------- c Set LSODA parameters. Use same value for absolute c and relative tolerance. c----------------------------------------------------------- itol = 1 rtol = tol atol = tol itask = 1 iopt = 0 jt = 2 c----------------------------------------------------------- c Initialize counters c----------------------------------------------------------- pos_concav=0 !initialize the positive concavity counter set_tail=.false. !when true fit a tail to phi configuration c----------------------------------------------------------- c Initialize the solution. c----------------------------------------------------------- a(1) = 1.0d0 ! a(0) = 1 - Regularity condition alpha(1) = 1.0d0 ! alpha(0) = 1 - Rescale along with w later phi(1) = ph0 ! phi(0) = ph0 - Central field pp(1) = 0.0d0 ! pp(0) = 0 - Regularity condition m(1) = r(1) * 0.5d0 * (1.0d0 - 1.0d0/a(1)**2) ! mass aspect function c----------------------------------------------------------- c Initialize the solution for comunication with lsoda. c----------------------------------------------------------- y(1) = a(1) ! a(0) = 1 - Regularity condition y(2) = alpha(1) ! alpha(0) = 1 - Rescale along with w later y(3) = phi(1) ! phi(0) = ph0 - Central field y(4) = pp(1) ! pp(0) = 0 - Regularity condition c----------------------------------------------------------- c Do the integration. c----------------------------------------------------------- c call xsetun (0) istate = 1 i = 2 do while(.not.set_tail) call lsoda(fcn,neq,y,r(i-1),r(i), & itol,rtol,atol,itask, & istate,iopt,rwork,lrw,iwork,liw,jac,jt) if( istate .lt. 0 ) then write(0,*) 'tail_polareal: Error return ', istate, & ' from LSODA ' write(0,*) 'tail_polareal: Current interval ', & r(i), r(i+1) return end if a(i) = y(1) alpha(i) = y(2) phi(i) = y(3) pp(i) = y(4) m(i) = r(i) * 0.5d0 * (1.0d0 - 1.0d0/a(i)**2) c----------------------------------------------------------- c The last part of this unit of code test if necessary c and fits a tail to the solution for phi and pp c----------------------------------------------------------- pp_prime(i) = 0.5d0*(3.0d0*pp(i) - 4.0d0*pp(i-1) + pp(i-2) )/ h c----------------------------------------------------------- c It finds rcut and calculate A and B to be used in the c exponencial fitting the tail of the scalar field, c T(r)=A exp (-B r), for example. c rcut correspond to the position when phi'' changes sign c from + to - and/or phi(r) goes below a prescribed tolerance. c----------------------------------------------------------- if(int(sign(1.0d0,pp_prime(i))).eq.1.and.(.not.set_tail) ) then pos_concav=1 endif c write(*,*)'phi sign = ', int(sign(1.0d0,y(3))) c write(*,*)'pp_prime sign = ', int(sign(1.0d0,pp_prime)) if((int(sign(1.0d0,pp_prime(i))).eq.-1.or.(phi(i).lt.phi_tol)) & .and. pos_concav.eq.1) then set_tail=.true. ifit = i write(*,*)cdnm,': ifit = ', ifit endif i = i + 1 end do c----------------------------------------------------------- c lsoda rewrite the value of r(i-1) at r(i). To avoid c having r(ifit-1)=r(ifit), values of r(i) are reassigned c before another calculation. c----------------------------------------------------------- do i = 1, n r(i) = (i-1)*h enddo r(n) = rmx write(*,*)cdnm, ': Should be printed at most once!!!!!!' if(tail_tp.eq.1)then write(*,*)cdnm, ': Tail to be fit: A/r + B/r^2 + C/r^3' B=-3.0d0*r(ifit-1)**2*phi(ifit-1)-5.0d0*r(ifit-1)**3 & *pp(ifit-1)-r(ifit-1)**4*pp_prime(ifit-1) AA=0.5d0*(r(ifit-1)**2*pp(ifit-1)+3.0d0*r(ifit-1)*phi(ifit-1) & -B/r(ifit-1)) C=r(ifit-1)**3*phi(ifit-1)-AA*r(ifit-1)**2-B*r(ifit-1) if(AA.lt.0.0d0.or.B.lt.0.0d0.or.C.lt.0.0d0)then write(*,*)cdnm, ': WARNING! Tail choice does not fit properly' tail_tp = 2 endif endif if (tail_tp.eq.2)then write(*,*)cdnm, ': Tail to be fit: A exp(-Br) + C/r' fc(ifit-1)=r(ifit-1)**2*pp_prime(ifit-1)-2.0d0*phi(ifit-1) gc(ifit-1)=r(ifit-1)*pp(ifit-1)+phi(ifit-1) hc(ifit-1)=fc(ifit-1)/gc(ifit-1) B=0.5d0/r(ifit-1)**2*(-r(ifit-1)*hc(ifit-1)+dsqrt(r(ifit-1)**2 & *hc(ifit-1)**2+4.0d0*r(ifit-1)**2*(2.0d0+hc(ifit-1)))) AA=gc(ifit-1)*dexp(B*r(ifit-1))/(1.0d0-B*r(ifit-1)) C=r(ifit-1)*(phi(ifit-1)-AA*dexp(-B*r(ifit-1))) if(AA.lt.0.0d0.or.B.lt.0.0d0.or.C.lt.0.0d0)then write(*,*)cdnm, ': WARNING! Tail choice does not fit properly' tail_tp = 3 endif endif if(tail_tp.eq.3)then write(*,*)cdnm, ': Tail to be fit: A exp(-Br) ' B=-pp(ifit-1)/phi(ifit-1) AA=phi(ifit-1)*dexp(B*r(ifit-1)) if(AA.lt.0.0d0.or.B.lt.0.0d0)then write(*,*)cdnm, ': WARNING! Tail choice does not fit properly' tail_tp = 3 endif endif write(*,*)cdnm,': tail_type = ',tail_tp write(*,*)cdnm,': AA= ',AA,'B= ',B,'C = ',C,'phi= ' & ,phi(ifit),'pp = ',pp(ifit) do i = ifit , n if(tail_tp.eq.1)then phi(i)=AA/r(i)+B/r(i)**2+C/r(i)**3 pp(i)=-AA/r(i)**2-2.0d0*B/r(i)**3-3.0d0*C/r(i)**4 else if (tail_tp.eq.2)then phi(i)=AA*dexp(-B*r(i))+C/r(i) pp(i)=-AA*B*dexp(-B*r(i))-C/r(i)**2 else phi(i)=AA*dexp(-B*r(i)) pp(i)=-B*phi(i) endif enddo return end c=========================================================== c polareal: 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 Integrates only the equations for a and alpha. It uses c the values for phi and pp with a tail fit. c c=========================================================== subroutine polareal(a,alpha,phi,pp,m,r,n) implicit none c----------------------------------------------------------- c Fields and geometric functions used in spherical symmetric c case. m is the derived ADM mass c----------------------------------------------------------- integer n real*8 phi(n), pp(n), alpha(n), a(n), & r(n), m(n) c----------------------------------------------------------- c LSODA Variables. c----------------------------------------------------------- external fcn, jac integer neq parameter ( neq = 2 ) real*8 y(neq) integer itol real*8 rtol, atol integer itask, istate, iopt integer lrw parameter ( lrw = 22 + neq * 16 ) real*8 rwork(lrw) integer liw parameter ( liw = 20 + neq ) integer iwork(liw) integer jt c----------------------------------------------------------- c Common communication with routine 'fcn' in 'fcn.f'. c----------------------------------------------------------- include 'fcn.inc' c----------------------------------------------------------- c Locals. c----------------------------------------------------------- integer i c----------------------------------------------------------- c Debug auxiliary variables c----------------------------------------------------------- character*(8) cdnm parameter ( cdnm = 'polareal' ) logical ltracet parameter ( ltracet = .true. ) logical ltracef parameter ( ltracef = .false. ) c----------------------------------------------------------- c Set LSODA parameters. Use same value for absolute c and relative tolerance. c----------------------------------------------------------- itol = 1 rtol = tol atol = tol itask = 1 iopt = 0 jt = 2 c----------------------------------------------------------- c Initialize the solution for comunication with lsoda. c----------------------------------------------------------- y(1) = a(1) ! a(0) = 1 - Regularity condition y(2) = alpha(1) ! alpha(0) = 1 - Rescale along with w later phj = phi(1) ! value of central field at origin ppj = pp(1) ! and its first derivative in r m(1) = r(1) * 0.5d0 * (1.0d0 - 1.0d0/a(1)**2) ! mass aspect function c----------------------------------------------------------- c Do the integration. c----------------------------------------------------------- c call xsetun (0) istate = 1 do i = 2, n call lsoda(fcn,neq,y,r(i-1),r(i), & itol,rtol,atol,itask, & istate,iopt,rwork,lrw,iwork,liw,jac,jt) if( istate .lt. 0 ) then write(0,*) 'polareal: Error return ', istate, & ' from LSODA ' write(0,*) 'polareal: Current interval ', & r(i), r(i+1) return end if a(i) = y(1) alpha(i) = y(2) phj = phi(i) ppj = pp(i) m(i) = r(i) * 0.5d0 * (1.0d0 - 1.0d0/a(i)**2) end do return end c----------------------------------------------------------- c rescale: rescales the values of alpha and w c----------------------------------------------------------- subroutine rescale(a,alpha,n) implicit none c----------------------------------------------------------- c Geometric functions used in spherical symmetric c case. c----------------------------------------------------------- integer n real*8 alpha(n), a(n) c----------------------------------------------------------- c Common communication with routine 'fcn' in 'fcn.f'. c----------------------------------------------------------- include 'fcn.inc' c----------------------------------------------------------- c Locals. c----------------------------------------------------------- integer i real*8 c c----------------------------------------------------------- c Debug auxiliary variables c----------------------------------------------------------- character*(7) cdnm parameter ( cdnm = 'rescale' ) logical ltracet parameter ( ltracet = .true. ) logical ltracef parameter ( ltracef = .false. ) c = 1.0d0 / (alpha(n)*a(n)) w = c * w do i = 1 , n alpha(i) = c * alpha(i) enddo return end