From 480c49dcf09fa9bee84dd978bede453dc1a8fa81 Mon Sep 17 00:00:00 2001
From: Yingjie Wang <phywyj@gmail.com>
Date: Sun, 7 Dec 2025 22:41:59 -0500
Subject: [PATCH] fix: fix intergrad

---
 src/common.h |  2 ++
 src/render.c | 14 +++++++++-----
 2 files changed, 11 insertions(+), 5 deletions(-)

diff --git a/src/common.h b/src/common.h
index 55feacb..f6ed152 100644
--- a/src/common.h
+++ b/src/common.h
@@ -5,6 +5,8 @@
 #include <gsl/gsl_errno.h>
 #include <math.h>
 
+#define POW2(x) (x*x)
+
 #define FLAT 0
 #define PI 3.1415926535897932384626433832795028841971693993751058
 #define Rs 1
diff --git a/src/render.c b/src/render.c
index 3f6d85b..a766411 100644
--- a/src/render.c
+++ b/src/render.c
@@ -6,6 +6,7 @@
 #include <math.h>
 #include <stdio.h>
 #include "render.h"
+#include "common.h"
 
 typedef struct {
   const System *system;
@@ -28,7 +29,8 @@ 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));
+  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) {
@@ -40,7 +42,8 @@ double integrand_g(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[1]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc));
+  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) {
@@ -52,7 +55,8 @@ double integrand_b(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[2]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc));
+  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];
 }
 
 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){
@@ -63,7 +67,7 @@ int MC_pixel_render(const System *system, int i, int j, const Spline_data spline
   double xu[2] = {x_lu + dw, y_lu};
 
   //TODO: fill the MC intergral
-  // linear rgb = ∫ rgb(χ, ϕ) sin^2(χ) |χ'(ψ)| dx dy
+  // linear rgb = ∫ rgb(χ, ϕ) sin^2(χ)/sin^2(ψ) |χ'(ψ)| dx dy
   // χ'(ψ) = χ'(cot(ψ)) * cot'(ψ) = - χ'(cot(ψ))/\sin^2(ψ)
   gsl_monte_miser_state *miser_state = gsl_monte_miser_alloc(2);
 
@@ -94,7 +98,7 @@ int MC_pixel_render(const System *system, int i, int j, const Spline_data spline
       if (ncalls > pixel_render_max) ncalls = pixel_render_max;
       gsl_monte_miser_integrate(&F, xl, xu, 2, ncalls, r, miser_state, &color, &err);
     }
-    rgb[c] = color;
+    rgb[c] = color/(dw*dw);
   }
   //printf(")\n");