From de6d3bc04a691750e6ca68ac9c35de9203c9a895 Mon Sep 17 00:00:00 2001
From: Yingjie Wang <phywyj@gmail.com>
Date: Sun, 7 Dec 2025 22:47:53 -0500
Subject: [PATCH] update: use a single integrand function

---
 src/render.c | 49 +++++++------------------------------------------
 1 file changed, 7 insertions(+), 42 deletions(-)

diff --git a/src/render.c b/src/render.c
index a766411..fdddb7d 100644
--- a/src/render.c
+++ b/src/render.c
@@ -9,6 +9,7 @@
 #include "common.h"
 
 typedef struct {
+  int c;
   const System *system;
   const Spline_data *spline_data;
 } Integrand_params;
@@ -20,7 +21,7 @@ inline double xy_to_b(double x, double y) {
   return b;
 }
 
-double integrand_r(double *xy, size_t dim, void *integrand_params_void) {
+double integrand(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;
@@ -29,34 +30,7 @@ double integrand_r(double *xy, size_t dim, void *integrand_params_void) {
   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))*POW2(sin(theta))*POW2(1+POW2(cotpsi));
-  return rgb[0];
-}
-
-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))*POW2(sin(theta))*POW2(1+POW2(cotpsi));
-  return rgb[1];
-}
-
-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))*POW2(sin(theta))*POW2(1+POW2(cotpsi));
-  return rgb[2];
+  return rgb[integrand_params->c]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc))*POW2(sin(theta))*POW2(1+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){
@@ -73,23 +47,14 @@ int MC_pixel_render(const System *system, int i, int j, const Spline_data spline
 
   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;
     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);
@@ -132,7 +97,7 @@ int render(System *system, double *buffer, int pixel_render_max, double pixel_re
     gsl_interp_accel *acc = gsl_interp_accel_alloc();
     Spline_data spline_data = {spline, acc};
     gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
-    
+ 
     #pragma omp for
     for(int j = 0; j < H; j++) {
       for (int i = 0; i < W; i++) {