Astro-Cosmology
view release on metacpan or search on metacpan
Cosmology.pm view on Meta::CPAN
225226227228229230231232233234235236237238239240241242243244245246247248249250251If you prefer the C<func(
$object
,...)> style, then you need to
import
the functions:
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
783784785786787788789790791792793794795796797798799800801802803# 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
979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188# 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
200201202203204205206207208209210211212213214215216217218219220
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
225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266
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
270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310$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
313314315316317318319320321322323324325326327328329330331332333pp_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
338339340341342343344345346347348349350351352353354355356357358
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.268 second using v1.01-cache-2.11-cpan-2b0bae70ee8 )