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

#ref 207: Adding in backscatter at a beamsplitter

parent c6621dd4
......@@ -1043,6 +1043,8 @@ typedef struct beamsplitter {
double R; //!< beam splitter power reflectance
double T; //!< beam splitter power transmittance
double backscatter;
double phi; //!< beam splitter tuning \see mirror_t::phi
double x_off;
double y_off;
......@@ -1076,14 +1078,18 @@ typedef struct beamsplitter {
//! matrix element of the component in the interferometer matrix
complex_t ***a12, ***a21, ***a13, ***a31;
complex_t ***a24, ***a42, ***a34, ***a43;
complex_t ***a11, ***a22, ***a33, ***a44;
complex_t *****a12f, *****a21f, *****a13f, *****a31f;
complex_t *****a24f, *****a42f, *****a34f, *****a43f;
complex_t *****a11f, *****a22f, *****a33f, *****a44f;
coupling_info_t a_cplng_12, a_cplng_21;
coupling_info_t a_cplng_13, a_cplng_31;
coupling_info_t a_cplng_24, a_cplng_42;
coupling_info_t a_cplng_34, a_cplng_43;
coupling_info_t a_cplng_11, a_cplng_22;
coupling_info_t a_cplng_33, a_cplng_44;
int x_rhs_idx; // RHS vector index for the mirror translational motion
......
......@@ -676,14 +676,16 @@ void set_coupling_info() {
for(i=0; i<inter.num_beamsplitters; i++){
beamsplitter_t *bs = &inter.bs_list[i];
coupling_info_t *v[8] = {&bs->a_cplng_12, &bs->a_cplng_21,
coupling_info_t *v[12] = {&bs->a_cplng_12, &bs->a_cplng_21,
&bs->a_cplng_34, &bs->a_cplng_43,
&bs->a_cplng_13, &bs->a_cplng_31,
&bs->a_cplng_24, &bs->a_cplng_42};
&bs->a_cplng_24, &bs->a_cplng_42,
&bs->a_cplng_11, &bs->a_cplng_22,
&bs->a_cplng_33, &bs->a_cplng_44};
int j;
for(j=0; j<8; j++){
for(j=0; j<12; j++){
if(!options.use_coupling_reduction) {
__set_no_reduce_coupling_info(v[j]);
} else {
......
......@@ -294,6 +294,7 @@
#define MECH_RX_TF 524288
#define MECH_RY_TF 1048576
#define HOM_ANGLE 2097152
#define BACKSCATTER 4194304
/* phase amplitude map atributes */
#define NO_MAP_SET 0 //!< Default map value
......
......@@ -1827,6 +1827,8 @@ void get_beamsplitter_elements(beamsplitter_t *bs, ifo_matrix_vars_t *matrix, in
int bs_rhs_idx = matrix->bs_rhs_idx[bs->comp_index];
int col = bs_rhs_idx;
bool isBackscattering = is_var_tuned(&bs->backscatter) || bs->backscatter > 0;
if(!node1->gnd_node) {
// port a1
for(fin=0; fin <matrix->num_frequencies; fin++){
......@@ -1836,9 +1838,28 @@ void get_beamsplitter_elements(beamsplitter_t *bs, ifo_matrix_vars_t *matrix, in
// set the value for the diagonal
set_diagonal_element(matrix, bs_rhs_idx, &col, &got_nnz, &bs->nnz_count, n1_port, fin, i, mode);
// Back scatter 1->1
if(isBackscattering) {
for(fout=0; fout < matrix->num_frequencies; fout++){
if((bs->dither_order == 0 && fin==fout) ^ (bs->dither_order !=0 && matrix->bs_f_couple_allocd[bs->dither_list_index][fin][fout])){
if(bs->a_cplng_11.coupling_off){
set_element(matrix, &(bs->a11f[f_offset+fin][f_offset+fout][i][i]), bs_rhs_idx, &col, &got_nnz, &bs->nnz_count, n1_port+1, fout, i, mode);
} else {
for(j=0; j < inter.num_fields; j++){
int nj, mj;
get_tem_modes_from_field_index(&nj, &mj, j);
bool couples = i == j || include_mode_coupling(&bs->a_cplng_11, ni, nj, mi, mj);
if(couples)
set_element(matrix, &(bs->a11f[f_offset+fin][f_offset+fout][i][j]), bs_rhs_idx, &col, &got_nnz, &bs->nnz_count, n1_port+1, fout, j, mode);
}
}
}
}
}
// couplings from a1 -> a4
if(!node2->gnd_node){
if(!node2->gnd_node) {
for(fout=0; fout < matrix->num_frequencies; fout++){
if((bs->dither_order == 0 && fin==fout) ^ (bs->dither_order !=0 && matrix->bs_f_couple_allocd[bs->dither_list_index][fin][fout])){
if(bs->a_cplng_12.coupling_off){
......@@ -1927,6 +1948,25 @@ void get_beamsplitter_elements(beamsplitter_t *bs, ifo_matrix_vars_t *matrix, in
set_diagonal_element(matrix, bs_rhs_idx, &col, &got_nnz, &bs->nnz_count, n2_port, fin, i, mode);
// Back scatter 2->2
if(isBackscattering) {
for(fout=0; fout < matrix->num_frequencies; fout++){
if((bs->dither_order == 0 && fin==fout) ^ (bs->dither_order !=0 && matrix->bs_f_couple_allocd[bs->dither_list_index][fin][fout])){
if(bs->a_cplng_22.coupling_off){
set_element(matrix, &(bs->a22f[f_offset+fin][f_offset+fout][i][i]), bs_rhs_idx, &col, &got_nnz, &bs->nnz_count, n2_port+1, fout, i, mode);
} else {
for(j=0; j < inter.num_fields; j++){
int nj, mj;
get_tem_modes_from_field_index(&nj, &mj, j);
bool couples = i == j || include_mode_coupling(&bs->a_cplng_22, ni, nj, mi, mj);
if(couples)
set_element(matrix, &(bs->a22f[f_offset+fin][f_offset+fout][i][j]), bs_rhs_idx, &col, &got_nnz, &bs->nnz_count, n2_port+1, fout, j, mode);
}
}
}
}
}
if(!node4->gnd_node){
for(fout=0; fout < matrix->num_frequencies; fout++){
if((bs->dither_order == 0 && fin==fout) ^ (bs->dither_order != 0 && matrix->bs_f_couple_allocd[bs->dither_list_index][fin][fout])){
......@@ -1995,6 +2035,25 @@ void get_beamsplitter_elements(beamsplitter_t *bs, ifo_matrix_vars_t *matrix, in
set_diagonal_element(matrix, bs_rhs_idx, &col, &got_nnz, &bs->nnz_count, n3_port, fin, i, mode);
// Back scatter 3->3
if(isBackscattering) {
for(fout=0; fout < matrix->num_frequencies; fout++){
if((bs->dither_order == 0 && fin==fout) ^ (bs->dither_order !=0 && matrix->bs_f_couple_allocd[bs->dither_list_index][fin][fout])){
if(bs->a_cplng_33.coupling_off){
set_element(matrix, &(bs->a33f[f_offset+fin][f_offset+fout][i][i]), bs_rhs_idx, &col, &got_nnz, &bs->nnz_count, n3_port+1, fout, i, mode);
} else {
for(j=0; j < inter.num_fields; j++){
int nj, mj;
get_tem_modes_from_field_index(&nj, &mj, j);
bool couples = i == j || include_mode_coupling(&bs->a_cplng_33, ni, nj, mi, mj);
if(couples)
set_element(matrix, &(bs->a33f[f_offset+fin][f_offset+fout][i][j]), bs_rhs_idx, &col, &got_nnz, &bs->nnz_count, n3_port+1, fout, j, mode);
}
}
}
}
}
if(!node4->gnd_node){
for(fout=0; fout < matrix->num_frequencies; fout++){
if((bs->dither_order == 0 && fin==fout) ^ (bs->dither_order != 0 && matrix->bs_f_couple_allocd[bs->dither_list_index][fin][fout])){
......@@ -2084,6 +2143,25 @@ void get_beamsplitter_elements(beamsplitter_t *bs, ifo_matrix_vars_t *matrix, in
// set the value for the diagonal
set_diagonal_element(matrix, bs_rhs_idx, &col, &got_nnz, &bs->nnz_count, n4_port, fin, i, mode);
// Back scatter 4->4
if(isBackscattering) {
for(fout=0; fout < matrix->num_frequencies; fout++){
if((bs->dither_order == 0 && fin==fout) ^ (bs->dither_order !=0 && matrix->bs_f_couple_allocd[bs->dither_list_index][fin][fout])){
if(bs->a_cplng_44.coupling_off){
set_element(matrix, &(bs->a44f[f_offset+fin][f_offset+fout][i][i]), bs_rhs_idx, &col, &got_nnz, &bs->nnz_count, n4_port+1, fout, i, mode);
} else {
for(j=0; j < inter.num_fields; j++){
int nj, mj;
get_tem_modes_from_field_index(&nj, &mj, j);
bool couples = i == j || include_mode_coupling(&bs->a_cplng_44, ni, nj, mi, mj);
if(couples)
set_element(matrix, &(bs->a44f[f_offset+fin][f_offset+fout][i][j]), bs_rhs_idx, &col, &got_nnz, &bs->nnz_count, n4_port+1, fout, j, mode);
}
}
}
}
}
if(has_optic_dof){
for(j=0; j<bs->num_motions; j++)
set_element_mech(matrix, &bs->a4i_x[j][fin][i], bs->x_rhs_idx+j, &col, &got_nnz, &bs->nnz_count, mode);
......@@ -4261,6 +4339,10 @@ void fill_ccs_matrix(ifo_matrix_vars_t *matrix){
bs->a31 = bs->a31f[f_offset+fin][f_offset+fin];
bs->a24 = bs->a24f[f_offset+fin][f_offset+fin];
bs->a42 = bs->a42f[f_offset+fin][f_offset+fin];
bs->a11 = bs->a11f[f_offset+fin][f_offset+fin];
bs->a22 = bs->a22f[f_offset+fin][f_offset+fin];
bs->a33 = bs->a33f[f_offset+fin][f_offset+fin];
bs->a44 = bs->a44f[f_offset+fin][f_offset+fin];
fill_matrix_beamsplitter_elements(bs, f_list[fin]->f, f_list[fin]->f, 0, do_conj);
} else {
......@@ -4274,7 +4356,11 @@ void fill_ccs_matrix(ifo_matrix_vars_t *matrix){
bs->a31 = bs->a31f[f_offset+fin][f_offset+fout];
bs->a24 = bs->a24f[f_offset+fin][f_offset+fout];
bs->a42 = bs->a42f[f_offset+fin][f_offset+fout];
bs->a11 = bs->a11f[f_offset+fin][f_offset+fout];
bs->a22 = bs->a22f[f_offset+fin][f_offset+fout];
bs->a33 = bs->a33f[f_offset+fin][f_offset+fout];
bs->a44 = bs->a44f[f_offset+fin][f_offset+fout];
fill_matrix_beamsplitter_elements(bs, f_list[fin]->f, f_list[fout]->f, matrix->bs_f_couple_order[bs->dither_list_index][fin][fout], do_conj);
}
}
......
......@@ -1759,6 +1759,19 @@ int allocate_memory_for_beamsplitter_list(long int *bytes) {
beamsplitter_t *bs = &inter.bs_list[i];
int err;
err = alloc_field_freq_ptrs_mod(&bs->a11f, bytes, 7, 7, 8);
if(err) return err;
err = alloc_field_freq_ptrs_mod(&bs->a22f, bytes, 7, 7, 8);
if(err) return err;
err = alloc_field_freq_ptrs_mod(&bs->a33f, bytes, 7, 7, 8);
if(err) return err;
err = alloc_field_freq_ptrs_mod(&bs->a44f, bytes, 7, 7, 8);
if(err) return err;
err = alloc_field_freq_ptrs_mod(&bs->a12f, bytes, 7, 7, 8);
if(err) return err;
......
......@@ -7320,6 +7320,8 @@ int get_attribute_type(const char *parameter_name,
attribute_type = MECH_RY_TF;
} else if (strncasecmp("homangle", parameter_name, 8) == 0) {
attribute_type = HOM_ANGLE;
} else if (strncasecmp("backscatter", parameter_name, 11) == 0) {
attribute_type = BACKSCATTER;
} else {
gerror("Line `%s':\nno attribute named '%s', "
"use Rap, map_deg, M, Rc, Rcx, Rcy, gx, gy, xbeta, ybeta, alpha, eta or rho!\n",
......@@ -7878,6 +7880,15 @@ void set_beamsplitter_attribute(component_attribute_t *component_attribute) {
double parameter_value = component_attribute->parameter_value;
switch (attribute_type) {
case BACKSCATTER:
if (parameter_value <= 0 && parameter_value < 1) {
gerror("Line `%s':\nbackscatter must be 0 > backscatter < 1\n",
component_attribute->command_string);
}
beamsplitter->backscatter = parameter_value;
break;
case MASS:
if (parameter_value <= 0) {
gerror("Line `%s':\nmass must be >0\n",
......@@ -10250,6 +10261,12 @@ int assign_beamsplitter_parameter(component_param_t *component_param) {
strcpy(component_param->unit, " [deg]");
*component_param->lborder = 0;
*component_param->uborder = 0;
} else if (strcasecmp("backscatter", component_param->parameter) == 0) {
*component_param->xparam = &(beamsplitter->backscatter);
strcpy(component_param->unit, " ");
*component_param->lborder = 0;
*component_param->uborder = 0;
} else {
server_gerror("Line `%s':\nno parameter '%s' at beamsplitter\n",
......
......@@ -2105,7 +2105,8 @@ void fill_matrix_beamsplitter_elements(beamsplitter_t *bs, double fin, double fo
else
Jk = bessn(K, TWOPI * bs->dither_m / init.lambda * (1.0 + fin / inter.f0 + 1.0 + fout / inter.f0));
double r = Jk*sqrt(bs->R);
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 phi = bs->phi * RAD * (1.0 + fin / inter.f0 + 1.0 + fout / inter.f0) * cos(alpha);
......@@ -2115,12 +2116,27 @@ void fill_matrix_beamsplitter_elements(beamsplitter_t *bs, double fin, double fo
double c0 = _r.re;
double s0 = _r.im;
// term for backscatter
complex_t _rb = z_by_xphr(pow_complex_i(K), rb, phi + K*bs->dither_phase * RAD);
// field amplitudes
int j, k;
for (j = 0; j < inter.num_fields; j++) { // outputs
for (k = 0; k < inter.num_fields; k++) { // frequencies
if (bs->a11[j][k] != NULL && bs->backscatter > 0) {
// Backscatter
bs->a11[j][k]->re = -_rb.re;
bs->a11[j][k]->im = -_rb.im;
// TODO !!! Coupling matrix is not consider yet for backscatter
if(conjugate){
bs->a11[j][k]->im *= -1;
}
}
if (bs->a12[j][k] != NULL) {
(bs->a12[j][k])->re = -c0;
(bs->a12[j][k])->im = -s0;
......@@ -2135,6 +2151,18 @@ void fill_matrix_beamsplitter_elements(beamsplitter_t *bs, double fin, double fo
//printf("after a12, a21 (%d, %d) %s %s\n",j,k,complex_form(*(bs->a12[j][k])), complex_form(*(bs->a21[j][k])));
}
if (bs->a22[j][k] != NULL && bs->backscatter > 0) {
// Backscatter
bs->a22[j][k]->re = -_rb.re;
bs->a22[j][k]->im = -_rb.im;
// TODO !!! Coupling matrix is not consider yet for backscatter
if(conjugate){
bs->a22[j][k]->im *= -1;
}
}
if (bs->a21[j][k] != NULL) {
(bs->a21[j][k])->re = -c0;
(bs->a21[j][k])->im = -s0;
......@@ -2173,6 +2201,18 @@ void fill_matrix_beamsplitter_elements(beamsplitter_t *bs, double fin, double fo
//printf("a13, a31 (%d, %d) %s %s\n",j,k,complex_form(*(bs->a13[j][k])), complex_form(*(bs->a31[j][k])));
}
if (bs->a33[j][k] != NULL && bs->backscatter > 0) {
// Backscatter
bs->a33[j][k]->re = -_rb.re;
bs->a33[j][k]->im = _rb.im;
// TODO !!! Coupling matrix is not consider yet for backscatter
if(conjugate){
bs->a33[j][k]->im *= -1;
}
}
if (bs->a34[j][k] != NULL) {
(bs->a34[j][k])->re = -c0;
(bs->a34[j][k])->im = s0;
......@@ -2186,6 +2226,18 @@ void fill_matrix_beamsplitter_elements(beamsplitter_t *bs, double fin, double fo
//printf("after a34, a43 (%d, %d) %s %s\n",j,k,complex_form(*(bs->a34[j][k])), complex_form(*(bs->a43[j][k])));
}
if (bs->a44[j][k] != NULL && bs->backscatter > 0) {
// Backscatter
bs->a44[j][k]->re = -_rb.re;
bs->a44[j][k]->im = _rb.im;
// TODO !!! Coupling matrix is not consider yet for backscatter
if(conjugate){
bs->a44[j][k]->im *= -1;
}
}
if (bs->a43[j][k] != NULL) {
(bs->a43[j][k])->re = -c0;
(bs->a43[j][k])->im = s0;
......
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