fixed albedo bug

c_src/_web.c

@@ -272,7 +272,7 @@ static PyObject *web_line_intersect(PyObject *self, PyObject *args)
 static PyObject *web_generate_planet(PyObject *self, PyObject *args)
 {
     int n_layers,n_1,n_2,bright_type,n_star;
-    double lambda0,phi0,p_u1,p_u2,rp;
+    double lambda0,phi0,p_u1,p_u2,rp,star_surface_bright,star_bright;
     PyObject *bright_obj,*teff_obj,*flux_obj;
 
     /* Parse the input tuple */
@@ -315,10 +315,13 @@ static PyObject *web_generate_planet(PyObject *self, PyObject *args)
     int n_temps=100;
     int n_bb_seg=20;
 
-    double star_bright = 1.0;
-
     double r2 = 1.0/rp; //invert planet radius ratio - planets always have radius 1 in this code
 
+
+    star_bright = 1.0;
+    star_surface_bright = star_bright/(M_PI*pow(r2,2));
+
+
     if(bright_type == 1 || bright_type == 3 || bright_type == 4 || bright_type == 8|| bright_type == 10){
         double l1 = brightness_params[1];
         double l2 = brightness_params[2];
@@ -329,16 +332,16 @@ static PyObject *web_generate_planet(PyObject *self, PyObject *args)
 //        star_bright = bb_flux(l1,l2,star_T,n_bb_seg);
         ypp = spline_cubic_set( n_star, stellar_teffs, stellar_fluxes, 0, 0, 0, 0 );
 
-        star_bright = spline_cubic_val( n_star, stellar_teffs, stellar_fluxes, ypp, star_T, &ypval, &yppval);
+        star_surface_bright = spline_cubic_val( n_star, stellar_teffs, stellar_fluxes, ypp, star_T, &ypval, &yppval);
         free(ypp);
 
-        star_bright = star_bright*M_PI*pow(r2,2);
+        star_bright = star_surface_bright*M_PI*pow(r2,2);
 
 
         bb_g = bb_grid(l1, l2, T_start, T_end,n_temps,n_bb_seg);
     }
 
-    map_model(planet_struct,n_layers,lambda0,phi0,p_u1,p_u2,bright_type,brightness_params,bb_g,star_bright);
+    map_model(planet_struct,n_layers,lambda0,phi0,p_u1,p_u2,bright_type,brightness_params,bb_g,star_surface_bright);
 
     /* Build the output tuple */
 
@@ -592,6 +595,8 @@ static PyObject *web_call_map_model(PyObject *self, PyObject *args)
 
     double star_bright = 1.0;
 
+    //NEED TO UPDATE THIS WITH CORRECT STAR BRIGHTNESS VALUES!//
+
     double *vals = call_map_model(la,lo,lambda0,phi0,bright_type,brightness_params,bb_g,0,0.0,0.0,0.0,0.0,star_bright);
 
     /* Build the output tuple */

c_src/brightness_maps.c

@@ -247,6 +247,7 @@ double lambertian(double lat, double lon, double insol, double albedo){
 
     if((-M_PI/2.0 <= lon) && (lon <= M_PI/2.0)){
         b = albedo*insol*cos(lat)*cos(lon)/M_PI;
+//        b = albedo*insol/M_PI;
     }
     else{
         b = 0.0;

c_src/web.c

@@ -14,6 +14,7 @@ double *lightcurve(int n_layers, int n_points, double *t, double tc, double per,
     double *coords;
     double p_blocked, p_bright,phase_z,phase_dz,phase_dt;
     double star_bright;
+    double star_surface_bright;
     double **bb_g;
     double *ypp;
 
@@ -37,6 +38,7 @@ double *lightcurve(int n_layers, int n_points, double *t, double tc, double per,
     double transit_z = transit_coords[3];
 
     star_bright = 1.0;
+    star_surface_bright = star_bright/(M_PI*pow(r2,2));
 
     // brightness model 1 is the Xi 2016 model, requires a stellar temperature
     if(brightness_model == 1 || brightness_model == 3 || brightness_model == 4 || brightness_model == 6 || brightness_model == 8|| brightness_model == 10){
@@ -49,11 +51,11 @@ double *lightcurve(int n_layers, int n_points, double *t, double tc, double per,
 //        star_bright = bb_flux(l1,l2,star_T,n_bb_seg);
 
         ypp = spline_cubic_set( nstars, stellar_teffs, stellar_fluxes, 0, 0, 0, 0 );
-        star_bright = spline_cubic_val( nstars, stellar_teffs, stellar_fluxes, ypp, star_T, &ypval, &yppval);
+        star_surface_bright = spline_cubic_val( nstars, stellar_teffs, stellar_fluxes, ypp, star_T, &ypval, &yppval);
         free(ypp);
 
 
-        star_bright = star_bright*M_PI*pow(r2,2);
+        star_bright = star_surface_bright*M_PI*pow(r2,2);
 
     // also requires the precomputation of the blackbody interpolation grid
 //        printf("%f %f %f %f %f %f\n",l1,l2,T_start,T_end,n_temps,n_bb_seg);
@@ -89,7 +91,7 @@ double *lightcurve(int n_layers, int n_points, double *t, double tc, double per,
 
 //        printf("bb_g 3 %f\n",bb_g[0][1]);
 
-        map_model(planet,n_layers,lambda0,phi0,u1,u2,brightness_model,brightness_params,bb_g,star_bright);
+        map_model(planet,n_layers,lambda0,phi0,u1,u2,brightness_model,brightness_params,bb_g,star_surface_bright);
 
         p_bright = 0.0;
         for (j = 0; j < pow(n_layers,2); j++) {