before bug

c_src/_web.c

@@ -33,6 +33,8 @@ static PyObject *web_blocked(PyObject *self, PyObject *args);
 static PyObject *web_zhang_2016(PyObject *self, PyObject *args);
 static PyObject *web_separation_of_centers(PyObject *self, PyObject *args);
 static PyObject *web_lightcurve(PyObject *self, PyObject *args);
+static PyObject *web_calc_phase(PyObject *self, PyObject *args);
+static PyObject *web_calc_substellar(PyObject *self, PyObject *args);
 
 static PyMethodDef module_methods[] = {
     {"heron", web_heron, METH_VARARGS, web_docstring},
@@ -47,6 +49,8 @@ static PyMethodDef module_methods[] = {
     {"zhang_2016", web_zhang_2016, METH_VARARGS, quad_docstring},
     {"separation_of_centers", web_separation_of_centers, METH_VARARGS, quad_docstring},
     {"lightcurve", web_lightcurve, METH_VARARGS, quad_docstring},
+    {"calc_phase", web_calc_phase, METH_VARARGS, quad_docstring},
+    {"calc_substellar", web_calc_substellar, METH_VARARGS, quad_docstring},
     {NULL, NULL, 0, NULL}
 };
 
@@ -280,6 +284,8 @@ static PyObject *web_generate_planet(PyObject *self, PyObject *args)
     /* Get pointers to the data as C-types. */
     double *brightness_params    = (double*)PyArray_DATA(bright_array);
 
+    printf("%f %f\n",lambda0,phi0);
+
     map_model(planet_struct,n_layers,lambda0,phi0,p_u1,p_u2,bright_type,brightness_params);
 
     /* Build the output tuple */
@@ -368,6 +374,52 @@ static PyObject *web_separation_of_centers(PyObject *self, PyObject *args)
     return ret;
 }
 
+static PyObject *web_calc_phase(PyObject *self, PyObject *args)
+{
+    double t,tc,per;
+
+    /* Parse the input tuple */
+    if (!PyArg_ParseTuple(args, "ddd", &t,&tc,&per))
+        return NULL;
+
+    /* Call the external C function to compute the area. */
+    double output = calc_phase(t,tc,per);
+
+    PyObject *ret = Py_BuildValue("d",output);
+
+    return ret;
+}
+
+static PyObject *web_calc_substellar(PyObject *self, PyObject *args)
+{
+    double phase;
+    PyObject *c_obj;
+
+    /* Parse the input tuple */
+    if (!PyArg_ParseTuple(args, "dO", &phase,&c_obj))
+        return NULL;
+
+    PyObject *c_array = PyArray_FROM_OTF(c_obj, NPY_DOUBLE, NPY_IN_ARRAY);
+
+    /* If that didn't work, throw an exception. */
+    if (c_array == NULL) {
+        Py_XDECREF(c_array);
+        return NULL;
+    }
+
+    double *c    = (double*)PyArray_DATA(c_array);
+
+    /* Call the external C function to compute the area. */
+    double *output = calc_substellar(phase,c);
+
+    PyObject *ret = Py_BuildValue("[d,d]",output[0],output[1]);
+
+    Py_DECREF(c_array);
+    free(c);
+
+    return ret;
+}
+
 static PyObject *web_lightcurve(PyObject *self, PyObject *args)
 {
     int n_layers, bright_type;

c_src/old_web.c

@@ -0,0 +1,123 @@
+#include "generate.h"
+#include "blocked.h"
+#include "ephemeris.h"
+#include "math.h"
+#include "blackbody.h"
+#include <stdlib.h>
+#include <stdio.h>
+
+double *lightcurve(int n_layers, int n_points, double *t, double tc, double per, double a, double inc, double ecc, double omega, double a_rs, double rp,double u1, double u2,int brightness_model,double *brightness_params){
+    int n,j;
+    double phase,lambda0,phi0;
+    double *coords;
+    double p_blocked, p_bright,phase_z,phase_dz,phase_dt;
+    double star_bright;
+
+    double r2 = 1.0/rp; //invert planet radius ratio - planets always have radius 1 in this code
+
+    double c = 299792458.0;
+
+    double l1 = 1.1e-6;
+    double l2 = 1.7e-6;
+    int n_bb_seg = 10;
+
+    double *output = malloc(sizeof(double) * n_points);
+
+    // generate the planet grid
+    double **planet = generate_planet(n_layers);
+
+    double *transit_coords = separation_of_centers(tc,tc,per,a,inc,ecc,omega,a_rs,r2);
+
+    double transit_z = transit_coords[3];
+
+    if(brightness_model == 0){
+        star_bright = 1.0;
+    }
+
+    // brightness model 1 is the Xi 2016 model, requires a stellar temperature
+    if(brightness_model == 1 || brightness_model == 3 || brightness_model == 4){
+        double star_T =brightness_params[0];
+        star_bright = bb_flux(l1,l2,star_T,n_bb_seg);
+        star_bright = star_bright*M_PI*pow(r2,2);
+    }
+
+    free(coords);
+
+    for (n = 0; n < n_points; n++) {
+
+        double *old_coords = separation_of_centers(t[n],tc,per,a,inc,ecc,omega,a_rs,r2);
+        phase = ((t[n]-tc)/per);
+
+        // make correction for finite light travel speed
+
+        phase_z = old_coords[3];
+        phase_dz = transit_z-phase_z;
+        phase_dt = (phase_dz/c)/(3600.0*24.0);
+
+        free(old_coords);
+        double *coords = separation_of_centers(t[n]-phase_dt,tc,per,a,inc,ecc,omega,a_rs,r2);
+
+        if(phase > 1){
+            phase = phase - floor(phase);
+        }
+        if(phase < 0){
+            phase = phase + ceil(phase) + 1;
+        }
+
+        lambda0 = (M_PI+(phase*2*M_PI));
+        if(lambda0 > 2*M_PI){
+            lambda0 = lambda0 - 2*M_PI;
+        }
+        if(lambda0 < -2*M_PI){
+            lambda0 = lambda0 + 2*M_PI;
+        }
+
+        phi0 = atan2(coords[1],coords[2]);
+
+        map_model(planet,n_layers,lambda0,phi0,u1,u2,brightness_model,brightness_params);
+
+        p_bright = 0.0;
+        for (j = 0; j < pow(n_layers,2); j++) {
+            p_bright = p_bright + planet[j][16]*planet[j][15];
+        }
+
+        if(coords[2] < 0){
+            p_blocked = blocked(planet,n_layers,coords[0],coords[1],r2);
+        }
+        else{
+            // PRIMARY TRANSIT SHOULD GO IN HERE!
+            p_blocked = 0.0;
+        }
+        output[n] = (star_bright + p_bright - p_blocked)/star_bright;
+
+        free(coords);
+    }
+
+    int n_segments = pow(n_layers,2);
+    for (int i = 0; i < n_segments; ++i) {
+      free(planet[i]);
+      // each i-th pointer is now pointing to dynamic array (size 10) of actual int values
+    }
+    free(planet);
+    free(coords);
+    free(transit_coords);
+
+    return output;
+}
+
+double *call_blocked(int n_layers, int n_points, double *x2, double *y2, double r2){
+    int n;
+
+    // generate the planet grid
+    double **planet = generate_planet(n_layers);
+
+    double *output = malloc(sizeof(double) * n_points);
+
+    for (n = 0; n < n_points; n++) {
+        output[n] = blocked(planet,n_layers,x2[n],y2[n],r2);
+    }
+
+    free(planet);
+
+    return output;
+}

c_src/web.c

@@ -3,6 +3,7 @@
 #include "ephemeris.h"
 #include "math.h"
 #include "blackbody.h"
+#include "web.h"
 #include <stdlib.h>
 #include <stdio.h>
 
@@ -46,33 +47,25 @@ double *lightcurve(int n_layers, int n_points, double *t, double tc, double per,
     for (n = 0; n < n_points; n++) {
 
         double *old_coords = separation_of_centers(t[n],tc,per,a,inc,ecc,omega,a_rs,r2);
-        phase = ((t[n]-tc)/per);
+
+        phase = calc_phase(t[n],tc,per);
+
 
         // make correction for finite light travel speed
 
         phase_z = old_coords[3];
         phase_dz = transit_z-phase_z;
         phase_dt = (phase_dz/c)/(3600.0*24.0);
 
-        free(old_coords);
-        double *coords = separation_of_centers(t[n]-phase_dt,tc,per,a,inc,ecc,omega,a_rs,r2);
+        double *substellar = calc_substellar(phase,coords);
 
-        if(phase > 1){
-            phase = phase - floor(phase);
-        }
-        if(phase < 0){
-            phase = phase + ceil(phase) + 1;
-        }
+        lambda0 = substellar[0];
+        phi0 = substellar[1];
 
-        lambda0 = (M_PI+(phase*2*M_PI));
-        if(lambda0 > 2*M_PI){
-            lambda0 = lambda0 - 2*M_PI;
-        }
-        if(lambda0 < -2*M_PI){
-            lambda0 = lambda0 + 2*M_PI;
-        }
+        free(old_coords);
+        free(substellar);
 
-        phi0 = atan2(coords[1],coords[2]);
+        double *coords = separation_of_centers(t[n]-phase_dt,tc,per,a,inc,ecc,omega,a_rs,r2);
 
         map_model(planet,n_layers,lambda0,phi0,u1,u2,brightness_model,brightness_params);
 
@@ -105,6 +98,41 @@ double *lightcurve(int n_layers, int n_points, double *t, double tc, double per,
     return output;
 }
 
+double calc_phase(double t, double t0, double per){
+    double phase;
+
+    phase = ((t-t0)/per);
+
+    if(phase > 1){
+        phase = phase - floor(phase);
+    }
+    if(phase < 0){
+        phase = phase + ceil(phase) + 1;
+    }
+
+    return phase;
+}
+
+double *calc_substellar(double phase, double *coords){
+    double lambda0,phi0;
+    double *output = malloc(sizeof(double) * 2);
+
+    lambda0 = (M_PI+(phase*2*M_PI));
+    if(lambda0 > 2*M_PI){
+        lambda0 = lambda0 - 2*M_PI;
+    }
+    if(lambda0 < -2*M_PI){
+        lambda0 = lambda0 + 2*M_PI;
+    }
+
+    phi0 = atan2(coords[1],coords[2]);
+
+    output[0] = lambda0;
+    output[1] = phi0;
+    return output;
+
+}
+
 double *call_blocked(int n_layers, int n_points, double *x2, double *y2, double r2){
     int n;
 

c_src/web.h

@@ -1,2 +1,4 @@
 double *lightcurve(int n_layers, int n_points, double *t, double tc, double per, double a, double inc, double ecc, double omega, double r_s, double r2,double u1, double u2,int brightness_model,double* brightness_params);
-double *call_blocked(int n_layers, int n_points, double *x2, double *y2, double r2);
+double *call_blocked(int n_layers, int n_points, double *x2, double *y2, double r2);
+double calc_phase(double t, double t0, double per);
+double *calc_substellar(double phase, double *coords);