光度距離を計算するプログラム（C言語）

光度距離については以下を参照。以下のリンク先にある松原さんの講義ノートと須藤さんの『一般相対論入門』とで共動距離などの定義が微妙に異なってて、松原さんの講義ノートにいたっては章が変わると表記法が変わってたりして、個人的には非常に煩わしかったですが、大事なのは定義の統一ではなくて意味の一貫性ということで。直観的な理解については、ボスのおかげでなんとなくわかった気になっています。

松原さんのサイト > 講義のページ > 過去の講義ページ > 観測的宇宙論 > 講義オリジナルテキスト

各パラメータの値は[km/s/Mpc]としています。

/*--------------------------------------------------------------
 *
 *                      lumi_distance.c
 *
 *--------------------------------------------------------------
 *
 * Description:
 * ===========
 *     光度距離を計算するためのプログラム。天体の赤方偏移を
 *  与えれば、その天体の光度距離の計算結果をターミナルに表示する。
 *
 *    宇宙論パラメータについては次のように定義している。
 *       OMEGA_M  =  0.3
 *       OMEGA_L  =  0.7
 *       H_0      =   70 [km /s /Mpc]
 *
 *    なお、数値積分にはシンプソン則を用いており、プログラム作成に
 *  あたって『C言語による最新アルゴリズム事典』を参考にした。
 *
 *   usage:
 * ===========
 *   usage: lumi_distance  redshift
 *   
 *   redshift  =  天体の赤方偏移z (0 < z < 1000)
 *
 * Revision history:
 * ================
 *   created   December   1  2006  Ono
 *
 *--------------------------------------------------------------
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define  OMEGA_M  0.3 				/* $\Omega_m$ */
#define  OMEGA_L  0.7 				/* $\Omega_\Lambda$ */
#define  H_0      70.0				/* ハッブル定数[km /s /Mpc] */
#define  C        2.99792458e+5	        	/* 光速[km /s] */

double f(double z)  				/* 被積分関数 $f(x)$ */
{
    return C / H_0 /sqrt(OMEGA_M * pow((1+z),3) + OMEGA_L );
}

main(argc, argv)
int  argc;
char *argv[];
{
    int i, n;
    double a = 0, 
	   z,
	   h, 
	   err = 1,
	   epsilon = 1.0e-7,
	   trapezoid, midpoint, simpson,
	   new_simpson,
	   lumi_distance;

    if(argc != 2){
	(void) printf(" usage: lumi_distance  redshift\n");
	(void) exit(-1);
    }

    z = atof(argv[1]);

    if(z <= 0.0 || z >= 1000.0){
	(void) printf("select a redshift 0 < z < 1000.\n");
	(void) exit(-1);
    }

    h = z - a;  
    trapezoid = h * (f(a) + f(z)) / 2;

    for (n = 1; err > epsilon ; n *= 2) {
        midpoint = 0;
        for (i = 1; i <= n; i++) midpoint += f(a + h * (i - 0.5));
        midpoint *= h;
        new_simpson = (trapezoid + 2 * midpoint) / 3;
	err         = fabs(new_simpson - simpson) / fabs(new_simpson);
	simpson     = new_simpson;
        h /= 2;  trapezoid = (trapezoid + midpoint) / 2;
    }

    lumi_distance = (1+z) * simpson;
    printf("value  =  %f [Mpc] \n", lumi_distance);
    return 0;
}