#ifndef __RAYTRACING_H__
#define __RAYTRACING_H__

#ifndef __UTILITIES_H__
#include "utilities.h"
#endif

#ifndef __SOURCE_H__
#include "source.h"
#endif

#define SNELL_N_ITERATION 100
#define TOL 1E-9

struct coordinate_system{
  vector *origin;
  vector *euler_angles;
  int N;
  rtd *R, *d__R;
  float *float_R;
  vector *d__origin;
  char tag[32];
  void setup(void);
  void setup(vector _origin_, vector _euler_angles_);
  void setup(vector *_origin_, vector *_euler_angles_, int _N_);
  void setup(vector *_origin_, vector *_euler_angles_, int _N_, char *_tag_);
  void cleanup(void);
  void info(void);
  void info_details(void);
  void update(void);
  void update(vector _origin_, vector _euler_angles_,int idx);
};

struct conic {

  coordinate_system ref_frame;
  vector origin, *d__origin;
  rtd c, k, refractive_index;
  int even_asphere_N;
  rtd *d__even_asphere_a;

  void setup(rtd _c_, rtd _k_);
  void setup(rtd _c_, rtd _k_, vector _origin_, vector _euler_angles_);
  void setup(rtd _c_, rtd _k_, vector _origin_, vector _euler_angles_, 
             vector conic_origin);
  void setup(rtd _c_, rtd _k_, vector _origin_, vector _euler_angles_, 
             vector conic_origin, rtd _refractive_index_);
  void setup(rtd _c_, rtd _k_, vector _origin_, vector _euler_angles_, 
             vector conic_origin, rtd _refractive_index_,
             int asphere_N, rtd *asphere_a);
  void cleanup(void);
  void trace(bundle *src);
  void info(void);
};

__host__ __device__ rtd even_asphere(rtd r, int N, rtd *a);
__host__ __device__ rtd partial_x_or_y_even_asphere(rtd x_or_y, rtd r2, int N, rtd *a);

struct zernikeS {
  int max_n;
  unsigned int j, n, m, n_mode;
  rtd *a, *d__a;
  rtd *bx, *by, *d__bx, *d__by;
  unsigned int *bx_row_idx, *bx_col_ptr,
               *by_row_idx, *by_col_ptr,
               bx_nnz, by_nnz;
  rtd *cx, *d__cx, *cy, *d__cy;
  int N;
  coordinate_system ref_frame;
  void setup(int max_n, rtd *a);
  void setup(int _max_n_, rtd *_a_, int _N_);
  void setup(int max_n, rtd *a, vector _origin_, vector _euler_angles_);
  void setup(int max_n, rtd *a, vector _origin_, vector _euler_angles_, int _N_);
  void cleanup(void);
  void surface(rtd *S, rtd *r, rtd *o, const int N);
  void surface(rtd *S, rtd *r, rtd *o, const int N, int surf_id);
  void update(rtd *a);
  void surface_derivative_coefs(void);
  void surface_derivatives(rtd *d__dSdx, rtd *d__dSdy, rtd *d__r, rtd *d__o, const int N);
  void surface_and_derivatives(rtd *d__S, rtd *d__dSdx, rtd *d__dSdy,
                               rtd *d__r, rtd *d__o, const int N);
  void surface_and_derivatives(rtd *d__S, rtd *d__dSdx, rtd *d__dSdy,
                               rtd *d__r, rtd *d__o, const int N, int surf_id);
  void projection(float *d__phase, rtd *d__r, rtd *d__o, const int N);
};

__host__ __device__  rtd zernike_surface(rtd r, rtd o, int max_n, rtd *a);
__host__ __device__  rtd zernike_surface_x_derivative(rtd r, rtd o, int max_n, rtd *a);
__host__ __device__  rtd zernike_surface_y_derivative(rtd r, rtd o, int max_n, rtd *a);
__host__ __device__ void zernike_surface_and_partials(rtd *ZS, rtd *dZSdx, rtd *dZSdy,
                                                      rtd r, rtd o, int max_n,
                                                      rtd *a, rtd *ax, rtd *ay);
int zern_dx_coef(rtd *b, int j, int n, int m, int jp, int np, int mp);
int zern_dy_coef(rtd *b,int j, int n, int m, int jp, int np, int mp);

__host__ __device__ rtd conic_equation(vector *v, vector *v0, 
                                       const rtd k, const rtd c,
                                       const int N, rtd *a);
__device__ rtd  conic_surface(vector *v, vector *v0, 
                              const rtd k, const rtd c,
                              const int N, rtd *a);
__device__ rtd partial_x_conic_surface(vector *v, vector *v0,
                                       const rtd k, const rtd c,
                                       const int N, rtd *a);
__device__ rtd partial_y_conic_surface(vector *v, vector *v0,
                                       const rtd k, const rtd c,
                                       const int N, rtd *a);
__device__ rtd partial_z_conic_surface(void);

struct aperture {
  int D_px;
  rtd D, ri, D_seg;
  int N;
  mask V;
  coordinate_system ref_frame;

  void setup(rtd _D_, int _D_px_,
             vector _origin_, vector _euler_angles_);
  void setup(rtd _D_, rtd _ri_, int _D_px_,
             vector _origin_, vector _euler_angles_);
  void setup(rtd _D_, rtd _ri_, rtd _D_seg_, int _D_px_,
             vector *_origin_, vector *_euler_angles_, int N_segment);
  void setup_GMT_M1(rtd _D_, int _D_px_);
  void cleanup(void);
  void vignetting(bundle *rays);

};
__global__ void intersection(char *mask, ray *d__ray, int N_RAY);

__host__ __device__  void forward_transform(vector *v_out, vector *v_in,
                                   rtd *d__R, vector *d__origin);
__host__ __device__  void forward_transform_centered(vector *v_out, vector *v_in,
                                            rtd *d__R);
__host__ __device__  void backward_transform(vector *v_out, vector *v_in,
                                    rtd *d__R, vector *d__origin);
__host__ __device__  void backward_transform_centered(vector *v_out, vector *v_in,
                                             rtd *d__R);
__device__ rtd aspheric_surface(vector *v, vector *v0,
                                const rtd k, const rtd c,
                                const int max_n, rtd *a, const rtd R);
__device__ rtd zernike_surface(vector *v, vector *v0,
                               const int max_n, rtd *a, const rtd R);
__device__ rtd partial_x_aspheric_surface(vector *v, vector *v0,
                                          const rtd k, const rtd c,
                                          const int max_n, rtd *cx,
                                          const rtd R);
__device__ rtd partial_x_zernike_surface(vector *v, vector *v0,
                                         const int max_n, rtd *cx,
                                         const rtd R);
__device__ rtd partial_y_aspheric_surface(vector *v, vector *v0,
                                          const rtd k, const rtd c,
                                          const int max_n, rtd *cy,
                                          const rtd R);
__device__ rtd partial_y_zernike_surface(vector *v, vector *v0,
                                         const int max_n, rtd *cy,
                                         const rtd R);
__device__ rtd partial_z_aspheric_surface(void);
__device__ rtd partial_z_zernike_surface(void);


void transform_to_S(bundle *rays, conic *F);
void transform_to_S(bundle *rays, aperture *A);
void transform_to_R(bundle *rays, conic *F);
void transform_to_R(bundle *rays, aperture *A);
void intersect(bundle *rays, conic *F);
void reflect(bundle *rays);
void refract(bundle *rays, rtd mu);
#endif // __RAYTRACING_H__