diff --git a/src/server/backends/SIM/rt_sim.c b/src/server/backends/SIM/rt_sim.c
index d6d95374cd4bd02993889882e092c5dcd0d2e4f1..7ebe31d662e6cbbd91efcfd306bd225210eae6f1 100644
--- a/src/server/backends/SIM/rt_sim.c
+++ b/src/server/backends/SIM/rt_sim.c
@@ -66,14 +66,21 @@
#define SIM_V_RED_KMS 400.0
#define SIM_V_BLU_KMS -400.0
#define SIM_V_RES_KMS 1.0
+#define SIM_HI_BINS 801
+#define SIM_HI_AMP_CAL 10.0 /* conversion from cK to mK */
+
+#define SIM_V_REF_HZ (SIM_V_REF_MHZ * 1e6)
/* default allowed HW ranges */
-#define SIM_FREQ_MIN_HZ DOPPLER_FREQ(SIM_V_RED_KMS, SIM_V_REF_MHZ * 1e6)
-#define SIM_FREQ_MAX_HZ DOPPLER_FREQ(SIM_V_BLU_KMS, SIM_V_REF_MHZ * 1e6)
-#define SIM_FREQ_STP_HZ DOPPLER_FREQ_REL(SIM_V_RES_KMS, SIM_V_REF_MHZ * 1e6)
+#define SIM_FREQ_MIN_HZ DOPPLER_FREQ(SIM_V_RED_KMS, SIM_V_REF_HZ)
+#define SIM_FREQ_MAX_HZ DOPPLER_FREQ(SIM_V_BLU_KMS, SIM_V_REF_HZ)
+#define SIM_FREQ_STP_HZ DOPPLER_FREQ_REL(SIM_V_RES_KMS, SIM_V_REF_HZ)
+
#define SIM_IF_BW_HZ (SIM_FREQ_MAX_HZ - SIM_FREQ_MIN_HZ)
-#define SIM_BINS ((SIM_V_RED_KMS - SIM_V_BLU_KMS) / SIM_V_RES_KMS)
+//#define SIM_BINS (SIM_IF_BW_HZ / SIM_V_RES_KMS)
+#define SIM_BINS ((SIM_V_RED_KMS - SIM_V_BLU_KMS) / SIM_V_RES_KMS)
+
#define SIM_BW_DIV_MAX 0
/* other properties */
@@ -82,6 +89,9 @@
#define SIM_TSYS 100. /* default system temperature */
#define SIM_EFF 0.6 /* default efficiency */
+#define SIM_SIG_NOISE 12. /* default noise */
+#define SIM_READOUT_HZ 1 /* default readout frequency */
+#define SIM_SUN_SFU 48. /* Sun @1415 as of Jun 11 2019 12 UTC */
#define SKY_GAUSS_INTG_STP 0.10 /* integration step for gaussian */
#define SKY_BASE_RES 0.5 /* resolution of the underlying data */
@@ -139,6 +149,9 @@ static struct {
gdouble r_beam; /* the beam radius */
gdouble tsys; /* system temperature */
gdouble eff; /* system efficiency (eta) */
+ gdouble sig_n; /* noise sigma */
+ gdouble readout_hz; /* spectrum readout frequency */
+ gdouble sun_sfu; /* sun solar flux units */
struct {
GdkPixbuf *pb_sky;
@@ -177,6 +190,9 @@ static struct {
.r_beam = SIM_BEAM_RADIUS,
.tsys = SIM_TSYS,
.eff = SIM_EFF,
+ .sig_n = SIM_SIG_NOISE,
+ .readout_hz = SIM_READOUT_HZ,
+ .sun_sfu = SIM_SUN_SFU,
};
/**
@@ -661,32 +677,37 @@ static gboolean sim_sky_draw(GtkWidget *w, cairo_t *cr, gpointer data)
-#define UNOISE (rand()/((double)RAND_MAX + 1.))
+
+
/**
- * @brief gaussian noise macro
+ * @brief get the half width of the gaussian scaled by a beam radius
+ * (== 1/2 FWHM) at given sigma
*/
-#define GNOISE (UNOISE + UNOISE + UNOISE + UNOISE + UNOISE + UNOISE + UNOISE + UNOISE + UNOISE + UNOISE + UNOISE + UNOISE - 6.)
+static double gauss_half_width_sigma_r(double sigma, double r)
+{
+ return r * sigma * sqrt(2.0 * log(2.0));
+}
/**
- * @brief returns the value of the width-normalised gaussian for a given FWHM
+ * @brief returns the value of the width-normalised gaussian for a given
+ * beam radius (== 1/2 FWHM)
*/
-static double gauss_norm(double x, double fwhm)
+static double gauss_norm(double sigma, double r)
{
/* The FWHM auf a gaussian with sigma == 1 is 2*sqrt(2*ln(2)) ~2.35482
* To normalise x, we multiply by the latter and divide by two, since
* we moved the factor from the divisor to the dividend.
* We need only half of the FWHM, i.e. 1 / (fwhm / 2) which == 2 / fwhm,
- * so our final scale factor is 2.35482/fwhm
+ * so our final scale factor is 2.35482/fwhm where fwhm == r
*/
+ sigma = sigma / r * 2.0 * sqrt(2.0 * log(2.0));
- x = x / fwhm * 2.*sqrt(2.*log(2.));
-
- return 1.0 / sqrt(2.0 * M_PI) * exp( -0.5 * x * x);
+ return 1.0 / sqrt(2.0 * M_PI) * exp(-0.5 * sigma * sigma);
}
/**
@@ -730,208 +751,900 @@ static gpointer sim_spec_extract_HI_survey(gdouble glat, gdouble glon)
if (!data)
data = fopen("../data/sky_vel.dat", "rb");
- if (!data) {
- GtkWidget *dia;
- GtkWindow * win;
+ if (!data) {
+ GtkWidget *dia;
+ GtkWindow * win;
+
+ dia = gtk_message_dialog_new(NULL,
+ GTK_DIALOG_MODAL,
+ GTK_MESSAGE_ERROR,
+ GTK_BUTTONS_CLOSE,
+ "Please place sky_vel.dat in the "
+ "directory path this program is executed in.");
+
+ gtk_dialog_run(GTK_DIALOG(dia));
+ gtk_widget_destroy(dia);
+
+
+ exit(0);
+ /** XXX **/
+ g_error("%s: error opening file", __func__);
+ return 0;
+ }
+
+
+
+
+ map = g_malloc(VEL * SKY_WIDTH * SKY_HEIGHT * sizeof(guint16));
+
+ fread(map, sizeof(guint16), VEL * SKY_WIDTH * SKY_HEIGHT, data);
+ fclose(data);
+ }
+
+ as = g_malloc(VEL * sizeof(guint16));
+ offset = get_offset(glat, glon);
+
+ memcpy(as, &map[offset], sizeof(guint16) * VEL);
+
+ return as;
+}
+
+
+/**
+ * @brief order two frequencies so f0 < f1
+ *
+ * @param f0 the first frequency
+ * @param f1 the second frequency
+ */
+
+static void sim_order_freq(gdouble *f0, gdouble *f1)
+{
+ double tmp;
+
+ if (!f0)
+ return;
+
+ if (!f1)
+ return;
+
+ if ((*f0) < (*f1))
+ return;
+
+ tmp = (*f1);
+ (*f1) = (*f0);
+ (*f0) = tmp;
+}
+
+
+/**
+ * @brief clamp a frequency to the simulator settings
+ *
+ * @param f the frequency in Hz to clamp into range
+ *
+ * @returns the clamped frequency in Hz
+ */
+
+static gdouble sim_clamp_freq_range(gdouble f)
+{
+ if (f < sim.radio.freq_min_hz)
+ return sim.radio.freq_min_hz;
+
+ if (f > sim.radio.freq_max_hz)
+ return sim.radio.freq_max_hz;
+
+ return f;
+}
+
+
+/**
+ * @brief round off a frequency to the simulator spectral resolution setting
+ *
+ * @param f the frequency in Hz to round off
+ *
+ * @returns the rounded off frequency in Hz
+ */
+
+static gdouble sim_round_freq_to_res(gdouble f)
+{
+ return round(f / sim.radio.freq_inc_hz) * sim.radio.freq_inc_hz;
+}
+
+
+/**
+ * @brief get the number of data bins for a frequency range given the
+ * simulator configuration
+ *
+ * @param f0 the lower bound frequency in Hz
+ * @param f1 the upper bound frequency in Hz
+ *
+ * @returns the number of bins needed
+ * @note f0, f1 must be in order, clamped and rounded
+ */
+
+static gsize sim_get_spec_bins(gdouble f0, gdouble f1)
+{
+ return (gsize) ((f1 - f0) / sim.radio.freq_inc_hz) + 1;
+}
+
+
+/**
+ * @brief create empty spectral data for a frequency range
+ *
+ * @param f0 the lower bound frequency in Hz
+ * @param f1 the upper bound frequency in Hz
+ *
+ * @returns the zeroed struct spec_data or NULL on error
+ *
+ * @note f0, f1 will be arranged so that always: f0 < f1
+ * @note the range will be clamped to the configured simulator parameters
+ * @note the step size will be fixed to the configured simulator parameters
+ * @note since claming may happen, refer to s->freq_* for further processing
+ */
+
+
+static struct spec_data *sim_create_empty_spectrum(gdouble f0, gdouble f1)
+{
+ gsize bins;
+
+ struct spec_data *s;
+
+ const gsize data_size = sizeof(typeof(*((struct spec_data *) 0)->spec));
+
+
+ sim_order_freq(&f0, &f1);
+
+ f0 = sim_clamp_freq_range(f0);
+ f1 = sim_clamp_freq_range(f1);
+
+ f0 = sim_round_freq_to_res(f0);
+ f1 = sim_round_freq_to_res(f1);
+
+ bins = sim_get_spec_bins(f0, f1);
+
+ if (bins > sim.radio.max_bins) {
+ g_warning(MSG "Computed bins %u > sim.radio.max_bins (%u)",
+ bins, sim.radio.max_bins);
+ }
+
+ s = g_malloc0(sizeof(struct spec_data) + bins * data_size);
+ if (!s)
+ return NULL;
+
+
+ s->n = (typeof(s->n)) bins;
+ s->freq_min_hz = (typeof(s->freq_min_hz)) f0;
+ s->freq_max_hz = (typeof(s->freq_max_hz)) f1;
+ s->freq_inc_hz = (typeof(s->freq_inc_hz)) sim.radio.freq_inc_hz;
+
+ return s;
+}
+
+
+/**
+ * @brief round off a velocity to the HI data resolution
+ *
+ * @param v the velocity (in km/s) to round off
+ *
+ * @returns the rounded off velocity in km/s
+ */
+
+static gdouble HI_round_vel_to_res(gdouble v)
+{
+ v = v / SIM_V_RES_KMS;
+
+ if (v > 0.0)
+ v = ceil(v);
+ else
+ v = floor(v);
+
+ return v * SIM_V_RES_KMS;
+}
+
+
+/**
+ * @brief clamp a velocity into the available HI data range
+ *
+ * @param vel the input velocity (in km/s relative to HI data rest freq)
+ *
+ * @returns vel the range-clamped velocituy
+ */
+
+static gdouble HI_clamp_vel_range(gdouble v)
+{
+ if (v > SIM_V_RED_KMS)
+ return SIM_V_RED_KMS;
+
+ if (v < SIM_V_BLU_KMS)
+ return SIM_V_BLU_KMS;
+
+ return v;
+}
+
+
+/**
+ * @brief get the vlsr-shifted velocity of the spectrum
+ *
+ * @param f the frequency
+ * @param gal the galactic coodrinates
+ *
+ * @returns the doppler velocity clamped into available HI data range
+ *
+ */
+
+static gdouble HI_get_vel(double f, struct coord_galactic gal)
+{
+ gdouble v;
+
+ v = -doppler_vel(f, SIM_V_REF_HZ);
+ v = v - vlsr(galactic_to_equatorial(gal), 0.0);
+
+ v = HI_round_vel_to_res(v);
+
+ return HI_clamp_vel_range(v);
+}
+
+
+/**
+ * @brief geht the bin (array) offset for a given velocity
+ *
+ * @param v the velocity in km/s
+ *
+ * @returns the bin offset in the HI data array for the given velocity
+ *
+ * @note v must be rounded and clamped to HI data range
+ */
+
+static gsize HI_get_vel_bin(gdouble v)
+{
+ gsize bins;
+
+
+ bins = (gsize) (SIM_V_RED_KMS - v);
+
+#if 1
+ /* account for zero-km/s bin */
+ if (v >= 0.0)
+ bins += 1;
+#endif
+ if (bins > SIM_HI_BINS) {
+ g_warning(MSG "something went wrong bins %d > %d max bins",
+ bins, SIM_HI_BINS);
+
+ bins = SIM_HI_BINS;
+ }
+
+ return bins;
+}
+
+
+/**
+ * @brief get the frequency given a velocity for the frequency reference
+ *
+ * @param v a velocity (in km/s)
+ *
+ * @returns the doppler frequency (in Hz)
+ *
+ * @note v must be in the valid HI data range
+ */
+
+static gdouble HI_get_freq(double v)
+{
+ return doppler_freq(v, SIM_V_REF_HZ);
+}
+
+
+/**
+ * @brief fold GLAT into a circular range and round to HI base resolution
+ *
+ * @param lat a galactic latitude
+ *
+ * @returns the galactic latitude folded in range -90.0 to +90.0
+ */
+
+static gdouble HI_fold_glat_round(gdouble lat)
+{
+ lat = fmod(lat + 90.0, 180.0) - 90.0;
+
+ /* the sky is a sphere :) */
+ if (lat < -90.0)
+ lat = lat + 180.0;
+
+ if (lat > 90.0)
+ lat = lat - 180.0;
+
+ /* map to HI data grid */
+ lat = round(( 1.0 / SKY_BASE_RES ) * lat) * SKY_BASE_RES;
+
+
+ return lat;
+}
+
+
+/**
+ * @brief fold GLON into a circular range and round to HI base resolution
+ *
+ * @param lon galactic longitude
+ *
+ * @returns the galactic longitude folded in range 0.0 to +360.0
+ */
+
+static gdouble HI_fold_glon_round(gdouble lon)
+{
+ lon = fmod(lon, 360.0);
+
+ if (lon < 0.0)
+ lon = lon + 360.0;
+ if (lon > 360.0)
+ lon = lon - 360.0;
+
+ /* map to HI data grid */
+ lon = round(( 1.0 / SKY_BASE_RES ) * lon) * SKY_BASE_RES;
+
+
+ return lon;
+}
+
+
+/**
+ * @brief apply scale calibration to HI data for use with spec data
+ *
+ * @param val the value to calibrate
+ *
+ * @returns the calibrated value
+ *
+ * @note spectra data amplitude unit is milliKelvins
+ */
+
+static gdouble HI_raw_amp_to_mKelvins(gdouble val)
+{
+ return val * SIM_HI_AMP_CAL;
+}
+
+
+/**
+ * @brief extract and compute get the spectrum for a given galatic lat/lon
+ * center for a given beam
+ *
+ * @param gal the galactic lat/lon
+ * @param red the red velocity
+ * @param blue the blue velocity
+ * @param beam an nxn array which describes the shape of the beam
+ * @param n the shape of the beam in data bins
+ * @param w the width of the area in degrees (== width of the beam)
+ * @param[out] n_elem the number of elements in the returned array
+ *
+ * @returns the requested spectrum or NULL on error
+ *
+ * @note all parameters must be valid for the given Hi data
+ */
+
+static gdouble *HI_gen_conv_spec(struct coord_galactic gal,
+ gdouble red, gdouble blue,
+ const gdouble *beam, gsize n, gdouble w,
+ gsize *n_elem)
+{
+ gsize i;
+ gsize r, b;
+ gsize bw, bh;
+ gsize bins;
+
+ gdouble x, y;
+ gdouble off, res;
+ gdouble lat, lon;
+ gdouble amp;
+
+ gdouble *spec;
+
+ gint16 *raw;
+
+
+ r = HI_get_vel_bin(red);
+ b = HI_get_vel_bin(blue);
+
+
+ bins = b - r + 1;
+
+ spec = g_malloc0(bins * sizeof(double));
+
+ if (!spec)
+ return NULL;
+
+
+ off = 0.5 * w;
+ res = w / ((gdouble) n - 1.0);
+
+
+ bh = 0;
+
+ for (x = -off; x <= off; x += res) {
+
+ bw = 0;
+
+ for (y = -off; y <= off; y += res) {
+
+ lon = HI_fold_glon_round(gal.lon + x);
+ lat = HI_fold_glat_round(gal.lat + y);
+
+ raw = sim_spec_extract_HI_survey(lat, lon);
+
+ for (i = r; i <= b; i++) {
+ amp = (gdouble) raw[SIM_HI_BINS - i - 1];
+ amp = amp * beam[bh * n + bw];
+ spec[i - r] += amp;
+ }
+
+ g_free(raw);
+
+ bw++; /* beam width index */
+ }
+
+ bh++; /* beam height index */
+ }
+
+ for (i = 0; i < bins; i++)
+ spec[i] = HI_raw_amp_to_mKelvins(spec[i]);
+
+ (*n_elem) = bins;
+
+
+ return spec;
+}
+
+/**
+ * @brief get the index for a given frequency within the target spectrum
+ *
+ * @param s the target spec_data
+ * @param gal the galactic lat/lon
+ * @parma f the frequency in Hz
+ *
+ * @returns the index
+ *
+ * @note all inputs must be valid for the particular configuration
+ */
+
+static gsize HI_get_spec_idx(struct spec_data *s, struct coord_galactic gal,
+ gdouble f)
+{
+ gsize idx;
+
+
+ f = f - (gdouble) s->freq_min_hz + vlsr(galactic_to_equatorial(gal), 0.0);;
+ f = f / (gdouble) s->freq_inc_hz;
+
+ idx = (gsize) f;
+
+
+ /* adjust for actual spectrum */
+ if (idx < 0)
+ idx = 0;
+
+ if (idx > s->n)
+ idx = s->n;
+
+ return (gsize) idx;
+}
+
+
+/**
+ * @brief stack a HI spectrum for a give GLAT/GLON and a beam on a base spectrum
+ *
+ * @param s the target spec_data
+ *
+ * @param gal the galactic lat/lon
+ * @param beam an nxn array which describes the shape of the beam
+ * @param n the shape of the beam in data bins
+ * @param w the width of the area in degrees (== width of the beam)
+ */
+
+static void HI_stack_spec(struct spec_data *s,
+ struct coord_galactic gal,
+ const gdouble *beam, gint n, gdouble w)
+{
+ gsize bins;
+ gsize i, i0, i1;
+
+ gdouble f0, f1;
+ gdouble red, blue;
+
+ gdouble *spec;
+
+
+
+ red = HI_get_vel((gdouble) s->freq_min_hz, gal);
+ blue = HI_get_vel((gdouble) s->freq_max_hz, gal);
+
+ spec = HI_gen_conv_spec(gal, red, blue, beam, n, w, &bins);
+
+ if (!spec) {
+ g_warning(MSG "could not extract HI spectrum");
+ return;
+ }
+
+ /* we need to know where we'll have to put the corresponding velocity
+ * bins within the spectrum
+ */
+
+ f0 = HI_get_freq(red);
+ f1 = HI_get_freq(blue);
+
+ f0 = sim_round_freq_to_res(f0);
+ f1 = sim_round_freq_to_res(f1);
+
+ /* find index offsets */
+ i0 = HI_get_spec_idx(s, gal, f0);
+ i1 = HI_get_spec_idx(s, gal, f1);
+
+
+ for (i = i0; i <= i1; i++)
+ s->spec[i] = (typeof(*s->spec)) ((gdouble) s->spec[i] +
+ spec[i - i0]);
+
+
+ g_free(spec);
+}
+
+
+static gdouble *gauss_2d(gdouble sigma, gdouble r, gdouble res, gint *n);
+
+
+/**
+ * @brief stack the CMB for a give GLAT/GLON and a beam on a base spectrum
+ *
+ * @param s the target spec_data
+ *
+ * @param gal the galactic lat/lon
+ * @param beam an nxn array which describes the shape of the beam
+ * @param n the shape of the beam in data bins
+ * @param w the width of the area in degrees (== width of the beam)
+ *
+ * @todo maybe add a data source (e.g. from Planck, such as found at
+ * https://irsa.ipac.caltech.edu/data/Planck/release_2/all-sky-maps/previews/COM_CMB_IQU-100-fgsub-sevem-field-Int_2048_R2.01_full/index.html
+ * for now, we just use a uniform CMB
+ */
+
+static void sim_stack_cmb(struct spec_data *s,
+ struct coord_galactic gal,
+ const gdouble *beam, gint n, gdouble w)
+{
+ gsize i;
+
+
+ /* add CMB in mK */
+ for (i = 0; i < s->n; i++)
+ s->spec[i] = (typeof(*s->spec)) ((double) s->spec[i] + 2725.);
+}
+
+
+/**
+ * @brief stack the system temperature on a base spectrum
+ *
+ * @param tsys a system temperature in Kelvins
+ *
+ * @param s the target spec_data
+ */
+
+static void sim_stack_tsys(struct spec_data *s, gdouble tsys)
+{
+ gsize i;
+
+
+ tsys = tsys * 1000.0; /* to mK */
+
+ for (i = 0; i < s->n; i++)
+ s->spec[i] = (typeof(*s->spec)) ((gdouble) s->spec[i] + tsys);
+}
+
+
+/**
+ * @brief simulate the sun for a given GLAT/GLON and a beam
+ *
+ * @param gal the galactic lat/lon
+ * @param freq a frequency in Hz
+ * @param beam an nxn array which describes the shape of the beam
+ * @param n the shape of the beam in data bins
+ * @param w the width of the area in degrees (== width of the beam)
+ *
+ * @param returns the solar radio temperature
+ *
+ * @note this is a very simple solar model, where we assume that the
+ * smallest beam is 0.5 deg, so we can use the beam itself as the
+ * sun, ignoring any frequency-radius dependency and limb brightening.
+ * If this simulator is ever to be operated for higher resolutions, or
+ * frequency ranges significantly outside the neutral hydrogen range,
+ * this generator must be updated to do an actual convolution on a
+ * better solar base emission model.
+ *
+ * @warn this assumes that the beam's coordinate system's basis vector is in
+ * the same rotation as the horizon coordinate system
+ *
+ * @todo maybe also add an online data source for up-to-date flux values, e.g.
+ * ftp://ftp.swpc.noaa.gov/pub/lists/radio/rad.txt
+ *
+ */
+
+static gdouble sim_sun(struct coord_galactic gal, gdouble freq,
+ const gdouble *beam, gint n, gdouble w)
+{
+ gsize i;
+
+ gssize x, y;
+
+ gdouble res;
+
+ gdouble D;
+ gdouble A;
+ gdouble l;
+ gdouble T;
+
+ struct coord_galactic sun;
+
+
+ sun = equatorial_to_galactic(sun_ra_dec(0.0));
- dia = gtk_message_dialog_new(NULL,
- GTK_DIALOG_MODAL,
- GTK_MESSAGE_ERROR,
- GTK_BUTTONS_CLOSE,
- "Please place sky_vel.dat in the "
- "directory path this program is executed in.");
+ /* force onto grid */
+ sun.lat = round(( 1.0 / SKY_BASE_RES ) * sun.lat) * SKY_BASE_RES;
+ sun.lon = round(( 1.0 / SKY_BASE_RES ) * sun.lon) * SKY_BASE_RES;
- gtk_dialog_run(GTK_DIALOG(dia));
- gtk_widget_destroy(dia);
+ res = w / (gdouble) n;
+ x = (gssize) ((sun.lat - gal.lat) / res + 0.5 * (gdouble) n);
+ y = (gssize) ((sun.lon - gal.lon) / res + 0.5 * (gdouble) n);
- exit(0);
- /** XXX **/
- g_error("%s: error opening file", __func__);
- return 0;
- }
+ if ((x < 0) || (y < 0))
+ return 0.0;
+ if ((x >= n) || (y >= n))
+ return 0.0;
+ l = 299792458.0 / freq;
- map = g_malloc(VEL * SKY_WIDTH * SKY_HEIGHT * sizeof(guint16));
+ D = 1.22 * l / atan(2.0 * sim.r_beam / 180.0 * M_PI);
- fread(map, sizeof(guint16), VEL * SKY_WIDTH * SKY_HEIGHT, data);
- fclose(data);
- }
+ A = 0.25 * D * D * M_PI;
- as = g_malloc(VEL * sizeof(guint16));
- offset = get_offset(glat, glon);
+ /* convert flux to temperature
+ * sun is in solar flux units, i.e. 10^-22 W/m^2/Hz
+ */
- memcpy(as, &map[offset], sizeof(guint16) * VEL);
+ T = A * 0.5 * 1e-22 * sim.sun_sfu / 1.38064852e-23;
- return as;
+
+ /* beam is normalised to integral 1, scale temperature to center bin */
+ T = T * beam[y * n + x] / beam[n/2 * n + n/2];
+
+ return T;
}
/**
- * @brief build a gaussian convolution kernel for a beam of radius r
+ * @brief stack the sun for a given GLAT/GLON and a beam on a base spectrum
+ *
+ * @param s the target spec_data
*
- * @param r the beam radius in degrees
+ * @param gal the galactic lat/lon
+ * @param beam an nxn array which describes the shape of the beam
+ * @param n the shape of the beam in data bins
+ * @param w the width of the area in degrees (== width of the beam)
*/
-static gpointer gauss_kernel_spec(double r, double glat, double glon)
+
+static void sim_stack_sun(struct spec_data *s, struct coord_galactic gal,
+ const gdouble *beam, gint n, gdouble w)
{
- static double x, x0, x1;
-static double bins;
+ gsize i;
+
+ gdouble T;
+
+
+ T = sim_sun(gal, s->freq_min_hz, beam, n, w);
- double binw;
+ /* convert to mK */
+ T = 1000. * T;
- int i, n;
+ for (i = 0; i < s->n; i++)
+ s->spec[i] = (typeof(*s->spec)) ((double) s->spec[i] + T);
+}
- double s;
-static double intg, prev, tmp, total;
+/**
+ * @brief simulate moon for a given GLAT/GLON and a beam
+ *
+ * @param gal the galactic lat/lon
+ * @param beam an nxn array which describes the shape of the beam
+ * @param n the shape of the beam in data bins
+ * @param w the width of the area in degrees (== width of the beam)
+ *
+ * @param returns the moon's radio temperature
+ *
+ * @note this is a very simple moon model, where we assume that the
+ * smallest beam is 0.5 deg, so we can use the beam itself as the
+ * moon
+ *
+ * @warn this assumes that the beam's coordinate system's basis vector is in
+ * the same rotation as the horizon coordinate system
+ *
+ */
-static double *conv;
-static double *conv2;
+static gdouble sim_moon(struct coord_galactic gal, const gdouble *beam,
+ gint n, gdouble w)
+{
+ gsize i;
- /* our width-normalised kernel runs from -3 to +3 sigma */
+ gssize x, y;
- x1 = 3.0 / 2.0 * sqrt(2.0 * log(2.0)) * r;
+ gdouble res;
- x0 = - x1;
+ gdouble E;
+ gdouble T;
- /* need an extra center bin, does not count towards stepsize */
- bins = round((x1 - x0) / SKY_BASE_RES) + 1.0;
- binw = (x1 - x0) / (bins - 1.0);
+ gdouble a, b;
+ gdouble costheta;
+ struct coord_galactic moon;
+ struct coord_galactic sun;
- if(conv)
- g_free(conv);
- conv = g_malloc(bins * sizeof(double));
+ moon = equatorial_to_galactic(moon_ra_dec(server_cfg_get_station_lat(),
+ server_cfg_get_station_lon(),
+ 0.0));
+ /* force onto grid */
+ moon.lat = round(( 1.0 / SKY_BASE_RES ) * moon.lat) * SKY_BASE_RES;
+ moon.lon = round(( 1.0 / SKY_BASE_RES ) * moon.lon) * SKY_BASE_RES;
- x = x0;
- total = 0.0;
- for (i = 0; i < bins; i++) {
- n = 0;
+ res = w / (gdouble) n;
- s = x - binw * 0.5;
+ x = (gssize) ((moon.lat - gal.lat) / res + 0.5 * (gdouble) n);
+ y = (gssize) ((moon.lon - gal.lon) / res + 0.5 * (gdouble) n);
- tmp = 0.0;
- intg = 0.0;
- prev = gauss_norm(s, r);
- do {
- s += SKY_GAUSS_INTG_STP;
- /* simple numerical integration */
- intg += SKY_GAUSS_INTG_STP * 0.5 * ((gauss_norm(s, r) + prev));
- prev = intg;
- } while (s < (x + binw * 0.5));
+ if ((x < 0) || (y < 0))
+ return 0.0;
- conv[i] = prev;
+ if ((x >= n) || (y >= n))
+ return 0.0;
- total += intg;
+ /* approximate the reflected solar energy per unit area via the
+ * solar constant in Earth's neighbourhood, reduced by the angle
+ * between the surface of the moon facing the earth and the sun
+ * (varname I is reserved due to complex.h)
+ */
- x += binw;
+ sun = equatorial_to_galactic(sun_ra_dec(0.0));
- }
+ a = sun.lat - moon.lat;
+ b = sun.lon - moon.lon;
- /* oder conv[i] dividieren! */
- gdouble norm = 2.0 * bins * total;
+ costheta = (sun.lon * sun.lon + sun.lat * sun.lat) / euclid_norm(a, b);
+ costheta = fmod(costheta - 180.0, 360.) / 180. * M_PI;
- if (conv2)
- g_free(conv2);
- conv2 = g_malloc(bins * bins * sizeof(double));
- int j;
- for (i = 0; i < (int) bins; i++)
- for (j = 0; j < (int) bins; j++) {
- conv2[i* (int) bins + j] = (conv[i] + conv[j]) / norm;
- }
+ /* average the moon temperature across all degrees of latitude by
+ * divding by pi
+ */
+ E = 1366.0 * costheta / M_PI;
- /* strategy: allocate columns * rows in flat array, fill with
- * convolution kernel rows (2d gaussian has radial symmetry and can be convolve
- * by application to rows+columns; add data in same fashion, 2d-indexed ragged type:
- * compute once, apply continuously;
- * once, convolved, sum along 2d direction for final spectrum
- *
+ /* for emissivity ~ reflectance, we get the Moon's surface temperature
+ * via the Stefan-Boltzmann equation
*/
+ T = pow((E / 5.67e-8), 0.25);
- gdouble *res;
- static gint16 *raw;
+ /* scale to beam/moon ratio (apparent diameter ~0.5 deg) */
- double l, b;
- double dl, db;
- gdouble lat, lon;
+ T = T * 0.25 * 0.25 / (sim.r_beam * sim.r_beam);
+ /* beam is normalised to integral 1, scale temperature to center bin */
+ T = T * beam[y * n + x] / beam[n/2 * n + n/2];
- res = g_malloc0(VEL * sizeof(double));
+ return T;
+}
- gint idxl=0;
- gint idxb=0;
- /* coordinates needed */
- for (l = x0; l <= x1; l+=binw) {
- idxb = 0;
- for (b = x0; b <= x1; b+=binw) {
+/**
+ * @brief stack the moon for a given GLAT/GLON and a beam on a base spectrum
+ *
+ * @param s the target spec_data
+ *
+ * @param gal the galactic lat/lon
+ * @param beam an nxn array which describes the shape of the beam
+ * @param n the shape of the beam in data bins
+ * @param w the width of the area in degrees (== width of the beam)
+ *
+ * @note this is a very simple moon model, where we assume that the
+ * smallest beam is 0.5 deg, so we can use the beam itself as the
+ * moon
+ *
+ * @warn this assumes that the beam's coordinate system's basis vector is in
+ * the same rotation as the horizon coordinate system
+ *
+ */
- /* get needed coordinate */
- lat = glat + b;
- lon = glon + l;
+static void sim_stack_moon(struct spec_data *s, struct coord_galactic gal,
+ const gdouble *beam, gint n, gdouble w)
+{
+ gsize i;
+ gdouble T;
- /* fold lat into 180 deg range */
- lat = fmod(lat + 90.0, 180.0) - 90.0;
- /* same for lon */
- lon = fmod(lon, 360.0);
+ T = sim_moon(gal, beam, n, w);
- /* the sky is a sphere :) */
- if (lat < -90.0)
- lat = lat + 180.0;
- if (lat > 90.0)
- lat = lat - 180.0;
+ /* convert to mK */
+ T = 1000. * T;
- if (lon < 0.0)
- lon = lon + 360.0;
- if (lon > 360.0)
- lon = lon - 360.0;
+ for (i = 0; i < s->n; i++)
+ s->spec[i] = (typeof(*s->spec)) ((double) s->spec[i] + T);
+}
- /* data is reverse? */
- //lon = 360.0 - lon;
+/**
+ * @brief stack system efficiency on a base spectrum
+ *
+ * @param s the target spec_data
+ * @param eff a efficiency multiplier
+ */
- /***
- * strategy for image creation: sum up vel range first,
- * then store. Apply convolution to stored data
- */
+static void sim_stack_eff(struct spec_data *s, gdouble eff)
+{
+ gsize i;
+ for (i = 0; i < s->n; i++)
+ s->spec[i] = (typeof(*s->spec)) ((gdouble) s->spec[i] * eff);
+}
- /* map to data grid */
- lat = round(( 1.0 / SKY_BASE_RES ) * lat) * SKY_BASE_RES;
- lon = round(( 1.0 / SKY_BASE_RES ) * lon) * SKY_BASE_RES;
+/**
+ * @brief generate approximate gaussian noise
+ *
+ * @note we want something decently fast, this Box-Muller transform
+ * appears to work quite well
+ */
- raw = sim_spec_extract_HI_survey(lat, lon);
- /* get the corresponding coordinate */
+static gdouble sim_rand_gauss(void)
+{
+ static gboolean gen;
+ static gdouble u, v;
- for (i = 0; i < VEL; i++)
- res[i] += ((gdouble) raw[i]) * conv2[idxl * (int) bins + idxb];
+ gen = !gen;
- g_free(raw);
- idxb++;
- }
- idxl++;
- }
+ if (!gen)
+ return sqrt(- 2.0 * log(u)) * cos(2.0 * M_PI * v);
+
+ u = (rand() + 1.0) / (RAND_MAX + 2.0);
+ v = rand() / (RAND_MAX + 1.0);
+
+ return sqrt(-2.0 * log(u)) * sin(2.0 * M_PI * v);
+}
+
+
+
+/**
+ * @brief stack system noise on a base spectrum
+ *
+ * @param s the target spec_data
+ * @param sig a sigma mutliplier for the noise
+ *
+ * @note the noise is scaled by the sqrt() of the signal amplitude
+ */
+
+static void sim_stack_gnoise(struct spec_data *s, gdouble sig)
+{
+ gsize i;
+
+ gdouble amp;
-#if 0
- for (i = 0; i < VEL; i++)
- res[i] /= norm;
-#endif
- return res;
+ for (i = 0; i < s->n; i++) {
+ amp = (gdouble) s->spec[i];
+ amp = amp + sqrt(amp) * sig * sim_rand_gauss();
+ s->spec[i] = (typeof(*s->spec)) amp;
+ }
}
+
/**
* @brief transpose a two-dimensional array of complex doubles
*
@@ -1188,100 +1901,6 @@ static gdouble *sim_rt_get_HI_img(gint vmin, gint vmax)
}
-/**
- * @brief returns an NxN normalised gaussian convolution kernel for a given
- * sigma cutoff and resolution in degrees
- *
- * @note bins is always uneven (need center bin)
- */
-
-static gdouble *gauss_conv_kernel(gdouble r, gdouble sigma, gdouble res,
- gint *n)
-{
- gint i, j;
- gint bins;
-
- gdouble x;
- gdouble x0;
- gdouble x1;
-
- gdouble binw;
-
- gdouble tot;
- gdouble norm;
-
- double *conv;
-
- double *conv2;
-
-
-
- /* our width-normalised kernel runs from -sigma to +sigma */
- x1 = sigma / 2.0 * sqrt(2.0 * log(2.0)) * r;
- x0 = - x1;
- x = x0;
-
- /* need an extra center bin, does not count towards stepsize */
- bins = (gint) (round((x1 - x0) / res));
-
- if (!(bins & 0x1))
- bins++;
-
-
- binw = (x1 - x0) / ((gdouble) (bins - 1));
-
-
- /* 1d convolution profile */
- conv = g_malloc(bins * sizeof(gdouble));
-
-
- tot = 0.0;
-
- for (i = 0; i < bins; i++) {
-
- gdouble s;
- gdouble intg;
- gdouble prev;
-
- s = x - binw * 0.5;
-
- intg = 0.0;
- prev = gauss_norm(s, r);
-
- /* simple numerical integration */
- do {
- s += SKY_GAUSS_INTG_STP;
- intg += SKY_GAUSS_INTG_STP * 0.5 * ((gauss_norm(s, r) + prev));
- prev = intg;
- } while (s < (x + binw * 0.5));
-
- conv[i] = prev;
-
- tot += intg;
-
- x += binw;
- }
-
-
- /* 2d kernel */
- conv2 = g_malloc(bins * bins * sizeof(double));
-
- norm = 2.0 * tot * bins;
-
- for (i = 0; i < bins; i++)
- for (j = 0; j < bins; j++)
- conv2[i * bins + j] = (conv[i] + conv[j]) / norm;
-
-
- g_free(conv);
-
-
- (*n) = bins;
-
- return conv2;
-}
-
-
/**
* @brief returns an NxN normalised gaussian convolution kernel for a given
* sigma cutoff and resolution in degrees
@@ -1295,7 +1914,7 @@ static gdouble *gauss_conv_kernel(gdouble r, gdouble sigma, gdouble res,
* (for no added visual benefit)
*/
-static gdouble *gauss_2d(gdouble r, gdouble sigma, gdouble res, gint *n)
+static gdouble *gauss_2d(gdouble sigma, gdouble r, gdouble res, gint *n)
{
gint i, j;
gint bins;
@@ -1311,7 +1930,7 @@ static gdouble *gauss_2d(gdouble r, gdouble sigma, gdouble res, gint *n)
/* our width-normalised kernel runs from -sigma to +sigma */
- x = sigma / 2.0 * sqrt(2.0 * log(2.0)) * r;
+ x = gauss_half_width_sigma_r(sigma, r);
/* need an extra center bin, does not count towards stepsize */
bins = (gint) (round((2.0 * x) / res));
@@ -1340,7 +1959,7 @@ static gdouble *gauss_2d(gdouble r, gdouble sigma, gdouble res, gint *n)
}
/* normalise kernel */
- for (i = 0; i < bins; i++)
+ for (i = 0; i < bins * bins; i++)
ker[i] /= sum;
@@ -1483,10 +2102,18 @@ static void sim_gen_sky(void)
vmax = v2;
}
+ kernel = gauss_2d(3.0, sim.r_beam, SKY_BASE_RES, &n);
+
if (!conv)
conv = g_malloc(SKY1_WIDTH * SKY1_HEIGHT * sizeof(double complex));
if (!raw_sky) {
+#if 0
+ gdouble T;
+ struct coord_galactic gal;
+
+ gdouble w = 2.0 * gauss_half_width_sigma_r(3.0, sim.r_beam);
+#endif
raw_sky = sim_rt_get_HI_img(vmin, vmax);
if (csky) {
@@ -1494,6 +2121,45 @@ static void sim_gen_sky(void)
csky =NULL;
}
+#if 0
+ gal = equatorial_to_galactic(moon_ra_dec(server_cfg_get_station_lat(),
+ server_cfg_get_station_lon(),
+ 0.0));
+
+ T = sim_moon(gal, kernel, n, w);
+
+ T = T * VEL * 100.;
+
+ /* force onto grid */
+ gal.lat = round(( 1.0 / SKY_BASE_RES ) * gal.lat) * SKY_BASE_RES;
+ gal.lon = round(( 1.0 / SKY_BASE_RES ) * gal.lon) * SKY_BASE_RES;
+
+ gal.lat = 90. - gal.lat;
+ gal.lon = gal.lon - 90.;
+
+ raw_sky[SKY1_WIDTH * ((int)gal.lat) * 2 + (int)gal.lon * 2] += T;
+
+
+
+ gal = equatorial_to_galactic(sun_ra_dec(0.0));
+
+ T = sim_sun(gal, SIM_V_REF_HZ, kernel, n, w);
+
+ T = T * VEL * 100.;
+ g_message("at %g %g T %g", gal.lat, gal.lon, T);
+
+ /* force onto grid */
+ gal.lat = round(( 1.0 / SKY_BASE_RES ) * gal.lat) * SKY_BASE_RES;
+ gal.lon = round(( 1.0 / SKY_BASE_RES ) * gal.lon) * SKY_BASE_RES;
+
+ gal.lat = gal.lat + 90.;
+ gal.lon = gal.lon - 90.;
+ g_message("at %g %g T %g", gal.lat, gal.lon, T);
+
+ raw_sky[SKY1_WIDTH * ((int)gal.lat) * 2 + (int)gal.lon * 2] += T;
+#endif
+
+
csky = g_malloc0(SKY1_WIDTH * SKY1_HEIGHT * sizeof(double complex));
double complex *tmp = g_malloc0(SKY1_WIDTH * SKY1_HEIGHT * sizeof(double complex));
@@ -1505,7 +2171,6 @@ static void sim_gen_sky(void)
#if 1
- kernel = gauss_2d(sim.r_beam, 3.0, SKY_BASE_RES, &n);
#else
kernel = gauss_conv_kernel(sim.r_beam, 3.0, SKY_BASE_RES, &n);
#endif
@@ -1635,6 +2300,23 @@ static gboolean sim_spb_eff_value_changed_cb(GtkSpinButton *sb, gpointer data)
sim.eff = gtk_spin_button_get_value(sb);
}
+static gboolean sim_spb_sig_n_value_changed_cb(GtkSpinButton *sb, gpointer data)
+{
+ sim.sig_n = gtk_spin_button_get_value(sb);
+}
+
+static gboolean sim_spb_readout_hz_value_changed_cb(GtkSpinButton *sb, gpointer data)
+{
+ sim.readout_hz = gtk_spin_button_get_value(sb);
+}
+
+static gboolean sim_sun_sfu_value_changed_cb(GtkSpinButton *sb, gpointer data)
+{
+ sim.sun_sfu = gtk_spin_button_get_value(sb);
+}
+
+
+
static void sim_rt_gui_defaults(void)
{
@@ -1699,12 +2381,30 @@ static GtkWidget *sim_rt_par_gui(void)
g_signal_connect(GTK_SPIN_BUTTON(w), "value-changed",
G_CALLBACK(sim_spb_tsys_value_changed_cb), NULL);
+ w = gtk_label_new("SIG");
+ gtk_widget_set_halign(w, GTK_ALIGN_START);
+ gtk_widget_set_halign(w, GTK_ALIGN_START);
+ gtk_label_set_xalign(GTK_LABEL(w), 0.0);
+ gtk_grid_attach(GTK_GRID(grid), w, 0, 2, 1, 1);
+
+ w = gtk_spin_button_new_with_range(0., 20.0, 0.1);
+ gtk_entry_set_alignment(GTK_ENTRY(w), 1.0);
+ gtk_spin_button_set_numeric(GTK_SPIN_BUTTON(w), TRUE);
+ gtk_spin_button_set_digits(GTK_SPIN_BUTTON(w), 1);
+ gtk_spin_button_set_value(GTK_SPIN_BUTTON(w), SIM_SIG_NOISE);
+ gtk_widget_set_halign(w, GTK_ALIGN_FILL);
+ gtk_widget_set_hexpand(w, FALSE);
+ gtk_grid_attach(GTK_GRID(grid), w, 1, 2, 1, 1);
+
+ g_signal_connect(GTK_SPIN_BUTTON(w), "value-changed",
+ G_CALLBACK(sim_spb_sig_n_value_changed_cb), NULL);
+
w = gtk_label_new("EFF");
gtk_widget_set_halign(w, GTK_ALIGN_START);
gtk_widget_set_halign(w, GTK_ALIGN_START);
gtk_label_set_xalign(GTK_LABEL(w), 0.0);
- gtk_grid_attach(GTK_GRID(grid), w, 0, 2, 1, 1);
+ gtk_grid_attach(GTK_GRID(grid), w, 0, 3, 1, 1);
w = gtk_spin_button_new_with_range(0., 1.0, 0.1);
gtk_entry_set_alignment(GTK_ENTRY(w), 1.0);
@@ -1713,7 +2413,7 @@ static GtkWidget *sim_rt_par_gui(void)
gtk_spin_button_set_value(GTK_SPIN_BUTTON(w), SIM_EFF);
gtk_widget_set_halign(w, GTK_ALIGN_FILL);
gtk_widget_set_hexpand(w, FALSE);
- gtk_grid_attach(GTK_GRID(grid), w, 1, 2, 1, 1);
+ gtk_grid_attach(GTK_GRID(grid), w, 1, 3, 1, 1);
g_signal_connect(GTK_SPIN_BUTTON(w), "value-changed",
G_CALLBACK(sim_spb_eff_value_changed_cb), NULL);
@@ -1724,7 +2424,7 @@ static GtkWidget *sim_rt_par_gui(void)
gtk_widget_set_halign(w, GTK_ALIGN_START);
gtk_widget_set_halign(w, GTK_ALIGN_START);
gtk_label_set_xalign(GTK_LABEL(w), 0.0);
- gtk_grid_attach(GTK_GRID(grid), w, 0, 3, 1, 1);
+ gtk_grid_attach(GTK_GRID(grid), w, 0, 4, 1, 1);
w = gtk_spin_button_new_with_range(-90., 90., 0.01);
gtk_entry_set_alignment(GTK_ENTRY(w), 1.0);
@@ -1733,7 +2433,7 @@ static GtkWidget *sim_rt_par_gui(void)
gtk_spin_button_set_value(GTK_SPIN_BUTTON(w), server_cfg_get_station_lat());
gtk_widget_set_halign(w, GTK_ALIGN_FILL);
gtk_widget_set_hexpand(w, FALSE);
- gtk_grid_attach(GTK_GRID(grid), w, 1, 3, 1, 1);
+ gtk_grid_attach(GTK_GRID(grid), w, 1, 4, 1, 1);
g_signal_connect(GTK_SPIN_BUTTON(w), "value-changed",
@@ -1743,7 +2443,7 @@ static GtkWidget *sim_rt_par_gui(void)
gtk_widget_set_halign(w, GTK_ALIGN_START);
gtk_widget_set_halign(w, GTK_ALIGN_START);
gtk_label_set_xalign(GTK_LABEL(w), 0.0);
- gtk_grid_attach(GTK_GRID(grid), w, 0, 4, 1, 1);
+ gtk_grid_attach(GTK_GRID(grid), w, 0, 5, 1, 1);
w = gtk_spin_button_new_with_range(-180., 180., 0.01);
gtk_entry_set_alignment(GTK_ENTRY(w), 1.0);
@@ -1752,12 +2452,49 @@ static GtkWidget *sim_rt_par_gui(void)
gtk_spin_button_set_value(GTK_SPIN_BUTTON(w), server_cfg_get_station_lon());
gtk_widget_set_halign(w, GTK_ALIGN_FILL);
gtk_widget_set_hexpand(w, FALSE);
- gtk_grid_attach(GTK_GRID(grid), w, 1, 4, 1, 1);
+ gtk_grid_attach(GTK_GRID(grid), w, 1, 5, 1, 1);
g_signal_connect(GTK_SPIN_BUTTON(w), "value-changed",
G_CALLBACK(sim_spb_lon_value_changed_cb), NULL);
+ w = gtk_label_new("Rate [Hz]");
+ gtk_widget_set_halign(w, GTK_ALIGN_START);
+ gtk_widget_set_halign(w, GTK_ALIGN_START);
+ gtk_label_set_xalign(GTK_LABEL(w), 0.0);
+ gtk_grid_attach(GTK_GRID(grid), w, 0, 6, 1, 1);
+
+ w = gtk_spin_button_new_with_range(0.1, 32.0, 1.0);
+ gtk_entry_set_alignment(GTK_ENTRY(w), 1.0);
+ gtk_spin_button_set_numeric(GTK_SPIN_BUTTON(w), TRUE);
+ gtk_spin_button_set_digits(GTK_SPIN_BUTTON(w), 1);
+ gtk_spin_button_set_value(GTK_SPIN_BUTTON(w), SIM_READOUT_HZ);
+ gtk_widget_set_halign(w, GTK_ALIGN_FILL);
+ gtk_widget_set_hexpand(w, FALSE);
+ gtk_grid_attach(GTK_GRID(grid), w, 1, 6, 1, 1);
+
+ g_signal_connect(GTK_SPIN_BUTTON(w), "value-changed",
+ G_CALLBACK(sim_spb_readout_hz_value_changed_cb), NULL);
+
+ w = gtk_label_new("Sun [SFU]");
+ gtk_widget_set_halign(w, GTK_ALIGN_START);
+ gtk_widget_set_halign(w, GTK_ALIGN_START);
+ gtk_label_set_xalign(GTK_LABEL(w), 0.0);
+ gtk_grid_attach(GTK_GRID(grid), w, 0, 7, 1, 1);
+
+ w = gtk_spin_button_new_with_range(1., 200., 1.0);
+ gtk_entry_set_alignment(GTK_ENTRY(w), 1.0);
+ gtk_spin_button_set_numeric(GTK_SPIN_BUTTON(w), TRUE);
+ gtk_spin_button_set_digits(GTK_SPIN_BUTTON(w), 1);
+ gtk_spin_button_set_value(GTK_SPIN_BUTTON(w), SIM_SUN_SFU);
+ gtk_widget_set_halign(w, GTK_ALIGN_FILL);
+ gtk_widget_set_hexpand(w, FALSE);
+ gtk_grid_attach(GTK_GRID(grid), w, 1, 7, 1, 1);
+
+ g_signal_connect(GTK_SPIN_BUTTON(w), "value-changed",
+ G_CALLBACK(sim_sun_sfu_value_changed_cb), NULL);
+
+
return frame;
}
@@ -1981,11 +2718,9 @@ static uint32_t sim_spec_acquire(struct observation *obs)
struct coord_horizontal hor;
struct coord_galactic gal;
- gdouble *rawspec;
gint16 *rawu;
- //sim_gen_sky();
#if 0
if (!obs->acq.acq_max)
return 0;
@@ -1994,100 +2729,75 @@ static uint32_t sim_spec_acquire(struct observation *obs)
hor.el = sim.el.cur;
gal = horizontal_to_galactic(hor, server_cfg_get_station_lat(), server_cfg_get_station_lon());
-#if 0
- gal.lat = 0;
- gal.lon = 180.0;
-#endif
+
st.busy = 1;
st.eta_msec = (typeof(st.eta_msec)) 20.0; /* whatever */
ack_status_rec(PKT_TRANS_ID_UNDEF, &st);
-#if 1
/* conv */
+ static gdouble *beam;
+ static gdouble r;
+ static gdouble sky_deg;
+ static gsize n_vel;
+ static gint n_beam;
+ static gdouble glat;
+ static gdouble glon;
- int v1 = (int) DOPPLER_VEL(g_obs.acq.freq_stop_hz, SIM_V_REF_MHZ * 1e6) + SIM_V_RED_KMS +2 + vlsr(galactic_to_equatorial(gal), 0.0);
- int v0 = (int) DOPPLER_VEL(g_obs.acq.freq_start_hz, SIM_V_REF_MHZ * 1e6) + SIM_V_RED_KMS + vlsr(galactic_to_equatorial(gal), 0.0);
-
-// g_message("vmin: %d vmax %d", v0, v1);
- if (v0 < 0)
- v0 = 0;
- if (v1 > VEL)
- v1 = VEL;
-// g_message("vmin: %d vmax %d\n", v0, v1);
+ gal.lat = round(( 1.0 / SKY_BASE_RES ) * gal.lat) * SKY_BASE_RES;
+ gal.lon = round(( 1.0 / SKY_BASE_RES ) * gal.lon) * SKY_BASE_RES;
+ if (r != sim.r_beam || !beam) {
+ r = sim.r_beam;
- /* prepare and send: allocate full length */
- s = g_malloc0(sizeof(struct spec_data) + (v1 - v0) * sizeof(uint32_t));
-
- rawspec = gauss_kernel_spec(sim.r_beam, gal.lat, gal.lon);
-#if 0
- s->freq_min_hz = (typeof(s->freq_min_hz)) DOPPLER_FREQ((v0 - SIM_V_RED_KMS ), SIM_V_REF_MHZ * 1e6);
- s->freq_max_hz = (typeof(s->freq_max_hz)) DOPPLER_FREQ((v1 - SIM_V_RED_KMS - 2 ), SIM_V_REF_MHZ * 1e6);
-#endif
- s->freq_min_hz = (typeof(s->freq_min_hz)) g_obs.acq.freq_start_hz;
- s->freq_max_hz = (typeof(s->freq_max_hz)) g_obs.acq.freq_stop_hz;
- s->freq_inc_hz = (typeof(s->freq_inc_hz)) SIM_FREQ_STP_HZ;
-#if 1
- s->freq_min_hz = round(s->freq_min_hz / s->freq_inc_hz) * s->freq_inc_hz;
- s->freq_max_hz = round(s->freq_max_hz / s->freq_inc_hz) * s->freq_inc_hz;
-#endif
-
- //srand(time(0));
- for (i = v0; i < v1; i++) {
- /* reverse or not? */
- s->spec[i- v0] = (uint32_t) ((double) rawspec[VEL - i - 1] * 10. * sim.eff + sim.tsys * 1000.); /* to mK fom cK + tsys*/
- // s->spec[i] = (uint32_t) rawspec[i] *10 + 200000; /* to mK fom cK + tsys*/
- s->spec[i-v0] += (int) (GNOISE * sqrt((double) s->spec[i-v0] * 15.));
+ if (beam) {
+ g_free(beam);
+ beam = NULL;
+ }
- s->n++;
+ sky_deg = 2.0 * gauss_half_width_sigma_r(3.0, sim.r_beam);
+ beam = gauss_2d(3.0, sim.r_beam, SKY_BASE_RES, &n_beam);
}
- g_free(rawspec);
+ if ((glat != gal.lat) || (glon != gal.lon))
+ glat = gal.lat;
- /* handover for transmission */
- ack_spec_data(PKT_TRANS_ID_UNDEF, s);
+ s = sim_create_empty_spectrum(g_obs.acq.freq_start_hz,
+ g_obs.acq.freq_stop_hz);
-#else
- /* noconv */
+ if (!s) {
+ g_warning(MSG "could not create spectral data");
+ return 0;
+ }
+ /* NOTE: values in s->spec must ALWAYS be >0, since it is (usually)
+ * a uin32_t. Adding the CMB temperature is a good place to start
+ */
+ sim_stack_cmb(s, gal, beam, n_beam, sky_deg);
- /* prepare and send: allocate full length */
- s = g_malloc0(sizeof(struct spec_data) + VEL * sizeof(uint32_t));
- rawu = sim_spec_extract_HI_survey(gal.lat, gal.lon);
+ HI_stack_spec(s, gal, beam, n_beam, sky_deg);
+ sim_stack_moon(s, gal, beam, n_beam, sky_deg);
- s->freq_min_hz = (typeof(s->freq_min_hz)) DOPPLER_FREQ((SIM_V_RED_KMS + vlsr(galactic_to_equatorial(gal), 0.0)), SIM_V_REF_MHZ * 1e6);
- s->freq_max_hz = (typeof(s->freq_max_hz)) DOPPLER_FREQ((SIM_V_BLU_KMS + vlsr(galactic_to_equatorial(gal), 0.0)), SIM_V_REF_MHZ * 1e6);
- s->freq_inc_hz = (typeof(s->freq_inc_hz)) SIM_FREQ_STP_HZ;
+ sim_stack_sun(s, gal, beam, n_beam, sky_deg);
- //srand(time(0));
- for (i = 0; i < VEL; i++) {
- s->spec[i] = (uint32_t) rawu[VEL - i - 1] * 10 + 300000; /* to mK from cK + tsys*/
- // s->spec[i] += (int) (GNOISE * 1000.);
+ sim_stack_eff(s, sim.eff);
- s->n++;
- }
+ sim_stack_tsys(s, sim.tsys);
- g_free(rawu);
+ sim_stack_gnoise(s, sim.sig_n);
/* handover for transmission */
ack_spec_data(PKT_TRANS_ID_UNDEF, s);
-#endif
-
-
-
-
-
st.busy = 0;
st.eta_msec = 0;
@@ -2098,7 +2808,7 @@ static uint32_t sim_spec_acquire(struct observation *obs)
cleanup:
g_free(s);
- g_usleep(G_USEC_PER_SEC / 32);
+ g_usleep(G_USEC_PER_SEC / sim.readout_hz);
return obs->acq.acq_max;
}