v02:

The purpose of this version was to break up polareal.f in 3 diffent subroutines:

1) tpolareal.f: argument parsing and set up of the parameters used in the problem

2) polareal.f: do the integration using a interface to lsoda called drvlsoda and 
               then output the results.

3) drvlsoda.f: the original idea was to separate any parameter set up and calls 
               to lsoda from the problem to be integrated itself and its output. 
               This way this subroutine would store the lsoda variables and call 
               lsoda to integrate on time step: r to r+dr.


Problems:

It's not possible to separate into drvlsoda.f without compromising simplicity. 
We would end up having to set up lsoda parameters as rwork(lrw) and  iwork(liw)
in polareal.f, breaking this way the original idea of separating lsoda interface
from output manipulation.

Once a parameter is set up in tpolareal.f and passed as argument to polareal.f,
it cannot be included in a common block. An error message appears in the compilation
phase telling that there are "Conflicting attributes or multiple declaration of name".
This problem was solved by giving different names to the parameters set up in tpolareal.f
and those in the common block. Later in polareal.f parameters set up in the 
tpolareal.f should be assigned to those in the common block. The general solution 
for this kind of problem is to give different names to those parameters that will
be used internally in common blocks.

The next version has the intention of implementing the boson star function 
completely in its final version modulo future adjustments.


v03:
        
This is almost the final version for the boson star function in
polar-areal coordinates. The main purpose of this version is to
test the subroutines for bracketing and binary searching as well
as leave a template to study the rate of convergence of these
procedures for a simple test function. Next version, instead of
the function ftest, polareal_brac is actually implemented as well
as polareal.


v04:

This version is almost the final version. The tail for phi and pp wasn't
implemented yet. It comes in the next version. This code serves as basis 
to several diffent implementations for tails, i.e. if you want experiment
with the boson star tail, start writing on top of this code.


v05:

This version fits a 3 parameter fucntion for phi and pp as well as
for a and alpha. It solves a kink it was showing up before when
tried just with an exponential tail, A*exp(-Br), since when one more
parameter is added, A*exp(-Br) + C/r, the second derivative of phi is
also matched at the fitting point. Several trials were done in order to
fit a and alpha as well. Successless! Despite the smoothness of the 
fitting, the control checking mass function, lacks conservation by 5%
from the fitting point on. Next version will reintegrate the equations
for a and alpha in order to solve this problem.



bsidpa 0.06 6 50 1.0e-8 1.0e-8 0.001 1 1.01


v06:

Finally the a and alpha re-integration is successful. However when doing convergence
tests, I realized the tailed was becoming negative, very small though. A more 
close look, showed that most of the time C in A*exp(-Br)+C/r is fitted a negative
number. As the exponential decays faster than 1/r, C/r ends up dominating the tail
behaviour. I played a little with phi_tol in order to avoid being too close of any
turning points before the fitting but couldn't improve the results. Next try will be
with the following tail function:

A/r+B/r^2+C/r^3


v07:

This is the back up for final version for solving the initial value
problem in spherical symmetry. I would say that the tail fitting is 
still an art. I got several times the function parameters fitted such
that the field becomes negative, small though. THe choices of parameters
for the last run seemed to contour this problem. I left 3 possibilities to
implement the tail in case a simple change of parameters doesn't prevent
of the field going negative. Despite this version is considered final for now,
it may be changed in the future in order to gain robustness. I am planning to
use it to find the most compact boson star as well as to prepare initial data 
for the evolution code. It's also going to be the basis for the maximal isotropic
coordinates initial data.

Try:

bsidpa 0.06 10 50 1.0e-10 1.0e-10 0.0001 1 1.01