update: allow -DFLAT=1 to set a flat spacetime without black hole during compilation

src/common.h

@@ -5,18 +5,24 @@
 #include <gsl/gsl_errno.h>
 #include <math.h>
 
+#define POW2(x) ((x)*(x))
+
+#ifndef FLAT
 #define FLAT 0
+#endif
 #define PI 3.1415926535897932384626433832795028841971693993751058
 #define Rs 1
 #define M (FLAT ? 0 : 0.5*Rs)
 #define R0 (15*Rs)
-#define f(r) ((1.0-(2*M)/((double)r)))
+#define f(r) ((1.0-(2*M)/((double)(r))))
 #define bmin (FLAT ? 1e-10 : 1.5*sqrt(3.0)+5e-10)
 #define cotpsi_max (sqrt(R0*R0/(bmin*bmin*f(R0)) - 1))
+#define tanpsi_min (1.0/cotpsi_max)
 #define THETAERROR 100000
 #define color_index(i,j,c) (j*W*3+i*3+c)
 
 #define cutperc (FLAT ? 1.0 : 0.99)
+#define SCALE 1
 
 typedef struct {
   int W;

src/init.c

@@ -63,8 +63,8 @@ static inline double find_root_5ord(double u1, double u2, double v1, double v2,
        pow(u2,2))*v1*v2);
 }
 
-// \chi(cotpsi) is the totol deflection angle with parameter cotpsi=-dr/(r √f dphi)
-double chi(double cotpsi){
+// \chi(tanpsi) is the totol deflection angle with parameter tanpsi=-(r √f dphi)/dr
+double chi(double tanpsi){
   gsl_odeiv2_system sys = {func, jac, 2, NULL};
   gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk8pd, 1e-3, 1e-12, 1e-10);
   
@@ -72,7 +72,7 @@ double chi(double cotpsi){
   double y[2] = {
     u0,
     //sqrt(1.0/(b*b) - u0*u0 + 2.0*M*u0*u0*u0)
-    u0*sqrt(f(R0))*cotpsi
+    u0*sqrt(f(R0))/tanpsi
   };
   double phi = 0;
   double phi_max = 8*PI;
@@ -142,14 +142,13 @@ int init(double rmax, double rela_err_limit, int *size, double **x, double **y){
   sample.size=1;
   sample.x = malloc(sizeof(double) * 1);
   sample.y = malloc(sizeof(double) * 1);
-  double cotpsi_min = 1/rmax;
-  sample.x[0] = cotpsi_min;
-  sample.y[0] = chi(cotpsi_min);
-  double chi_cotpsimax = chi(cotpsi_max);
+  sample.x[0] = tanpsi_min;
+  sample.y[0] = chi(tanpsi_min);
+  double chi_tanpsimax = chi(rmax);
 
-  refine_interval(cotpsi_min, cotpsi_max, sample.y[0], chi_cotpsimax, rela_err_limit, &sample, 1);
+  refine_interval(tanpsi_min, rmax, sample.y[0], chi_tanpsimax, rela_err_limit, &sample, 1);
 
-  sample_push(&sample, cotpsi_max, chi_cotpsimax);
+  sample_push(&sample, rmax, chi_tanpsimax);
 
   //sort
 

src/quick_select.c

@@ -46,7 +46,7 @@ double percentile(const double *arr, size_t n, double perc) {
   }
 
   double *arr_copy = malloc(sizeof(double)*n);
-  memcpy(arr_copy, arr, n);
+  memcpy(arr_copy, arr, n*sizeof(double));
   // nearest-rank: ceil(0.9*n) - 1
   size_t k = (size_t)(perc*n) - 1;
   double ans = select_k(arr_copy, n, k);

src/render.c

@@ -6,8 +6,10 @@
 #include <math.h>
 #include <stdio.h>
 #include "render.h"
+#include "common.h"
 
 typedef struct {
+  int c;
   const System *system;
   const Spline_data *spline_data;
 } Integrand_params;
@@ -19,40 +21,18 @@ inline double xy_to_b(double x, double y) {
   return b;
 }
 
-double integrand_r(double *xy, size_t dim, void *integrand_params_void) {
-  double cotpsi = 1.0/hypot(xy[0], xy[1]);
-  if (cotpsi >= cotpsi_max) return 0;
+double integrand(double *xy, size_t dim, void *integrand_params_void) {
+  double tanpsi = hypot(xy[0], xy[1]);
+  double cotpsi = 1.0/tanpsi;
+  if (tanpsi <= tanpsi_min) return 0;
   Integrand_params *integrand_params = (Integrand_params *)integrand_params_void;
-  double theta = gsl_spline_eval(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc);
+  double theta = gsl_spline_eval(integrand_params->spline_data->spline, tanpsi, integrand_params->spline_data->acc);
   double phi = atan2(xy[1], xy[0]);
   double rgb[3] = {0};
   double theta_phi[2] = {theta, phi};
   integrand_params->system->angle_to_pixel(theta_phi, rgb);
-  return rgb[0]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc));
-}
-
-double integrand_g(double *xy, size_t dim, void *integrand_params_void) {
-  double cotpsi = 1.0/hypot(xy[0], xy[1]);
-  if (cotpsi >= cotpsi_max) return 0;
-  Integrand_params *integrand_params = (Integrand_params *)integrand_params_void;
-  double theta = gsl_spline_eval(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc);
-  double phi = atan2(xy[1], xy[0]);
-  double rgb[3] = {0};
-  double theta_phi[2] = {theta, phi};
-  integrand_params->system->angle_to_pixel(theta_phi, rgb);
-  return rgb[1]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc));
-}
-
-double integrand_b(double *xy, size_t dim, void *integrand_params_void) {
-  double cotpsi = 1.0/hypot(xy[0], xy[1]);
-  if (cotpsi >= cotpsi_max) return 0;
-  Integrand_params *integrand_params = (Integrand_params *)integrand_params_void;
-  double theta = gsl_spline_eval(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc);
-  double phi = atan2(xy[1], xy[0]);
-  double rgb[3] = {0};
-  double theta_phi[2] = {theta, phi};
-  integrand_params->system->angle_to_pixel(theta_phi, rgb);
-  return rgb[2]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc));
+  //printf("(x,y)=(%g,%g), tanpsi=%g, theta=%g, %g\n", xy[0], xy[1], tanpsi, theta, POW2(sin(theta))*(2+POW2(cotpsi)+POW2(tanpsi)));
+  return rgb[integrand_params->c]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, tanpsi, integrand_params->spline_data->acc))*POW2(sin(theta))*(2+POW2(tanpsi)+POW2(cotpsi));
 }
 
 int MC_pixel_render(const System *system, int i, int j, const Spline_data spline_data, double *rgb, int pixel_render_max, double pixel_render_err, gsl_rng *r){
@@ -61,37 +41,42 @@ int MC_pixel_render(const System *system, int i, int j, const Spline_data spline
   double y_lu = -j*dw + (system->H)*dw/2;
   double xl[2] = {x_lu, y_lu - dw};
   double xu[2] = {x_lu + dw, y_lu};
+  double err_dw = pixel_render_err*dw*dw;
+
+  //printf("integrating for (i,j)=(%d, %d)\n", i, j);
 
   //TODO: fill the MC intergral
-  // linear rgb = ∫ rgb(χ, ϕ) sin^2(χ) |χ'(ψ)| dx dy
-  // χ'(ψ) = χ'(cot(ψ)) * cot'(ψ) = - χ'(cot(ψ))/\sin^2(ψ)
+  // linear rgb = ∫ rgb(χ, ϕ) sin^2(χ)/sin^2(ψ) |χ'(ψ)| dx dy
+  // χ'(ψ) = χ'(tan(ψ)) * tan'(ψ) = - χ'(tan(ψ))/cos^2(ψ) = - χ'(tan(ψ))(1+tan^2(ψ))
   gsl_monte_miser_state *miser_state = gsl_monte_miser_alloc(2);
 
   double color;
   double err;
-  Integrand_params integrand_params = {system, &spline_data};
+  Integrand_params integrand_params = {0, system, &spline_data};
   gsl_monte_function F;
+  F.f = integrand;
   F.params = &integrand_params;
   F.dim = 2;
   //printf("i=%d, j=%d, rgb=( ", i, j);
   for (int c = 0; c < 3; c++) {
-    switch (c) {
-      case 0:
-        F.f = integrand_r;
-        break;
-      case 1:
-        F.f = integrand_g;
-        break;
-      case 2:
-        F.f = integrand_b;
-        break;
+    integrand_params.c = c;
+    if (i==0 && j==0) {
+      //printf("(x,y) = (%g, %g), intergrand[%d]=%g\n", xl[0], xl[1], c, integrand(xl, 2, &integrand_params));
     }
     gsl_monte_miser_init(miser_state);
     gsl_monte_miser_integrate(&F, xl, xu, 2, 100, r, miser_state, &color, &err);
     //printf("%g +- %g ", color, err);
-    rgb[c] = color;
+    if (err > err_dw) {
+      //int ncalls = 100*POW2(err/err_dw);
+      //printf("(i,j)=(%d,%d), err=%g, err_dw=%g, ncalls=%d\n", i, j, err, err_dw, ncalls);
+      //if (ncalls < 0 || ncalls > pixel_render_max) ncalls = pixel_render_max;
+      int ncalls = pixel_render_max;
+      //printf("(i,j)=(%d,%d), err=%g, err_dw=%g, ncalls=%d\n", i, j, err, err_dw, ncalls);
+      gsl_monte_miser_integrate(&F, xl, xu, 2, ncalls, r, miser_state, &color, &err);
     }
-  //printf(")\n");
+    rgb[c] = color/(dw*dw);
+  }
+  //printf("finished for (i,j)=(%d, %d)\n", i, j);
 
   return 0;
 }
@@ -134,5 +119,6 @@ int render(System *system, double *buffer, int pixel_render_max, double pixel_re
     gsl_interp_accel_free(acc);
     gsl_rng_free(r);
   }
+  printf("rendering finished.\n");
   return 0;
 }

src/write_png.c

@@ -14,22 +14,11 @@ static inline double linear_to_srgb(double x){
 }
 
 int buffer_normalize_srgb(double *buffer, int W, int H) {
-  double max = 0;
-  //for (int j = 0; j < H; j++) {
-  //  for (int i = 0; i < W; i++) {
-  //    for (int c = 0; c < 3; c++) {
-  //      if (max < buffer[color_index(i, j, c)]) max = buffer[color_index(i, j, c)];
-  //    }
-  //  }
-  //}
-  max = percentile(buffer, W*H*3, cutperc);
-  printf("max in the buffer: %g\n", max);
-  if (max == 0) return 1;
+  #pragma omp parallel for
   for (int j = 0; j < H; j++) {
     for (int i = 0; i < W; i++) {
       for (int c = 0; c < 3; c++) {
-        buffer[color_index(i, j, c)] /= max;
-        buffer[color_index(i, j, c)] = linear_to_srgb(buffer[color_index(i, j, c)]);
+        buffer[color_index(i, j, c)] = linear_to_srgb(SCALE*buffer[color_index(i, j, c)]);
       }
     }
   }
@@ -39,7 +28,7 @@ int buffer_normalize_srgb(double *buffer, int W, int H) {
 int write_png(char *filename, double *buffer, int W, int H) {
   buffer_normalize_srgb(buffer, W, H);
   size_t bufsize = W*H*3;
-  printf("bufsize: %d\n", bufsize);
+  //printf("bufsize: %d\n", bufsize);
   uint8_t *img = (uint8_t *)malloc(bufsize*sizeof(uint8_t));
   if (!img) {return 1;}