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

trying to fix PI hom coupling

parent 66ec0b24
......@@ -1330,8 +1330,8 @@ int set_k_mirror(int mirror_index) {
mirror->knm.k22[n][m] = rev_gouy(mirror->knm_no_rgouy.k22[n][m], n1, m1, n2, m2, knm_q.qxt1_22, knm_q.qxt2_22, knm_q.qyt1_22, knm_q.qyt2_22);
for(k=0; k<mirror->num_surface_motions; k++){
mirror->knm_surf_x_a_1[k][n][m] = rev_gouy(mirror->knm_surf_x_a_1[k][n][m], n1, m1, n2, m2, knm_q.qxt1_11, knm_q.qxt2_11, knm_q.qyt1_11, knm_q.qyt2_11);
mirror->knm_surf_x_a_2[k][n][m] = rev_gouy(mirror->knm_surf_x_a_2[k][n][m], n1, m1, n2, m2, knm_q.qxt1_22, knm_q.qxt2_22, knm_q.qyt1_22, knm_q.qyt2_22);
//mirror->knm_surf_x_a_1[k][n][m] = rev_gouy(mirror->knm_surf_x_a_1[k][n][m], n1, m1, n2, m2, knm_q.qxt1_11, knm_q.qxt2_11, knm_q.qyt1_11, knm_q.qyt2_11);
//mirror->knm_surf_x_a_2[k][n][m] = rev_gouy(mirror->knm_surf_x_a_2[k][n][m], n1, m1, n2, m2, knm_q.qxt1_22, knm_q.qxt2_22, knm_q.qyt1_22, knm_q.qyt2_22);
}
}
}
......@@ -1369,6 +1369,7 @@ int set_k_mirror(int mirror_index) {
for (n = 0; n < inter.num_fields; n++) {
for (m = 0; m < inter.num_fields; m++) {
get_tem_modes_from_field_index(&n1, &m1, n);
get_tem_modes_from_field_index(&n2, &m2, m);
if (inter.set_tem_phase_zero & 1) {
......@@ -1393,6 +1394,10 @@ int set_k_mirror(int mirror_index) {
mirror->knm_surf_x_a_2[k][n][m] = cminus(mirror->knm_surf_x_a_2[k][n][m]);
}
}
// for(k=0; k<mirror->num_surface_motions; k++){
// warn("knmnm %i%i->%i%i = %s\n", n1, m1, n2, m2, complex_form15(mirror->knm_surf_x_a_1[k][n][m]));
// }
}
}
......
......@@ -1139,6 +1139,11 @@ void fill_light_signal_rhs(signal_t *signal, ifo_matrix_vars_t* M_car, ifo_matri
}
for (j = 0; j < inter.num_fields; j++) {
// the factor of 2 here comes from the fact we only model
// the positive frequency components. This wasn't a problem until
// radiation pressure was added as there were no effects where effects
// were proportional to the power. With RP though we need to use the
// proper scaling.
complex_t input_amplitude = z_by_ph(co(sqrt(2*light_input->I0 / init.epsilon_0_c), 0.0), light_input->phase);
input_amplitude = z_by_z(input_amplitude, light_input->power_coeff_list[j]);
......@@ -1776,6 +1781,12 @@ void data_point_new(){
bug_error("Node did not have direction set when filling in RHS vector");
for (j = 0; j < inter.num_fields; j++) {
// the factor of 2 here comes from the fact we only model
// the positive frequency components. This wasn't a problem until
// radiation pressure was added as there were no effects where effects
// were proportional to the power. With RP though we need to use the
// proper scaling.
// !! Must match up with fill_light_signal_rhs too!!
complex_t input_amplitude = z_by_ph(co(sqrt(2*light_input->I0 / init.epsilon_0_c), 0.0), light_input->phase);
input_amplitude = z_by_z(input_amplitude, light_input->power_coeff_list[j]);
......@@ -3374,11 +3385,12 @@ void minimize(){
if(inter.minimizer == NULL) return;
double m = *inter.minimizer->variable;
double a = inter.minimizer->A;
double b = inter.minimizer->B;
double a = m+inter.minimizer->A;
double b = m+inter.minimizer->B;
local_min(a, b, inter.minimizer->abserr, 1e-15, &fn_minimizer, &m);
//glomin (a, b, m, 10, 2e-16, inter.minimizer->abserr, inter.minimizer->abserr, &fn_minimizer, &m);
inter.minimizer->x_result = m;
inter.minimizer->result_output->signal.re = m;
......
......@@ -2257,7 +2257,7 @@ void print_help() {
" gr[n] name d node1 node2 [node3 [node4]] - grating\n"
" isol name S node1 node2 [node3] - isolator\n"
" mod name f midx order am/pm [phase] node1 node2 - modulator\n"
" lens f node1 node2 - thin lens\n"
" lens name f node1 node2 - thin lens\n"
" sq name f r angle node - squeezed input\n"
"** Detectors:\n"
" pd[n] name [f1 [phase1 [f2... ]]] node[*] - photodetector [mixer]\n"
......
......@@ -954,7 +954,7 @@ void calc_mirror_knm_surf_motions_map(mirror_t *m, double nr1, double nr2, compl
// the msign(n1+n2) is the difference between the incoming beam
// coordinate system and that of the maps, the x-axis is inverted.
m->knm_surf_motion_1o[k][in][out] = knm;
m->knm_surf_motion_1i[k][in][out] = cconj(z_by_x(knm, msign(p.n1+p.n2)));
m->knm_surf_motion_1i[k][in][out] = knm; //cconj(z_by_x(knm, msign(p.n1+p.n2)));
//warn("%i%i->%i%i = %.15g\n",p.n1, p.m1, p.n2, p.m2, zabs(knm));
......
......@@ -12110,7 +12110,7 @@ void read_debug(const char *command_string) {
void read_minimize(const char *command_string) {
const char *usage = "maximize/minimize name output parameter target parameter lower upper abserr iters";
const char *usage = "maximize/minimize name output parameter target parameter lower upper abserr";
char command_name[LINE_LEN] = {0};
char name[LINE_LEN] = {0};
......@@ -12122,12 +12122,11 @@ void read_minimize(const char *command_string) {
char upper[LINE_LEN] = {0};
char abserr[LINE_LEN] = {0};
char rest_string[LINE_LEN] = {0};
int iters = 1000;
int num_vars_read = sscanf(command_string, "%s %s %s %s %s %s %s %s %s %i %[^\n]",
command_name, name, output_name, output_param, target_name, target_param, lower, upper, abserr, &iters, rest_string);
int num_vars_read = sscanf(command_string, "%s %s %s %s %s %s %s %s %s %[^\n]",
command_name, name, output_name, output_param, target_name, target_param, lower, upper, abserr, rest_string);
if (num_vars_read != 10) {
if (num_vars_read != 9) {
gerror("Line `%s':\nexpected '%s'\n", command_string, usage);
}
......@@ -12196,11 +12195,6 @@ void read_minimize(const char *command_string) {
gerror("Line `%s':\nMinimize error not a valid floating point number\n", command_string);
}
if(iters <= 0) {
gerror("Line `%s':\nNumber of iterations should be > 0\n", command_string);
}
output_data_t *result = &inter.output_data_list[inter.num_outputs];
strcpy(result->name, name);
......
......@@ -1388,7 +1388,9 @@ void fill_optic_motion_coupling(){
complex_t factor_x_a = co(0, - 1.0 / (init.clight) * (TWOPI*(inter.f0 + fcar->f)*init.x_scale));
// factor from field to motion coupling
double factor_a_x = init.epsilon_0_c / (init.clight * init.x_scale);
// extra factor of 2 due to the fact we throw away the
// negative frequency term of the power fluctuations
double factor_a_x = 2 * 0.5*init.epsilon_0_c / (init.clight * init.x_scale);
// TODO add index of refraction to these later
// phase that the carrier and sideband will pick up due to tuning
......@@ -1613,7 +1615,9 @@ void fill_optic_motion_coupling(){
complex_t factor_x_a = co(0, - alpha / (init.clight) * (TWOPI*(inter.f0 + fcar->f)*init.x_scale));
// factor from field to motion coupling
double factor_a_x = alpha * init.epsilon_0_c / (init.clight * init.x_scale);
// extra factor of 2 due to the fact we throw away the
// negative frequency term of the power fluctuations
double factor_a_x = alpha * 2 * 0.5 * init.epsilon_0_c / (init.clight * init.x_scale);
// TODO add index of refraction to these later
// phase that the carrier and sideband will pick up due to tuning
......@@ -2090,13 +2094,14 @@ void dump_matrix_space_elements(FILE *fp) {
void fill_matrix_beamsplitter_elements(beamsplitter_t *bs, double fin, double fout, int K, bool conjugate) {
double Jk;
if(bs->dither_m == 0.0 && K == 0)
if(bs->dither_m == 0.0 && K == 0) {
Jk = 1.0;
else
} else {
if (K < 0)
Jk = pow_neg_1(K) * bessn(-K, TWOPI * bs->dither_m / init.lambda * (1.0 + fin / inter.f0 + 1.0 + fout / inter.f0));
else
Jk = bessn(K, TWOPI * bs->dither_m / init.lambda * (1.0 + fin / inter.f0 + 1.0 + fout / inter.f0));
}
double rb = Jk*sqrt(bs->backscatter);
double r = Jk*sqrt(bs->R - bs->backscatter);
......@@ -2110,7 +2115,7 @@ void fill_matrix_beamsplitter_elements(beamsplitter_t *bs, double fin, double fo
double s0 = _r.im;
// term for backscatter
complex_t _rb = z_by_xphr(pow_complex_i(K), rb, phi + K*bs->dither_phase * RAD);
complex_t _rb = z_by_xphr(pow_complex_i(K), -rb, phi + K*bs->dither_phase * RAD);
// field amplitudes
int j, k;
......
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