test_sbb_ran_planck.c (4555B)
1 /* Copyright (C) 2018-2023 |Méso|Star> (contact@meso-star.com) 2 * 3 * This program is free software: you can redistribute it and/or modify 4 * it under the terms of the GNU General Public License as published by 5 * the Free Software Foundation, either version 3 of the License, or 6 * (at your option) any later version. 7 * 8 * This program is distributed in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU General Public License for more details. 12 * 13 * You should have received a copy of the GNU General Public License 14 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 15 16 #include "sbb.h" 17 #include <rsys/math.h> 18 #include <math.h> 19 20 static double 21 rand_canonic(void) 22 { 23 return (double)rand() / (double)((size_t)RAND_MAX+1); 24 } 25 26 static void 27 check_planck(void) 28 { 29 struct sbb_ran_planck_create_args args = SBB_RAN_PLANCK_CREATE_ARGS_DEFAULT; 30 struct sbb_ran_planck* planck = NULL; 31 32 args.range[0] = 3.6e-6; 33 args.range[1] = 12.3e-6; 34 args.ref_temperature = 500; 35 args.nbands = 0; 36 37 CHK(sbb_ran_planck_create(NULL, &planck) == RES_BAD_ARG); 38 CHK(sbb_ran_planck_create(&args, NULL) == RES_BAD_ARG); 39 CHK(sbb_ran_planck_create(&args, &planck) == RES_OK); 40 41 CHK(sbb_ran_planck_ref_get(NULL) == RES_BAD_ARG); 42 CHK(sbb_ran_planck_ref_get(planck) == RES_OK); 43 CHK(sbb_ran_planck_ref_put(NULL) == RES_BAD_ARG); 44 CHK(sbb_ran_planck_ref_put(planck) == RES_OK); 45 CHK(sbb_ran_planck_ref_put(planck) == RES_OK); 46 47 args.range[0] = 1; 48 args.range[1] = 0; 49 CHK(sbb_ran_planck_create(&args, &planck) == RES_BAD_ARG); 50 args.range[0] = 0; 51 args.range[1] = 1; 52 args.ref_temperature = -100; 53 CHK(sbb_ran_planck_create(&args, &planck) == RES_BAD_ARG); 54 args.range[0] = 3.6e-6; 55 args.range[1] =-12.3e-6; 56 args.ref_temperature = 300; 57 CHK(sbb_ran_planck_create(&args, &planck) == RES_BAD_ARG); 58 } 59 60 static void 61 planck_integration 62 (const double lambda_lo, /* [m] */ 63 const double lambda_hi, /* [m] */ 64 const double Tref, /* [K] */ 65 const double T0, /* [K] */ 66 const size_t nbands) /* > 0 <=> Speed up Planck sampling */ 67 { 68 struct sbb_ran_planck_create_args args = SBB_RAN_PLANCK_CREATE_ARGS_DEFAULT; 69 struct sbb_ran_planck* planck = NULL; 70 71 const double delta_lambda = lambda_hi - lambda_lo; /* [m] */ 72 73 /* Variables of the Monte Carlo integration */ 74 const size_t N = 10000; /* #realisations */ 75 size_t irealisation = 0; 76 double sum = 0; /* Sum of weights */ 77 double sum2 = 0; /* Sum of weights squared */ 78 double E = 0; /* Expected value */ 79 double V = 0; /* Variance */ 80 double SE = 0; /* Standard Error */ 81 82 double ref = 0; 83 84 /* Setup the planck Distribution */ 85 args.range[0] = lambda_lo; 86 args.range[1] = lambda_hi; 87 args.ref_temperature = Tref; 88 args.nbands = nbands; 89 CHK(sbb_ran_planck_create(&args, &planck) == RES_OK); 90 91 FOR_EACH(irealisation, 0, N) { 92 const double r0 = rand_canonic(); 93 const double r1 = rand_canonic(); 94 double lambda = 0; 95 double pdf = 0; 96 double w = 0; 97 98 lambda = sbb_ran_planck_sample(planck, r0, r1, &pdf); 99 w = sbb_planck_monochromatic(lambda, T0) / pdf; 100 sum += w; 101 sum2 += w*w; 102 } 103 104 E = sum / (double)N; 105 V = MMAX(sum2 / (double)N - E*E, 0); 106 SE = sqrt(V/(double)N); 107 108 ref = sbb_planck(lambda_lo, lambda_hi, T0) * delta_lambda; 109 110 printf( 111 "planck(%g m, %g m, %g K) = %g ~ %g +/- %g [W/m^2/sr] " 112 "(Tref = %g K; nbands = %lu)\n", 113 lambda_lo, lambda_hi, T0, ref, E, SE, Tref, (unsigned long)nbands); 114 CHK(eq_eps(ref, E, 3*SE)); 115 116 CHK(sbb_ran_planck_ref_put(planck) == RES_OK); 117 } 118 119 int 120 main(int argc, char** argv) 121 { 122 /* Input parameters */ 123 double lambda_lo = 0; 124 double lambda_hi = 0; 125 double delta_lambda = 0; 126 size_t nbands = 0; 127 (void)argc, (void)argv; 128 129 check_planck(); 130 131 lambda_lo = 3.6e-6; /* [m] */ 132 lambda_hi = 12.3e-6; /* [m] */ 133 delta_lambda = lambda_hi - lambda_lo; 134 nbands = (size_t)(delta_lambda * 1e9); 135 136 planck_integration(lambda_lo, lambda_hi, 500, 550, 0); 137 planck_integration(lambda_lo, lambda_hi, 500, 550, nbands); 138 planck_integration(lambda_lo, lambda_hi, 300, 550, 0); 139 planck_integration(lambda_lo, lambda_hi, 300, 550, nbands); 140 planck_integration(lambda_lo, lambda_hi, 300, 550, nbands/4); 141 142 lambda_lo = 3.7e-6; /* [m] */ 143 lambda_hi = 4.5e-6; /* [m] */ 144 delta_lambda = lambda_hi - lambda_lo; 145 nbands = (size_t)(delta_lambda * 1e9); 146 147 planck_integration(lambda_lo, lambda_hi, 300, 550, 0); 148 planck_integration(lambda_lo, lambda_hi, 300, 550, nbands); 149 return 0; 150 }