Commit b7fbea7b authored by Daniel Brown's avatar Daniel Brown
Browse files

making it more clear in the code that BS alpha is just for side 1, also...

making it more clear in the code that BS alpha is just for side 1, also storing side 2 alpha. Fixing error where fsig was not using correct alphas on each side
parent f7848973
......@@ -1083,10 +1083,12 @@ typedef struct beamsplitter {
double phi; //!< beam splitter tuning \see mirror_t::phi
double x_off;
double y_off;
double alpha; //!< angle of incidence
double ir1; //!< refractive index (where???)
double ir2; //!< refractive index
double alpha_1; //!< angle of incidence set by user [deg]
double alpha_2; //!< angle of incidence on side 2 calculated during ABCD calculations [deg]
double ir1; //!< refractive index side 1 computed during ABCD calculations
double ir2; //!< refractive index side 2 computed during ABCD calculations
transfer_func_t *mech_tf; // index of mechanical transfer function applied to this mirror
......
......@@ -3391,7 +3391,6 @@ void set_ABCD_beamsplitter(int bs_index) {
node_t node1, node2, node3, node4;
double n1, n2;
double co1, co2, si1, si2, coco;
double alpha2;
double ir[2] = {0};
beamsplitter_t *bs;
......@@ -3464,25 +3463,25 @@ void set_ABCD_beamsplitter(int bs_index) {
bs->ir1 = n1;
bs->ir2 = n2;
si1 = sin(bs->alpha * RAD);
si1 = sin(bs->alpha_1 * RAD);
si2 = n1 / n2*si1;
if (fabs(si2) > 1) { // fabs 290402
gerror("total reflection at beamsplitter %s\n", bs->name);
}
alpha2 = DEG * asin(si2);
bs->alpha_2 = DEG * asin(si2);
co1 = cos(bs->alpha * RAD);
co2 = cos(alpha2 * RAD);
co1 = cos(bs->alpha_1 * RAD);
co2 = cos(bs->alpha_2 * RAD);
if (inter.debug & 256) {
fprintf(stdout, "Setting ABCD for bs %s\n", bs->name);
fprintf(stdout, "%s: alpha %g , alpha2 %g\n", bs->name, bs->alpha, alpha2);
fprintf(stdout, "%s: alpha %g , alpha2 %g\n", bs->name, bs->alpha_1, bs->alpha_2);
fprintf(stdout, "n1 %g n2 %g c1 %g c2 %g si %g s2 %g\n", n1, n2, co1, co2, si1, si2);
}
if (co1 == 0.0 || co2 == 0.0) {
gerror("beamsplitter '%s': alpha(%s) out of range\n", bs->name,
double_form(bs->alpha));
double_form(bs->alpha_1));
}
coco = n2 * co2 - n1*co1;
......
......@@ -1641,7 +1641,6 @@ void check_for_long_line(FILE *fp, char *s) {
}
}
//! Prepare input line for processing
/*!
......
......@@ -649,14 +649,7 @@ complex_t xyamp(node_t *sig_node, complex_t q) {
}
void calc_beamsplitter_signal_vars( beamsplitter_t *beamsplitter, double *factor, double *ph1, double f, double f_c){
// TODO check refractive index
*factor *= 2.0 * (1.0 + f_c / inter.f0)
* sqrt(beamsplitter->R)
* cos(beamsplitter->alpha * RAD);
*ph1 = beamsplitter->phi
* ((1.0 + f / inter.f0) + (1.0 + f_c / inter.f0))
* cos(beamsplitter->alpha * RAD);
}
complex_t calc_reflect_tilt_signal(signal_t *signal, complex_t **knm, int in, int out, complex_t X_bar){
......@@ -766,15 +759,19 @@ void fill_bs_signal_rhs(signal_t *signal, ifo_matrix_vars_t* M_car, ifo_matrix_v
for(f=0; f<M_sig->num_frequencies; f++){
frequency_t *f_sig = M_sig->frequencies[f];
double phase = (90.0 + f_sig->order * signal->phase);
double factor = signal->amplitude * EPSILON * 0.5; //Epsilon/2 for fsig
double turn, ph1;
double _factor = signal->amplitude * EPSILON * 0.5; //Epsilon/2 for fsig
double turn, _ph1;
if(signal->type == SIG_AMP)
phase -= 90;
calc_beamsplitter_signal_vars(bs, &factor, &ph1, f_sig->f, f_sig->carrier->f);
// TODO check refractive index
_factor *= 2.0 * (1.0 + f_sig->carrier->f / inter.f0) * sqrt(bs->R);
_ph1 = bs->phi * ((1.0 + f_sig->f / inter.f0) + (1.0 + f_sig->carrier->f / inter.f0));
if(!node1->gnd_node && !node2->gnd_node){
double factor = _factor * cos(bs->alpha_1 * RAD);
double ph1 = _ph1 * cos(bs->alpha_1 * RAD);
double phase1 = phase + ph1;
// get the index for the signal field of frequency f and mode i and
......@@ -818,6 +815,8 @@ void fill_bs_signal_rhs(signal_t *signal, ifo_matrix_vars_t* M_car, ifo_matrix_v
}
if(!node3->gnd_node && !node4->gnd_node){
double factor = _factor * cos(bs->alpha_2 * RAD);
double ph1 = _ph1 * cos(bs->alpha_2 * RAD);
double phase2 = phase - ph1 + 180.0 * f_sig->order;
// get the index for the signal field of frequency f and mode i and
......@@ -877,13 +876,7 @@ void fill_bs_signal_rhs(signal_t *signal, ifo_matrix_vars_t* M_car, ifo_matrix_v
double phase = f_sig->order * signal->phase;
double factor = signal->amplitude * EPSILON * 0.5; //Epsilon/2 for fsig
double turn = 1.0;
// additional scaling from coordinate system misalignments at non-normal incidence
// see section 4.11 (Misalignment angles at a beam splitter) in manual.
// This only affects pitch signals
if (signal->type == SIG_Y)
factor *= cos(bs->alpha * RAD);
complex_t txs1 = complex_0;
complex_t txs2 = complex_0;
complex_t txs3 = complex_0;
......@@ -1069,31 +1062,44 @@ void fill_bs_signal_rhs(signal_t *signal, ifo_matrix_vars_t* M_car, ifo_matrix_v
}
}
}
double AoI_side_1 = 1;
double AoI_side_2 = 1;
// additional scaling from coordinate system misalignments at non-normal incidence
// see section 4.11 (Misalignment angles at a beam splitter) in manual.
// This only affects pitch signals
if (signal->type == SIG_Y){
// the angle of incidence is different on each side depending on
// the refractive indices
AoI_side_1 *= cos(bs->alpha_1 * RAD);
AoI_side_2 *= cos(bs->alpha_2 * RAD);
}
if (!node1->gnd_node) {
sidx = M_sig->node_rhs_idx_1[node1->list_index] + get_node_rhs_idx(2, f, n, M_sig->num_frequencies);
srhs[sidx] = z_by_xph(a21, turn * factor, phase);
srhs[sidx] = z_by_xph(a21, turn * factor * AoI_side_1, phase);
// check if we need to conjugate the value for the lower sideband
if(f_sig->order == -1) srhs[sidx].im *= -1;
}
if (!node2->gnd_node) {
sidx = M_sig->node_rhs_idx_1[node2->list_index] + get_node_rhs_idx(2, f, n, M_sig->num_frequencies);
srhs[sidx] = z_by_xph(a12, turn * factor, phase);
srhs[sidx] = z_by_xph(a12, turn * factor * AoI_side_1, phase);
// check if we need to conjugate the value for the lower sideband
if(f_sig->order == -1) srhs[sidx].im *= -1;
}
if (!node3->gnd_node) {
sidx = M_sig->node_rhs_idx_1[node3->list_index] + get_node_rhs_idx(2, f, n, M_sig->num_frequencies);
srhs[sidx] = z_by_xph(a43, turn * factor, phase + ph2);
srhs[sidx] = z_by_xph(a43, turn * factor * AoI_side_2, phase + ph2);
// check if we need to conjugate the value for the lower sideband
if(f_sig->order == -1) srhs[sidx].im *= -1;
}
if (!node4->gnd_node) {
sidx = M_sig->node_rhs_idx_1[node4->list_index] + get_node_rhs_idx(2, f, n, M_sig->num_frequencies);
srhs[sidx] = z_by_xph(a34, turn * factor, phase + ph2);
srhs[sidx] = z_by_xph(a34, turn * factor * AoI_side_2, phase + ph2);
// check if we need to conjugate the value for the lower sideband
if(f_sig->order == -1) srhs[sidx].im *= -1;
}
......
......@@ -820,7 +820,7 @@ void dump_beamsplitter(FILE *fp, beamsplitter_t *bs) {
fprintf(fp, "angle (ybeta) : %.15g\n", bs->beta_y);
}
fprintf(fp, "incidence (alpha) : %.15g\n", bs->alpha);
fprintf(fp, "incidence (alpha) : %.15g\n", bs->alpha_1);
fprintf(fp, "node1 : %d\n", bs->node1_index);
fprintf(fp, "node2 : %d\n", bs->node2_index);
fprintf(fp, "node3 : %d\n", bs->node3_index);
......@@ -1341,7 +1341,7 @@ void dump_all_beamsplitters(FILE *fp) {
node_print(inter.bs_list[i].node2_index),
node_print(inter.bs_list[i].node3_index),
node_print(inter.bs_list[i].node4_index));
fprintf(fp, " alpha= %.15g", inter.bs_list[i].alpha);
fprintf(fp, " alpha= %.15g", inter.bs_list[i].alpha_1);
if (inter.bs_list[i].attribs == 0) {
fprintf(fp, "\n");
......
......@@ -310,8 +310,8 @@ bool compute_bs_bayer_helms_knm(beamsplitter_t *bs, double nr1, double nr2) {
// angle coefficients that scale depending on the beam rotation
double co1, co2;
co1 = cos(bs->alpha * RAD);
co2 = cos(asin(sin(bs->alpha * RAD) * (nr1 / nr2)));
co1 = cos(bs->alpha_1 * RAD);
co2 = cos(asin(sin(bs->alpha_1 * RAD) * (nr1 / nr2)));
// now loop through each component, check whether we should calculate
// the knm, then interpolate etc.
......
......@@ -540,7 +540,7 @@ void get_bs_int_limit(bs_knm_map_int_params_t *p, int knum) {
*usermessages = *usermessages | USING_POLAR;
}
double ca = 1./cos(p->bs->alpha * RAD);
double ca = 1./cos(p->bs->alpha_1 * RAD);
if (p->usingMap && !(*usermessages & MAP_TOO_SMALL)) {
surface_merged_map_t *map = p->merged_map;
......@@ -580,7 +580,7 @@ void get_bs_int_limit(bs_knm_map_int_params_t *p, int knum) {
, p->xmax[1] - p->xmin[1]
, (map->cols - 1) * map->xstep
, (map->rows - 1) * map->ystep
, p->bs->alpha);
, p->bs->alpha_1);
*usermessages = *usermessages | MAP_TOO_SMALL;
......@@ -701,7 +701,7 @@ void compute_bs_knm_integral(beamsplitter_t *bs, double nr1, double nr2, int int
recompute = recompute | (map->_nr2 != nr2);
recompute = recompute | (map->_angle != bs->map_rotation);
recompute = recompute | (map->_AoI != bs->alpha);
recompute = recompute | (map->_AoI != bs->alpha_1);
recompute = recompute | (map->_x_off != bs->x_off);
recompute = recompute | (map->_y_off != bs->x_off);
......@@ -1139,7 +1139,7 @@ void compute_bs_knm_integral(beamsplitter_t *bs, double nr1, double nr2, int int
map->bs_knm_q.qyt2_42 = bs->knm_q.qyt2_42;
map->_angle = bs->map_rotation;
map->_AoI = bs->alpha;
map->_AoI = bs->alpha_1;
map->_r_aperture = bs->r_aperture;
map->_xbeta = bs->beta_x;
map->_ybeta = bs->beta_y;
......
......@@ -515,14 +515,14 @@ void integrand_get_variables(void *p, KNM_COMPONENT_t knmcmp,
void calc_bs_AOI(bs_knm_map_int_params_t *pbs, KNM_BS_NODE_DIRECTION_t k, double *cos_alpha){
if(pbs->bs->alpha != 0.0) {
if(pbs->bs->alpha_1 != 0.0) {
// if we are computing a reflection from the n1/n2 side we are in nr1
// and the angle of incidence is the same as defined. For the side n3/n4
// though the angle of incidence is defined by Snells law and the refractive
// index difference
if(eq(pbs->nr1,pbs->nr2)){
*cos_alpha = cos(pbs->bs->alpha*RAD);
*cos_alpha = cos(pbs->bs->alpha_1*RAD);
} else {
switch ((KNM_BS_NODE_DIRECTION_t) k) {
......@@ -531,14 +531,14 @@ void calc_bs_AOI(bs_knm_map_int_params_t *pbs, KNM_BS_NODE_DIRECTION_t k, double
case BS21:
case BS31:
case BS42:
*cos_alpha = cos(pbs->bs->alpha*RAD);
*cos_alpha = cos(pbs->bs->alpha_1*RAD);
break;
// these are the directions that end up on the n3/n4 side
case BS13:
case BS24:
case BS34:
case BS43:
*cos_alpha = cos(asin(sin(pbs->bs->alpha*RAD) * pbs->nr1 / pbs->nr2));
*cos_alpha = cos(asin(sin(pbs->bs->alpha_1*RAD) * pbs->nr1 / pbs->nr2));
break;
}
}
......@@ -633,7 +633,7 @@ static int integrand_cuba_serial(const int *ndim, const double xx[], const int *
// When at a non-normal incidence we also need to increase the depth of a
// phase map. As the larger angle of incidence causes a longer path length
// difference
if (pbs->bs->alpha != 0) {
if (pbs->bs->alpha_1 != 0) {
phi *= cos_alpha;
}
} else
......@@ -752,7 +752,7 @@ static int integrand_cuba_para(const int *ndim, const double xx[], const int *nc
// When at a non-normal incidence we also need to increase the depth of a
// phase map. As the larger angle of incidence causes a longer path length
// difference
if (pbs->bs->alpha != 0) {
if (pbs->bs->alpha_1 != 0) {
phi *= cos_alpha;
}
} else
......@@ -1387,7 +1387,7 @@ void do_newton_cotes_int(void *userdata, KNM_COMPONENT_t knmcmp, complex_t *resu
// The scaling for beamsplitters went tilted has not been implemented yet
// for the Newton-Cotes integration method. The beam or map will need
// scaling as is done in the function do_riemann_sum_new_int, for example.
if (pbs->bs->alpha != 0) {
if (pbs->bs->alpha_1 != 0) {
warn("Newton-Cotes solver does not support non-normal beamsplitters with maps");
}
......@@ -1589,7 +1589,7 @@ void do_riemann_sum_new_int(void *userdata, KNM_COMPONENT_t knmcmp, complex_t *r
map_rotation = pmr->mirror->angle;
} else if (knmcmp == BEAMSPLITTER_CMP) {
map_rotation = pbs->bs->map_rotation;
incident_angle = pbs->bs->alpha;
incident_angle = pbs->bs->alpha_1;
cos_alpha = cos(incident_angle * RAD);
} else
gerror("knmcmp parameter not recognised");
......
......@@ -713,7 +713,7 @@ void check_beamsplitter(beamsplitter_t *bs) {
gerror("%s reflectivity must be 0 <= r <= 1\n", bs->name);
}
if (bs->alpha >= 90 || bs->alpha <= -90) {
if (bs->alpha_1 >= 90 || bs->alpha_1 <= -90) {
gerror("%s alpha must be: -90 <= alpha <= 90\n", bs->name);
} else if (bs->T > 1 || bs->T < 0) {
gerror("%s transmittivity must be: 0 <= t <= 1\n", bs->name);
......
......@@ -1895,7 +1895,7 @@ void read_beamsplitter(const char *command_string, int mode) {
bug_error("incorrect mode");
}
bs->phi = tuning_value;
bs->alpha = incidence_angle; // angle of incidence
bs->alpha_1 = incidence_angle; // angle of incidence
// grab the node indices
bs->node1_index = update_node_index_for(node1_name, IN_OUT);
......@@ -2015,7 +2015,7 @@ void read_beamsplitter2(const char *command_string) {
if (atod(incidence_angle_string, &incidence_angle)) {
gerror("Line `%s':\nUnable to read incidence angle\n");
} else {
bs->alpha = incidence_angle;
bs->alpha_1 = incidence_angle;
}
// grab the node indices
......@@ -10362,7 +10362,7 @@ int assign_beamsplitter_parameter(component_param_t *component_param) {
*component_param->min = 0;
} else if (strcasecmp("alpha", component_param->parameter) == 0 ||
strcasecmp("Alpha", component_param->parameter) == 0) {
*component_param->xparam = &(beamsplitter->alpha);
*component_param->xparam = &(beamsplitter->alpha_1);
strcpy(component_param->unit, " [deg]");
beamsplitter->rebuild = 2;
inter.rebuild |= 2;
......
......@@ -1611,7 +1611,7 @@ void fill_optic_motion_coupling(){
// we need to store the mirror motion in units of something other than meters
// using the init.x_scale
double alpha = cos(bs->alpha*RAD);
double alpha = cos(bs->alpha_1*RAD);
// factor for motion to field coupling for upper sideband, lower one is just the conjugate of this
complex_t factor_x_a = co(0, - alpha / (init.clight) * (TWOPI*(inter.f0 + fcar->f)*init.x_scale));
......@@ -2107,7 +2107,7 @@ void fill_matrix_beamsplitter_elements(beamsplitter_t *bs, double fin, double fo
double rb = Jk*sqrt(bs->backscatter);
double r = Jk*sqrt(bs->R - bs->backscatter);
double t = sqrt(bs->T);
double alpha = bs->alpha * RAD;
double alpha = bs->alpha_1 * RAD;
double phi = bs->phi * RAD * (1.0 + fin / inter.f0 + 1.0 + fout / inter.f0) * cos(alpha);
complex_t _r = z_by_xphr(pow_complex_i(K), r, phi + K*bs->dither_phase * RAD);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment