Astro-Cosmology

 view release on metacpan or  search on metacpan

Cosmology.pm  view on Meta::CPAN

If you prefer the C<func($object,...)> style, then you need to
import the functions:

  use Astro::Cosmology qw( :Func );

Most functions have two names; a short one and a (hopefully) more
descriptive one, such as C<pmot_dist()> and C<proper_motion_distance()>.

Most of the routines below include a C<sig:> line in their documentation.
This is an attempt to say how they
`L<thread|PDL::indexing>' (in the L<PDL|PDL> sense of the word).
So, for routines like C<lum_dist> - which have a sig line of
C<dl() = $cosmo-E<gt>lum_dist( z() )> - the return value has the
same format as the input C<$z> value; supply a scalar, get a scalar back,
send in a piddle and get a piddle of the same dimensions back.
For routines like C<abs_mag> - with a sig line of
C<absmag() = $cosmo-E<gt>abs_mag( appmag(), z() )> - you can thread over
either of the two input values,
in this case the apparent magnitude and redshift.

=head1 SUMMARY

=head2 Utility routines

=over 4

=item * new

Cosmology.pm  view on Meta::CPAN

# note:
#  one parameter:  age of z = 0  -> $1
#  two parameters: age of z = $1 -> $2
#
# ie lookback_time ( [ $z_low ], $z_high )
#
# note: we call the C code directly
#
sub lookback_time ($$;$) {
    my $self   = shift;
    my $z_low  = $#_ == 1 ? shift : 0;   # let PDL do the threading if z_high is a piddle
    my $z_high = shift;

    return Astro::Cosmology::Internal::_lookback_time( $z_low, $z_high,
			   $self->{OMEGA_MATTER}, $self->{OMEGA_LAMBDA},
			   $self->{ABSTOL}, $self->{TCONV} );

} # sub: lookback_time

####################################################################
#

Internal/internal.pd  view on Meta::CPAN

# not convinced what the correct thing to do with OtherPars and float/double
# parameters
#
#my $opar = "$TFD(float,double)";  # NO - unsurprisingly, OtherPars is not checked for macros
#my $opar = "double";
my $opar = "float";

#
# the luminosity distance calculation is just a big switch statement
# to select the cosmology (defined by the p_method parameter)
# and then a threadloop on the redshifts to calculate
# the results
#
pp_def( '_lum_dist',
        Doc => undef,
        Pars => 'z(); [o] d();',
        OtherPars => "int p_method; $opar p_om; $opar p_ol; $opar p_kappa; $opar p_abstol; $opar p_convert;",
        GenericTypes => ['D', 'F'],
        Code => '$GENERIC() zz, zzp1;
                 $GENERIC() int_fn, int_fn_err;

                 $GENERIC() om      = $COMP( p_om );
                 $GENERIC() ol      = $COMP( p_ol );
                 $GENERIC() kappa   = $COMP( p_kappa );
                 $GENERIC() abstol  = $COMP( p_abstol );
                 $GENERIC() convert = $COMP( p_convert );

                 int method  = $COMP( p_method );

                 switch( method ) {
                   case EMPTY:
                     threadloop %{
                        zz   = $z();
                        $d() = 0.5 * convert * zz * ( 2.0 +  zz );
                     %}
                     break;

                   case MATTER_FLAT:
                     threadloop %{
                        zzp1 = 1.0 + $z();
                        $d() = 2.0 * convert * ( zzp1 - sqrt( zzp1 ) );
                     %}
                     break;

                   case MATTER_OPEN:
                   case MATTER_CLOSED:
                     threadloop %{
                        zz   = $z();
                        $d() = 2.0 * convert *
                              ( om * zz + (om - 2.0) * ( sqrt( 1.0 + om * zz ) - 1.0 ) ) /
                              ( om * om );
                     %}
                     break;

                   case LAMBDA_CLOSED:
                     threadloop %{
                        zz   = $z();
                        zzp1 = 1.0 + zz;
                        if ( ! $TFD(romberg_f,romberg_d) (
                                 $TFD(comov_dist_fn_f,comov_dist_fn_d),
                                 om, ol, 0.0, zz, abstol, &int_fn, &int_fn_err ) ) {
			    croak( "Unable to integrate comoving distance.\n" );
			}
                        $d() = convert * zzp1 * sin( kappa * int_fn ) / kappa;
                     %}
                     break;

                   case LAMBDA_OPEN:
                     threadloop %{
                        zz   = $z();
                        zzp1 = 1.0 + zz;
                        if ( ! $TFD(romberg_f,romberg_d) (
                                 $TFD(comov_dist_fn_f,comov_dist_fn_d),
                                 om, ol, 0.0, zz, abstol, &int_fn, &int_fn_err ) ) {
			    croak( "Unable to integrate comoving distance.\n" );
			}
                        $d() = convert * zzp1 * sinh( kappa * int_fn ) / kappa;
                     %}
                     break;

                   case LAMBDA_FLAT:
                     threadloop %{
                        zz   = $z();
                        zzp1 = 1.0 + zz;
                        if ( ! $TFD(romberg_f,romberg_d) (
                                 $TFD(comov_dist_fn_f,comov_dist_fn_d),
                                 om, ol, 0.0, zz, abstol, &int_fn, &int_fn_err ) ) {
			    croak( "Unable to integrate comoving distance.\n" );
			}
                        $d() = convert * zzp1 * int_fn;
                     %}
                     break;

Internal/internal.pd  view on Meta::CPAN

        Pars => 'z(); [o] d();',
        OtherPars => "$opar p_om; $opar p_ol; $opar p_abstol; $opar p_convert;",
        GenericTypes => ['D', 'F'],
        Code => '$GENERIC() int_fn, int_fn_err;

                 $GENERIC() om      = $COMP( p_om );
                 $GENERIC() ol      = $COMP( p_ol );
                 $GENERIC() abstol  = $COMP( p_abstol );
                 $GENERIC() convert = $COMP( p_convert );

                 threadloop %{
                    if ( ! $TFD(romberg_f,romberg_d) (
                             $TFD(comov_dist_fn_f,comov_dist_fn_d),
                             om, ol, 0.0, $z(), abstol, &int_fn, &int_fn_err ) ) {
		        croak( "Unable to integrate comoving distance.\n" );
		    }
                    $d() = convert * int_fn;
                 %}'
        );

# the differential proper motion distance stuff is currently only used

Internal/internal.pd  view on Meta::CPAN

        Pars => 'z(); dm(); [o] answer();',
        OtherPars => "$opar p_cosmo; $opar p_om; $opar p_ol; $opar p_ok;",
        GenericTypes => ['D', 'F'],
        Code => '$GENERIC() cosmo   = $COMP( p_cosmo );
                 $GENERIC() om      = $COMP( p_om );
                 $GENERIC() ol      = $COMP( p_ol );
                 $GENERIC() absok   = fabs( $COMP( p_ok ) );

                 $GENERIC() ddm;

                 threadloop %{
                    $answer() = $TFD(comov_dist_fn_f,comov_dist_fn_d) ( $z(), om, ol );
                 %}

                 if ( cosmo != MATTER_FLAT && cosmo != LAMBDA_FLAT  ) {
                    threadloop %{
                       ddm        = $dm();
                       $answer() *= sqrt( 1.0 + absok * ddm * ddm );
                    %}
                 }'
        );

####################################################################
#
# volume measures: PP code
#
####################################################################

#
# the comoving volume calculation is just a big switch statement
# to select the cosmology (defined by the p_method parameter)
# and then a threadloop on the proper-motion distance (which is
# in units of the Hubble distance) to calculate the volumes
#
pp_def( '_comov_vol',
     Doc => undef,
     Pars => 'dm(); [o] v();',
     OtherPars => "int p_method; $opar p_om; $opar p_ol; $opar p_ok; $opar p_kappa; $opar p_convert;",
     GenericTypes => ['D', 'F'],
     Code => '$GENERIC() ddm;

              $GENERIC() om      = $COMP( p_om );

Internal/internal.pd  view on Meta::CPAN

              $GENERIC() convert = $COMP( p_convert );

              int    method  = $COMP( p_method );

              /* NOTE: here dm is in units of the hubble distance */

              switch( method ) {
                case MATTER_FLAT:
                case LAMBDA_FLAT:
                  convert /= 3.0;
                  threadloop %{
                     ddm  = $dm();
                     $v() = convert * ddm * ddm * ddm;
                  %}
                  break;

                case EMPTY:
                case MATTER_OPEN:
                case LAMBDA_OPEN:
                  convert /= ( 2.0 * ok );
                  threadloop %{
                     ddm  = $dm();
                     $v() = convert *
                            ( ddm * sqrt(1.0 + ok*ddm*ddm) - asinh( ddm * kappa ) / kappa );
                  %}
                  break;

                case MATTER_CLOSED:
                case LAMBDA_CLOSED:
                  convert /= ( 2.0 * ok );
                  threadloop %{
                     ddm  = $dm();
                     $v() = convert *
                            ( ddm * sqrt(1.0 + ok*ddm*ddm) - asin( ddm * kappa ) / kappa );
                  %}
                  break;

                default:
                  croak( "Currently can not handle method == %2d\n", method ); /* should not occur */

              } /* end of switch */'

Internal/internal.pd  view on Meta::CPAN

pp_def( '_dcomov_vol',
     Doc => undef,
     Pars => 'dm(); ddmdz(); [o] dv();',
     OtherPars => "$opar p_ok; $opar p_convert;",
     GenericTypes => ['D', 'F'],
     Code => '$GENERIC() ddm;

              $GENERIC() ok      = $COMP( p_ok );
              $GENERIC() convert = $COMP( p_convert );

              threadloop %{
                 ddm   = $dm();
                 $dv() = convert * ddm * ddm * $ddmdz() / sqrt( 1.0 + ok*ddm*ddm );
              %}'
     );

####################################################################
#
# time measures: PP code
#
####################################################################

Internal/internal.pd  view on Meta::CPAN

        OtherPars => "$opar p_om; $opar p_ol; $opar p_abstol; $opar p_convert;",
        GenericTypes => ['D', 'F'],
        Code => '
                 $GENERIC() om      = $COMP( p_om );
                 $GENERIC() ol      = $COMP( p_ol );
                 $GENERIC() abstol  = $COMP( p_abstol );
                 $GENERIC() convert = $COMP( p_convert );

                 $GENERIC() int_fn, int_fn_err;

                 threadloop %{
                    if ( ! $TFD(romberg_f,romberg_d) (
                                $TFD(lookback_time_fn_f,lookback_time_fn_d),
                                om, ol, $zl(), $zh(), abstol, &int_fn, &int_fn_err ) ) {
                       croak( "ERROR: Unable to integrate lookback time.\n" );
                    }
                    $t() = convert * int_fn;
                 %}'
        );

# end



( run in 0.743 second using v1.01-cache-2.11-cpan-3cd7ad12f66 )