From b9954ff8fad8cc4080fe7261ee0d305b9b0e5687 Mon Sep 17 00:00:00 2001
From: Yingjie Wang <phywyj@gmail.com>
Date: Mon, 10 Nov 2025 15:14:21 -0500
Subject: [PATCH] fix bug and cleanup

---
 src/init.c | 32 ++++++++++++++++++--------------
 1 file changed, 18 insertions(+), 14 deletions(-)

diff --git a/src/init.c b/src/init.c
index d68a381..e7c51a7 100644
--- a/src/init.c
+++ b/src/init.c
@@ -1,8 +1,11 @@
 #include <gsl/gsl_matrix_double.h>
 #include <gsl/gsl_odeiv2.h>
 #include <gsl/gsl_matrix.h>
+#include <stdio.h>
 #include "common.h"
 
+#define MAX_DEPTH 15
+
 typedef struct {
   int capacity;
   int size;
@@ -24,11 +27,11 @@ static int jac(double phi, const double y[], double *dfdy, double dfdphi[], void
   gsl_matrix *m = &dfdy_mat.matrix;
   gsl_matrix_set(m, 0, 0, 0.0);
   gsl_matrix_set(m, 0, 1, 1.0);
-  gsl_matrix_set(m, 1, 0, -3*M*y[0]-1.0);
+  gsl_matrix_set(m, 1, 0, 63*M*y[0]-1.0);
   gsl_matrix_set(m, 1, 1, 0.0);
 
   dfdphi[0] = 0.0;
-  dfdphi[0] = 0.0;
+  dfdphi[1] = 0.0;
 
   return GSL_SUCCESS;
 }
@@ -61,8 +64,7 @@ static inline double find_root_5ord(double u1, double u2, double v1, double v2,
 
 // \chi(b) is the totol deflection angle with parameter b
 double chi(double b){
-  double x;
-  gsl_odeiv2_system sys = {func, jac, 2, &x};
+  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);
   
   double u0 = 1.0/R0;
@@ -85,7 +87,7 @@ double chi(double b){
       printf("Error: GSL driver returns %d\n", status);
       break;
     }
-    printf("phi: %g, y[0]: %g, y[1]: %g\n",phi, y[0], y[1]);
+    //printf("phi: %g, y[0]: %g, y[1]: %g\n",phi, y[0], y[1]);
   }
 
   //printf("%g, %g \n",find_root(y_previous, phi_previous),find_root(y, phi));
@@ -106,25 +108,26 @@ void sample_push(SampleData *sample, double x, double y){
     sample->x = realloc(sample->x, sizeof(double) * newcap);
     sample->y = realloc(sample->y, sizeof(double) * newcap);
     sample->capacity = newcap;
+    //printf("enlarge sample capacity: %d\n", newcap);
   }
   sample->x[sample->size] = x;
   sample->y[sample->size] = y;
   sample->size += 1;
 }
 
-static void refine_interval(double xl, double xr, double yl, double yr, double rela_err_limit, SampleData *sample){
-  printf("calculating interval (%g, %g)\n", xl, xr);
+static void refine_interval(double xl, double xr, double yl, double yr, double rela_err_limit, SampleData *sample, int depth){
+  //printf("calculating interval (%.16g, %.16g)\n", xl, xr);
   double xm = (xl + xr)/2;
   double ym = chi(xm);
-  int need_refine = (fabs(ym-(yl+yr)/2) >= fabs(rela_err_limit * ym));
+  int need_refine = (fabs(ym-(yl+yr)/2) >= fabs(rela_err_limit * ym)) && depth <= MAX_DEPTH;
   if (need_refine) {
-    printf(" ");
-    refine_interval(xl, xm, yl, ym, rela_err_limit, sample);
+    //printf(" ");
+    refine_interval(xl, xm, yl, ym, rela_err_limit, sample, depth + 1);
   }
   sample_push(sample, xm, ym);
   if (need_refine) {
-    printf(" ");
-    refine_interval(xm, xr, ym, yr, rela_err_limit, sample);
+    //printf(" ");
+    refine_interval(xm, xr, ym, yr, rela_err_limit, sample, depth + 1);
   }
   return ;
 }
@@ -139,14 +142,15 @@ int init(double bmax, double rela_err_limit){
   sample.y[0] = chi(bmin);
   double chi_bmax = chi(bmax);
 
-  refine_interval(bmin, bmax, sample.y[0], chi_bmax, rela_err_limit, &sample);
+  refine_interval(bmin, bmax, sample.y[0], chi_bmax, rela_err_limit, &sample, 1);
 
   sample_push(&sample, bmax, chi_bmax);
 
   //sort
 
   //output
-  printf("Total sample points number: %d\n", sample.size);
+  //printf("Total sample points capacity: %d\n", sample.capacity);
+  //printf("Total sample points number: %d\n", sample.size);
   for (int i = 0; i < sample.size; i++) {
     printf("%.16g %.16g\n", sample.x[i], sample.y[i]);
   }