#ifndef __SOURCE_H__
#define __SOURCE_H__

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

typedef struct {
  vector coordinates;
  vector directions;
  vector surface_normal;
  rtd optical_path_length;
  rtd optical_path_difference;
  char v;
  int n_iteration;
  rtd throughput;
} ray;
struct bundle {
  int N_RAY;
  ray *d__ray;
  int N_BUNDLE;
  int N_RAY_TOTAL;
  vector *d__origin;
  double rot_angle;
  ray *d__chief_ray;
  vector *d__chief_origin;
  mask V;
  char geom[8];
  int N_RADIUS, N_THETA;
  int N_L;
  rtd L;
  rtd *d__sphere_distance, *d__sphere_radius;
  vector *d__sphere_origin;
  int *d__piston_mask;
  rtd refractive_index;
  double *d__Vx, *d__Vy;
  void setup(rtd RADIUS, int N_RADIUS, int N_THETA, vector origin, int N_SRC);
  void setup(rtd L, int N_L, vector origin, int N_SRC);
  void setup(rtd L, int N_L, vector origin, vector chief_origin, int N_SRC);
  void setup_free(int _N_RAY_, double *x, double *y, vector origin);
  void setup_free(double zenith, double azimuth,
                  int _N_RAY_, double *x, double *y, vector origin);
  void cleanup(void);
  void to_z_plane(rtd z_in);
  void to_focal_plane(rtd z_chief_on_axis, rtd rho_focal_plane);
  void to_sphere(vector sphere_origin);
  void to_sphere(rtd z_chief_on_axis, rtd rho_focal_plane);
  const void get_coordinates(double *d__coord);
  void get_chief_coordinates(double *d__coord);
  void get_sphere_origins(double *d__coord);
  void get_directions(double *d__dir);
  void get_chief_directions(double *d__dir);
  void get_chief_optical_path_length(double *d__opl);
  void get_optical_path_length(double *d__opl);
  void get_optical_path_difference(double *d__opd);
  void get_optical_path_difference(double *d__opd,
                                   float const delta_x, int N_x,
                                   float const delta_y, int N_y);
  void get_vignetting(double *d__v);
  void get_n_iteration(int *n_iteration);
  void gmt_truss_onaxis(const double scale);
  void gmt_truss_onaxis_1(const double scale);
  void gmt_truss_onaxis_2(const double scale);
  void gmt_truss_onaxis_3(const double scale);
  void gmt_pma_onaxis(void);
  void gmt_pma_plate_onaxis(void);
  void gmt_m2_baffle(const double baffle_radius, const double scale);
};

struct complex_amplitude {

  int N_PX, N;
  float *amplitude, *phase;
  mask *M;
  cublasHandle_t handle;
  float *buffer;

  void setup(int n_pixel);
  void setup(int n_pixel, int n_src);
  void cleanup(void);
  void reset(void);
  void reset(complex_amplitude *wavefront);
  void reset_amplitude(void);
  void reset_phase(void);
  void reset_phase(complex_amplitude *wavefront_prime);
  void add_phase(float alpha, float *phase_prime);
  void add_same_phase(float alpha, float *phase_prime);
  void masked(void);
  void masked(mask *M_in);
  void rms(float *rms);
  void finite_difference(float *sx, float *sy, int NL, float d);
  void finite_difference(float *sx, float *sy, int NL, float d, mask *valid_lenslet);
  void gradient_average(float *sx, float *sy, int NL, float d);
  void gradient_average(float *sx, float *sy, float d);
  void segments_gradient_average(float *sx, float *sy, float D, int *segment_markers);
  void segments_gradient_averageFast(float *sx, float *sy, float D, int *segment_markers);
  void show_phase(char *filename);
  void show_phase(char *filename, int N_SRC);
  void show_amplitude(char *filename);
  void show_amplitude(char *filename, int N, int M);
};

struct source {

  int N_SRC;
  float zenith, azimuth, height, theta_x, theta_y;
  rtd _zenith_64_, _azimuth_64_, _height_64_, _theta_x_64_, _theta_y_64_;
  const char *photometric_band;
  float magnitude, N_PHOTON;
  float fwhm;
  complex_amplitude wavefront;
  source *dev_ptr;
  char tag[8];
  char rays_exist;
  bundle rays;

  void setup(const char *_photometric_band,
             float zenith, float azimuth, float height);
  void setup(const char *_photometric_band,
             float zenith, float azimuth, float height,
             int resolution);
  void setup(const char *_photometric_band,
             float zenith, float azimuth, float height,
             const char *tag_in);
  void setup(const char *_photometric_band,
             float zenith, float azimuth, float height,
             int resolution, const char *tag_in);
  void setup(const char *_photometric_band,
             float *_zenith, float *_azimuth, float _height,
             int _N_SRC);
  void setup(const char *_photometric_band,
             float *_zenith, float *_azimuth, float _height,
             int _N_SRC, int resolution);
  void setup(const char *_photometric_band,
             float *_zenith, float *_azimuth, float _height,
             int _N_SRC, rtd _L_, int _N_L_, vector origin);
  void setup(const char *_photometric_band, float *magnitude,
             float *_zenith, float *_azimuth, float _height,
             int _N_SRC, rtd _L_, int _N_L_, vector origin);
  void setup(const char *_photometric_band, float *magnitude,
             rtd *_zenith, rtd *_azimuth, float _height,
             int _N_SRC, rtd _L_, int _N_L_, vector origin);
  void setup_chief(const char *_photometric_band, float *magnitude,
             rtd *_zenith, rtd *_azimuth, float _height,
             int _N_SRC, rtd _L_, int _N_L_,
             vector origin, vector chief_origin);
  void setup(const char *_photometric_band, float *_magnitude,
             rtd *_zenith, rtd *_azimuth, float _height, int _N_SRC,
             int _N_RAY_, double *x, double *y, vector origin);
   void cleanup(void);

  void reset_rays(void);
  void reset_rays(int RESET_RAYS_MASK);

  void opd2phase(void);
  void opd2phase(int RESET_RAYS_MASK);

  void info(void);

  void phase2file(const char *filename);

  float wavelength(void);
  float wavelength_micron(void);
  float spectral_bandwidth(void);
  float n_photon(void);
  float n_photon(float _magnitude_);
  float n_background_photon(float backgroundMagnitude);
  float wavenumber(void);
  void update_directions(double *zenith, double *azimuth, int N_DIR);
  void update_magnitude(float *magnitude, int N_MAG);
  void copy_magnitude(source *other_src);
  void optical_transfer_function(float2 *d__otf);
};

struct pssn {
  int N_O, N_O0, n_byte;
  float2 *d__O, *d__O0, *buffer;
  float2 *d__C;
  int N_PX, N;  
  float2 *d__W;
  int N_OTF, N_OTF2, NN;
  cufftHandle plan;
  cublasHandle_t handle;
  float num, *denom;
  void setup(source *src, const float r0, const float L0);
  void cleanup();
  void __otf__(source *src, float2 *_d__O_);
  void otf(source *src);
  void atm_otf(float d, float r0, float L0);
  float eval(void);
  float oeval(void);
  void eval(float *results);
  void oeval(float *results);
  void xotf(float *buffer);
  void O(float *buffer);
  void O0(float *buffer);
  void C(float *buffer);
  void B(float *buffer);
};

#endif // __SOURCE_H__