Commit 912a04c6 authored by Daniel Brown's avatar Daniel Brown

Adding in code to enable beamsplitters to do radius of curvature fsig...

Adding in code to enable beamsplitters to do radius of curvature fsig modulations. Still a work in progress. I also have made the functions get_node_q_to and get_node_q_from as a convenience function for more visibly getting the q in a particular direction at a node, relative to some component.
parent 856e80f2
......@@ -522,6 +522,8 @@ int get_signal_type(char *signal_type_string) {
return SIG_X;
} else if (string_matches_exactly(signal_type_string, "y")) {
return SIG_Y;
} else if (string_matches_exactly(signal_type_string, "Rc")) {
return SIG_RC;
} else if (string_matches_exactly(signal_type_string, "sagnac")) {
return SIG_SAGNAC;
} else if (signal_type_string[0] == 's' && isdigit(signal_type_string[1])) {
......@@ -2277,6 +2279,57 @@ void set_all_tems(void) {
}
}
/***
* Gets the complex beam parameter at a node going into a component.
* Doesn't do any checks to see if the node is connected to the component.
*
* node: ptr to node struct
* component_idx: unique component index
* qx, qy: ptr to complex_t to store beam params in, if NULL it won't set
*
* Returns 0 if ground/dump node
*/
int get_node_q_to(node_t *node, int component_idx, complex_t *qx, complex_t *qy){
if (NOT node->gnd_node) {
// If q was traced away from a component, a node it will have the component index
// set to it, so here we reverse it
if (node->component_index == component_idx) {
if(qx) *qx = cminus(cconj(node->qx));
if(qy) *qy = cminus(cconj(node->qy));
} else {
if(qx) *qx = node->qx;
if(qy) *qy = node->qy;
}
} else {
return 0;
}
}
/***
* Gets the complex beam parameter at a node going from a component.
* Doesn't do any checks to see if the node is connected to the component.
*
* node: ptr to node struct
* component_idx: unique component index
* qx, qy: ptr to complex_t to store beam params in, if NULL it won't set
*
* Returns 0 if ground/dump node
*/
int get_node_q_from(node_t *node, int component_idx, complex_t *qx, complex_t *qy){
if (NOT node->gnd_node) {
// If q was traced "from" a node it will have the component index
// set to it
if (node->component_index == component_idx) {
if(qx) *qx = node->qx;
if(qy) *qy = node->qy;
} else {
if(qx) *qx = cminus(cconj(node->qx));
if(qy) *qy = cminus(cconj(node->qy));
}
} else {
return 0;
}
}
//! Returns factorial for n<=100
......
......@@ -145,6 +145,9 @@ int get_function_idx(char* name);
int strhash(char *key);
int get_node_q_from(node_t *node, int component_idx, complex_t *qx, complex_t *qy);
int get_node_q_to(node_t *node, int component_idx, complex_t *qx, complex_t *qy);
#if OS == __WIN_32__ || OS == __WIN_64__
const char *strcasestr(const char *s1, const char *s2);
#endif
......
......@@ -632,8 +632,9 @@ void fill_signal_rhs() {
* \param q Gaussian beam q parameter
*
* \todo sig_node -> signal_node
*/
**/
complex_t xyamp(node_t *sig_node, complex_t q) {
complex_t txs;
double nr, w0;
......@@ -693,6 +694,51 @@ complex_t calc_reflect_tilt_signal(signal_t *signal, complex_t **knm, int in, in
return mult;
}
/**
* Computes the scattering matrix element for a small dioptre modulation of a components focusing power.
* Only considered on reflection, but probably also applicable on transmission of a lens. This is
* computed by asuming some perturbing time oscillating C-matrix-term in the ABCD.
*
* Essentially returns a matrix that scatters carrier modes to n+-2 modes of a signal sideband.
* The direct coupling n->n pick up some phase shift, due to the modulating waist position.
*
* Coupling coefficient returned is not reversed gouy'd.
*
* qin: beam basis being modulated
* amplitude, phase: amp and phs of signal, phs in degrees
* P: Focal power of the object being modulated.
*
* returns the coupling term for mode n -> n_
*/
complex_t calc_dioptre_modulation_knn(int n, int n_, complex_t qin, double P,
double amplitude, double phase) {
int dn = n - n_;
complex_t knn_ = complex_0;
if (dn == 2 || dn == -2 || dn == 0){
double zr = qin.im;
double z = qin.re;
double fac = sqrt((n+1) * (n_+2));
double eps_z = zr * (2/(P*P*zr*zr + (P*z-1)*(P*z-1)) - 1) - (z*z)/(zr);
double eps_zr = (2*z*(P*z-1) + 2*P*zr*zr) / (P*P*zr*zr + (P*z-1)*(P*z-1));
complex_t eps_q = co(eps_z, eps_zr);
if (dn == 2)
knn_ = z_by_z( z_by_x(complex_i, -fac / 4.0), eps_q);
else if(dn == -2)
knn_ = z_by_z( z_by_x(complex_i, -fac / 4.0), cconj(eps_q));
else
knn_ = co(1, 0.25 * (2*n+1)*eps_z);
} else {
knn_ = complex_0;
}
return z_by_xph(knn_, amplitude, phase);
}
void fill_bs_signal_rhs(signal_t *signal, ifo_matrix_vars_t* M_car, ifo_matrix_vars_t* M_sig){
assert(M_sig != NULL);
assert(M_car != NULL);
......@@ -705,7 +751,7 @@ void fill_bs_signal_rhs(signal_t *signal, ifo_matrix_vars_t* M_car, ifo_matrix_v
node_t *node2 = &inter.node_list[bs->node2_index];
node_t *node3 = &inter.node_list[bs->node3_index];
node_t *node4 = &inter.node_list[bs->node4_index];
for(i=0; i<bs->num_motions; i++){
// If a signal is applied and we have some motions available we apply
// the signal to the motion term, which then generates the sidebands
......@@ -1106,6 +1152,132 @@ void fill_bs_signal_rhs(signal_t *signal, ifo_matrix_vars_t* M_car, ifo_matrix_v
}
}
}
} else if((signal->type == SIG_RC)) {
for(f=0; f<M_sig->num_frequencies; f++){
frequency_t *f_sig = M_sig->frequencies[f];
double phase = f_sig->order * signal->phase;
double factor = signal->amplitude * EPSILON * 0.5; //Epsilon/2 for fsig
double turn = 1.0;
int n, i, j;
int cidx, sidx;
complex_t *crhs = (complex_t*)M_car->rhs_values;
complex_t *srhs = (complex_t*)M_sig->rhs_values;
// compute purely reflected carrier field, do this by taking the incoming beam and computing matrix multiplication
// with the reflection elements of the matrix. Store reflected field in memory already allocated for this in the
// signal structure
int ni, mi, nj, mj;
for(i=0; i < inter.num_fields; i++){
get_tem_modes_from_field_index(&ni, &mi, i);
signal->ar1[i] = complex_0;
signal->ar2[i] = complex_0;
signal->ar3[i] = complex_0;
signal->ar4[i] = complex_0;
for(j=0; j < inter.num_fields; j++){
get_tem_modes_from_field_index(&nj, &mj, j);
if(!node1->gnd_node && !node2->gnd_node) {
if(include_mode_coupling(&bs->a_cplng_21, ni,nj,mi,mj)){
cidx = M_car->node_rhs_idx_1[node2->list_index] + get_node_rhs_idx(1, f_sig->carrier->index, j, M_car->num_frequencies);
signal->ar1[i] = z_pl_z(signal->ar1[i], z_by_z(*bs->a21f[f_sig->carrier->index][f_sig->carrier->index][j][i], crhs[cidx]));
}
if(include_mode_coupling(&bs->a_cplng_12, ni,nj,mi,mj)){
cidx = M_car->node_rhs_idx_1[node1->list_index] + get_node_rhs_idx(1, f_sig->carrier->index, j, M_car->num_frequencies);
signal->ar2[i] = z_pl_z(signal->ar2[i], z_by_z(*bs->a12f[f_sig->carrier->index][f_sig->carrier->index][j][i], crhs[cidx]));
}
}
if(!node3->gnd_node && !node4->gnd_node) {
if(include_mode_coupling(&bs->a_cplng_43, ni,nj,mi,mj)){
cidx = M_car->node_rhs_idx_1[node4->list_index] + get_node_rhs_idx(1, f_sig->carrier->index, j, M_car->num_frequencies);
signal->ar3[i] = z_pl_z(signal->ar3[i], z_by_z(*bs->a43f[f_sig->carrier->index][f_sig->carrier->index][j][i], crhs[cidx]));
}
if(include_mode_coupling(&bs->a_cplng_34, ni,nj,mi,mj)){
cidx = M_car->node_rhs_idx_1[node3->list_index] + get_node_rhs_idx(1, f_sig->carrier->index, j, M_car->num_frequencies);
signal->ar4[i] = z_pl_z(signal->ar4[i], z_by_z(*bs->a34f[f_sig->carrier->index][f_sig->carrier->index][j][i], crhs[cidx]));
}
}
}
}
complex_t dq_fac = complex_1;
double sig_sb_sign;
if(f_sig->order == -1)
sig_sb_sign = 180;
else
sig_sb_sign = 0;
int ii = 0;
// Different reflections to compute
node_t *nodes[4][2] = {
{node1, node2},
{node2, node1},
{node3, node4},
{node4, node3}
};
for(ii=0; ii<4; ii++){
node_t *node_from = nodes[ii][0];
node_t *node_to = nodes[ii][1];
if (!node_from->gnd_node && !node_to->gnd_node) {
// Once this is done we then modulate the outgoing modes
for (n = 0; n < inter.num_fields; n++) {
int tn, tm;
get_tem_modes_from_field_index(&tn, &tm, n);
int tem_m2 = ((tn - 2) >= 0 && (tn - 2) < inter.num_fields) ? get_field_index_from_tem(tn - 2, tm) : -1;
int tem_p2 = ((tn + 2) < inter.num_fields) ? get_field_index_from_tem(tn + 2, tm) : -1;
double P = 2.0*(*node_from->n)/(bs->Rcx);
complex_t qin = complex_0;
complex_t qout = complex_0;
if(tem_m2 >= 0 && tem_m2 < inter.num_fields) { // check if this is a valid mode
// if so, take this carrier mode and scatter it into a higher order mode
get_node_q_to(node_from, bs->comp_index, &qin, NULL);
get_node_q_from(node_to, bs->comp_index, &qout, NULL);
complex_t k_n_n_2 = calc_dioptre_modulation_knn(tn, tn-2, qin, P, factor, phase);
k_n_n_2 = z_by_phr(k_n_n_2, -1.0 * (((n-2 + 0.5) * gouy(qout) - ((n + 0.5) * gouy(qout)))));
// rhs index for node first port + index offset for second port field
sidx = M_sig->node_rhs_idx_1[node_to->list_index] + get_node_rhs_idx(2, f, get_field_index_from_tem(tn-2, tm), M_sig->num_frequencies);
if(f_sig->order == -1) // lower sideband is conjugated
z_inc_zc(&srhs[sidx], z_by_z(signal->ar2[n], k_n_n_2));
else
z_inc_z (&srhs[sidx], z_by_z(signal->ar2[n], k_n_n_2));
}
if(tem_p2 >= 0 && tem_p2 < inter.num_fields) { // check if this is a valid mode
// if so, take this carrier mode and scatter it into a higher order mode
get_node_q_to(node_from, bs->comp_index, &qin, NULL);
get_node_q_from(node_to, bs->comp_index, &qout, NULL);
complex_t k_n_n2 = calc_dioptre_modulation_knn(tn, tn+2, qin, P, factor, phase);
k_n_n2 = z_by_phr(k_n_n2, -1.0 * (((n+2 + 0.5) * gouy(qout) - ((n + 0.5) * gouy(qout)))));
// rhs index for node first port + index offset for second port field
sidx = M_sig->node_rhs_idx_1[node_to->list_index] + get_node_rhs_idx(2, f, get_field_index_from_tem(tn+2, tm), M_sig->num_frequencies);
if(f_sig->order == -1) // lower sideband is conjugated
z_inc_zc(&srhs[sidx], z_by_z(signal->ar2[n], k_n_n2));
else
z_inc_z (&srhs[sidx], z_by_z(signal->ar2[n], k_n_n2));
}
}
}
}
}
}
}
......
......@@ -362,6 +362,7 @@
#define SIG_Z 11
#define NOSIGNAL 12 /** No signal applied, just a name given to signal frequency */
#define SIG_SAGNAC 13
#define SIG_RC 14
/* phasemap types */
#define PM_REFL 1 //!< Reflection phasemap from file
......
......@@ -6803,7 +6803,7 @@ void init_beamsplitter_signal(
if (type_is_set) {
if (signal->type != SIG_AMP && signal->type != SIG_PHS &&
signal->type != SIG_X && signal->type != SIG_Y && signal->type != SIG_Z
&& signal->type != SIG_FZ && signal->type != SIG_FRY && signal->type != SIG_FRX) {
&& signal->type != SIG_FZ && signal->type != SIG_FRY && signal->type != SIG_FRX && signal->type != SIG_RC) {
gerror("Line `%s':\ninvalid signal type for a beam splitter\n",
command_string);
}
......@@ -6865,8 +6865,9 @@ void init_mirror_signal(
gerror("Line `%s':\ninvalid signal type for a mirror\n", command_string);
}
} else if (signal->type != SIG_AMP && signal->type != SIG_PHS &&
signal->type != SIG_X && signal->type != SIG_Y && signal->type != SIG_Z
&& signal->type != SIG_FZ && signal->type != SIG_FRY && signal->type != SIG_FRX) {
signal->type != SIG_X && signal->type != SIG_Y && signal->type != SIG_Z
&& signal->type != SIG_FZ && signal->type != SIG_FRY && signal->type != SIG_FRX
&& signal->type != SIG_RC) {
gerror("Line `%s':\ninvalid signal type for a mirror\n", command_string);
}
} else {
......
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