Commit 628102d9 authored by Daniel Brown's avatar Daniel Brown
Browse files

adding in knm matrix loss calculation and noise injections for all components

parent b28af276
......@@ -137,6 +137,17 @@ added features:
by editting the GNUPLOT terminal settings in kat.ini with the relevant
command. When using v5.0 make sure you update the `gnuversion'
variable in kat.ini.
o Changed quantum noise injection to take advantage of the loss information
in the coupling matrices. We now inject the noise lost due to couplings
such as mode-mismatch, tits, maps, etc. by assuming all the power that
isn't coupled is replaced with vacuum noise. This means a large number
of modes to get quantum noise converging shouldn't be needed. You only
need enough to get the classical sidebands/carriers converging. This
is assuming the higher order modes not included do not contain any
quantum correllations, if they do your quantum noise values may still
require more modes to converge completely.
bug fixes:
o Fixed bug with scaling of qshotS and qnoisedS commands, also fixed
......
......@@ -318,6 +318,7 @@ typedef struct mirror_knm {
zmatrix k21;
//! coupling coefficient for TEM mode from mode 2 into mode 2
zmatrix k22;
bool IsIdentities;
} mirror_knm_t;
......@@ -819,6 +820,11 @@ typedef struct mirror {
mirror_knm_t knm_bayer_helms;
mirror_knm_t knm_aperture;
double *k12_sqrd_sum; // effective per mode loss due to scattering from node 1 to 2
double *k21_sqrd_sum; // effective per mode loss due to scattering from node 2 to 1
double *k22_sqrd_sum; // effective per mode loss due to scattering from node 2 to 2
double *k11_sqrd_sum; // effective per mode loss due to scattering from node 1 to 1
//! Transformation matrix for q of mirror; reflection in tangential plane (node 1 -> node1)
ABCD_t qqr1t;
//! Transformation matrix for q of mirror; reflection in sagittal plane (node 1 -> node1)
......@@ -992,11 +998,11 @@ typedef struct space {
coupling_info_t a_cplng_12, a_cplng_21;
//! coupling coefficient for TEM mode from mode 1 into mode 2
complex_t **k12;
//! coupling coefficient for TEM mode from mode 2 into mode 1
complex_t **k21;
complex_t **k12; //! coupling coefficient for TEM mode from mode 1 into mode 2
complex_t **k21; //! coupling coefficient for TEM mode from mode 2 into mode 1
double *k12_sqrd_sum; // effective per mode loss due to scattering from node 1 to 2
double *k21_sqrd_sum; // effective per mode loss due to scattering from node 2 to 1
int node1_index; //!< index of node1
int node2_index; //!< index of node2
......@@ -1190,6 +1196,15 @@ typedef struct beamsplitter {
//!< Stores the knm values computed from bayer-helms computations
bs_knm_t knm_bayer_helms;
double *k12_sqrd_sum; // effective per mode loss due to scattering from node 1 to 2
double *k21_sqrd_sum; // effective per mode loss due to scattering from node 2 to 1
double *k34_sqrd_sum; // effective per mode loss due to scattering from node 3 to 4
double *k43_sqrd_sum; // effective per mode loss due to scattering from node 4 to 3
double *k13_sqrd_sum; // effective per mode loss due to scattering from node 1 to 3
double *k31_sqrd_sum; // effective per mode loss due to scattering from node 3 to 1
double *k24_sqrd_sum; // effective per mode loss due to scattering from node 2 to 4
double *k42_sqrd_sum; // effective per mode loss due to scattering from node 4 to 2
//int map_type_set; //!< Which map types have been set
surface_map_t * map[MAX_MAPS]; //!< Surface map storage
int num_maps; //!< number of surface maps
......@@ -1727,6 +1742,9 @@ typedef struct modulator {
//! coupling coefficient for TEM mode from mode 2 into mode 1
complex_t **k21;
double *k12_sqrd_sum;
double *k21_sqrd_sum;
ABCD_t qq; //!< transformation matrix for q of modulator
complex_t ***q12; //!< quantum noise matrix element mode 1 -> mode 2
......@@ -1781,6 +1799,11 @@ typedef struct diode {
complex_t **k32; //! coupling coefficient for TEM mode from node 3 into mode 2
complex_t **k23; //! coupling coefficient for TEM mode from node 2 into mode 3
double *k12_sqrd_sum; // effective per mode loss due to scattering from node 1 to 2
double *k21_sqrd_sum; // effective per mode loss due to scattering from node 2 to 1
double *k32_sqrd_sum; // effective per mode loss due to scattering from node 3 to 2
double *k23_sqrd_sum; // effective per mode loss due to scattering from node 3 to 2
ABCD_t qq; //!< transformation matrix for q of diode
complex_t ***q12; //!< quantum noise matrix element mode 1 -> mode 2
......@@ -1827,7 +1850,10 @@ typedef struct lens {
complex_t **k12;
//! coupling coefficient for TEM mode from mode 2 into mode 1
complex_t **k21;
double *k12_sqrd_sum; // effective per mode loss due to scattering from node 1 to 2
double *k21_sqrd_sum; // effective per mode loss due to scattering from node 2 to 1
ABCD_t qq; //!< transformation matrix for q of lens
complex_t ***q12; //!< quantum noise matrix element mode 1 -> mode 2
......
......@@ -1388,6 +1388,20 @@ int set_k_mirror(int mirror_index) {
}
}
memset(mirror->k11_sqrd_sum, 0, sizeof(double) * inter.num_fields);
memset(mirror->k12_sqrd_sum, 0, sizeof(double) * inter.num_fields);
memset(mirror->k21_sqrd_sum, 0, sizeof(double) * inter.num_fields);
memset(mirror->k22_sqrd_sum, 0, sizeof(double) * inter.num_fields);
for (n = 0; n < inter.num_fields; n++) {
for (m = 0; m < inter.num_fields; m++) {
mirror->k11_sqrd_sum[n] += zabs(mirror->knm.k11[n][m]);
mirror->k12_sqrd_sum[n] += zabs(mirror->knm.k12[n][m]);
mirror->k21_sqrd_sum[n] += zabs(mirror->knm.k21[n][m]);
mirror->k22_sqrd_sum[n] += zabs(mirror->knm.k22[n][m]);
}
}
if(inter.debug & 4096){
double k_sqrd_sum[4] = {0};
......@@ -2024,10 +2038,31 @@ int set_k_beamsplitter(int bs_index) {
bs->knm.k21[n][m] = z_by_x(bs->knm.k21[n][m], turn);
bs->knm.k34[n][m] = z_by_x(bs->knm.k34[n][m], turn);
bs->knm.k43[n][m] = z_by_x(bs->knm.k43[n][m], turn);
}
}
memset(bs->k12_sqrd_sum, 0, sizeof(double) * inter.num_fields);
memset(bs->k21_sqrd_sum, 0, sizeof(double) * inter.num_fields);
memset(bs->k34_sqrd_sum, 0, sizeof(double) * inter.num_fields);
memset(bs->k43_sqrd_sum, 0, sizeof(double) * inter.num_fields);
memset(bs->k13_sqrd_sum, 0, sizeof(double) * inter.num_fields);
memset(bs->k31_sqrd_sum, 0, sizeof(double) * inter.num_fields);
memset(bs->k24_sqrd_sum, 0, sizeof(double) * inter.num_fields);
memset(bs->k42_sqrd_sum, 0, sizeof(double) * inter.num_fields);
for (n = 0; n < inter.num_fields; n++) {
for (m = 0; m < inter.num_fields; m++) {
bs->k12_sqrd_sum[n] += zabs(bs->knm.k12[n][m]);
bs->k21_sqrd_sum[n] += zabs(bs->knm.k21[n][m]);
bs->k34_sqrd_sum[n] += zabs(bs->knm.k34[n][m]);
bs->k43_sqrd_sum[n] += zabs(bs->knm.k43[n][m]);
bs->k13_sqrd_sum[n] += zabs(bs->knm.k13[n][m]);
bs->k31_sqrd_sum[n] += zabs(bs->knm.k31[n][m]);
bs->k24_sqrd_sum[n] += zabs(bs->knm.k24[n][m]);
bs->k42_sqrd_sum[n] += zabs(bs->knm.k42[n][m]);
}
}
if ((bs->knm_flags & PRINT_COEFFS) && bs->map_merged.save_knm_matrices) {
sprintf(buf, "%s_mergedalt", bs->map_merged.filename);
write_bs_knm_to_matrix_file(buf, &(bs->knm));
......@@ -2233,6 +2268,18 @@ int set_k_space(int space_index) {
}
}
}
int n, m;
memset(space->k12_sqrd_sum, 0, sizeof(double) * inter.num_fields);
memset(space->k21_sqrd_sum, 0, sizeof(double) * inter.num_fields);
for (n = 0; n < inter.num_fields; n++) {
for (m = 0; m < inter.num_fields; m++) {
space->k12_sqrd_sum[n] += zabs(space->k12[n][m]);
space->k21_sqrd_sum[n] += zabs(space->k21[n][m]);
}
}
return (space->mismatching);
}
......@@ -2436,6 +2483,18 @@ int set_k_lens(int lens_index) {
}
}
int n, m;
memset(lens->k12_sqrd_sum, 0, sizeof(double) * inter.num_fields);
memset(lens->k21_sqrd_sum, 0, sizeof(double) * inter.num_fields);
for (n = 0; n < inter.num_fields; n++) {
for (m = 0; m < inter.num_fields; m++) {
lens->k12_sqrd_sum[n] += zabs(lens->k12[n][m]);
lens->k21_sqrd_sum[n] += zabs(lens->k21[n][m]);
}
}
return (lens->mismatching);
}
......@@ -2623,7 +2682,18 @@ int set_k_modulator(int modulator_index) {
message("\n");
}
}
int n, m;
memset(modulator->k12_sqrd_sum, 0, sizeof(double) * inter.num_fields);
memset(modulator->k21_sqrd_sum, 0, sizeof(double) * inter.num_fields);
for (n = 0; n < inter.num_fields; n++) {
for (m = 0; m < inter.num_fields; m++) {
modulator->k12_sqrd_sum[n] += zabs(modulator->k12[n][m]);
modulator->k21_sqrd_sum[n] += zabs(modulator->k21[n][m]);
}
}
return (modulator->mismatching);
}
......@@ -2932,6 +3002,22 @@ int set_k_diode(int diode_index) {
}
}
int n, m;
memset(diode->k12_sqrd_sum, 0, sizeof(double) * inter.num_fields);
memset(diode->k21_sqrd_sum, 0, sizeof(double) * inter.num_fields);
memset(diode->k23_sqrd_sum, 0, sizeof(double) * inter.num_fields);
memset(diode->k32_sqrd_sum, 0, sizeof(double) * inter.num_fields);
for (n = 0; n < inter.num_fields; n++) {
for (m = 0; m < inter.num_fields; m++) {
diode->k12_sqrd_sum[n] += zabs(diode->k12[n][m]);
diode->k21_sqrd_sum[n] += zabs(diode->k21[n][m]);
diode->k23_sqrd_sum[n] += zabs(diode->k32[n][m]);
diode->k32_sqrd_sum[n] += zabs(diode->k23[n][m]);
}
}
return (diode->mismatching);
}
......
......@@ -1612,8 +1612,8 @@ void check_map_size(surface_map_t *map, complex_t qx1, complex_t qy1, complex_t
diam_beam = 2 * max(r0x, r0y);
//fprintf(stdout,"cols = %d, center =%g\n",map->cols, map->x0);
r0x = min(fabs(map->cols - (map->x0 + 0.5)), abs(map->x0 + 0.5)) * map->xstep;
r0y = min(fabs(map->cols - (map->y0 + 0.5)), abs(map->y0 + 0.5)) * map->ystep;
r0x = min(fabs(map->cols - (map->x0 + 0.5)), fabs(map->x0 + 0.5)) * map->xstep;
r0y = min(fabs(map->cols - (map->y0 + 0.5)), fabs(map->y0 + 0.5)) * map->ystep;
diam_map = 2 * min(r0x, r0y);
if (diam_map < 4 * diam_beam) {
......
......@@ -1465,6 +1465,11 @@ int allocate_memory_for_mirror_list(long int *bytes) {
mirror_t *m = &inter.mirror_list[i];
int err;
m->k11_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
m->k12_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
m->k21_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
m->k22_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
err = alloc_field_freq_ptrs_mod(&m->a11f, bytes, 4, 4, 5);
if(err) return err;
......@@ -1789,6 +1794,15 @@ int allocate_memory_for_beamsplitter_list(long int *bytes) {
beamsplitter_t *bs = &inter.bs_list[i];
int err;
bs->k12_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
bs->k21_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
bs->k34_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
bs->k43_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
bs->k13_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
bs->k31_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
bs->k24_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
bs->k42_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
err = alloc_field_freq_ptrs_mod(&bs->a11f, bytes, 7, 7, 8);
if(err) return err;
......@@ -2069,7 +2083,10 @@ int allocate_memory_for_space_list(long int *bytes) {
space_t *s = &inter.space_list[i];
int err;
s->k12_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
s->k21_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
err = alloc_field_freq_ptrs(&s->a12f, bytes, 10, 10, 11);
if(err) return err;
......@@ -2297,7 +2314,10 @@ int allocate_memory_for_modulator_list(long int *bytes) {
modulator_t *m = &inter.modulator_list[i];
int err;
m->k12_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
m->k21_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
err = alloc_field_freq_ptrs_mod(&m->a12f, bytes, 17, 17, 18);
if(err) return err;
......@@ -2433,7 +2453,10 @@ int allocate_memory_for_diode_list(long int *bytes) {
diode_t *d = &inter.diode_list[i];
int err;
d->k12_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
d->k21_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
err = alloc_field_freq_ptrs(&d->a12f, bytes, 24, 24, 25);
if(err) return err;
......@@ -2477,6 +2500,9 @@ int allocate_memory_for_lens_list(long int *bytes) {
lens_t *l = &inter.lens_list[i];
int err;
l->k12_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
l->k21_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
err = alloc_field_freq_ptrs(&l->a12f, bytes, 27, 27, 28);
if(err) return err;
......
......@@ -339,7 +339,10 @@ void find_quantum_components(){
// if the mirror is lossy or has an open port, add it.
// TODO - need to check for maps here too as they are also lossy technically
if((1.0 - m->R - m->T > 1e-14) || is_var_tuned(&m->R) || is_var_tuned(&m->T) || (is_open_port(m->node1_index) ^ is_open_port(m->node2_index))){
if((1.0 - m->R - m->T > 1e-14) || is_var_tuned(&m->R)
|| is_var_tuned(&m->T)
|| is_open_port(m->node1_index) || is_open_port(m->node2_index)
|| inter.tem_is_set){
inter.quantum_components[inter.num_quantum_components++] = get_overall_component_index(MIRROR, m->comp_index);
}
}
......@@ -349,7 +352,13 @@ void find_quantum_components(){
// if the mirror is lossy or has an open port, add it.
// TODO - need to check for maps here too as they are also lossy technically
if((bs->R+bs->T < 1) || is_var_tuned(&bs->R) || is_var_tuned(&bs->T) || is_open_port(bs->node1_index) || is_open_port(bs->node2_index) || is_open_port(bs->node3_index) || is_open_port(bs->node4_index)){
if((bs->R+bs->T < 1) || is_var_tuned(&bs->R)
|| is_var_tuned(&bs->T)
|| is_open_port(bs->node1_index)
|| is_open_port(bs->node2_index)
|| is_open_port(bs->node3_index)
|| is_open_port(bs->node4_index)
|| inter.tem_is_set){
inter.quantum_components[inter.num_quantum_components++] = get_overall_component_index(BEAMSPLITTER, bs->comp_index);
}
}
......@@ -417,7 +426,8 @@ void __check_quant_noise_for_node(bool count, int node_index, size_t *num, const
(*num)++;
}
if(loss != 0 || loss_tuned){
// of any loss or if HOM are used as coupling matrix might have losses in
if(loss != 0 || loss_tuned || inter.tem_is_set){
if(!count) {
set_noise_port(&inter.qnoise_in_list[*num], comp_list_index, comp_type, COMPONENT_LOSS, node_index);
......@@ -743,6 +753,80 @@ void create_qnoise_input_matrix(matrix_ccs_t *M){
//dump_matrix_ccs(M, "qnoise_input_%i.dat");
}
void get_knm_loss(qnoise_input_t *qin, double *loss) {
int k;
if(qin->component_type == MIRROR){
mirror_t *mirror = &inter.mirror_list[qin->component_index];
for (k = 0; k < inter.num_fields; k++) {
if(mirror->node1_index == qin->node->list_index){
loss[k] = mirror->R*(1 - mirror->k11_sqrd_sum[k]) + mirror->T*(1 - mirror->k21_sqrd_sum[k]);
} else {
loss[k] = mirror->R*(1 - mirror->k22_sqrd_sum[k]) + mirror->T*(1 - mirror->k12_sqrd_sum[k]);
}
}
} else if(qin->component_type == BEAMSPLITTER){
beamsplitter_t *bs = &inter.bs_list[qin->component_index];
for (k = 0; k < inter.num_fields; k++) {
if(bs->node1_index == qin->node->list_index){
loss[k] = bs->R*(1 - bs->k21_sqrd_sum[k]) + bs->T*(1 - bs->k31_sqrd_sum[k]);
} else if(bs->node2_index == qin->node->list_index){
loss[k] = bs->R*(1 - bs->k12_sqrd_sum[k]) + bs->T*(1 - bs->k42_sqrd_sum[k]);
} else if(bs->node3_index == qin->node->list_index){
loss[k] = bs->R*(1 - bs->k43_sqrd_sum[k]) + bs->T*(1 - bs->k13_sqrd_sum[k]);
} else if(bs->node4_index == qin->node->list_index){
loss[k] = bs->R*(1 - bs->k34_sqrd_sum[k]) + bs->T*(1 - bs->k24_sqrd_sum[k]);
}
}
} else if(qin->component_type == MODULATOR){
modulator_t *mod = &inter.modulator_list[qin->component_index];
for (k = 0; k < inter.num_fields; k++) {
if(mod->node1_index == qin->node->list_index){
loss[k] = 1 - mod->k21_sqrd_sum[k];
} else if(mod->node2_index == qin->node->list_index){
loss[k] = 1 - mod->k12_sqrd_sum[k];
}
}
} else if(qin->component_type == LENS){
lens_t *lens = &inter.lens_list[qin->component_index];
for (k = 0; k < inter.num_fields; k++) {
if(lens->node1_index == qin->node->list_index){
loss[k] = 1 - lens->k21_sqrd_sum[k];
} else if(lens->node2_index == qin->node->list_index){
loss[k] = 1 - lens->k12_sqrd_sum[k];
}
}
} else if(qin->component_type == SPACE){
space_t *s = &inter.space_list[qin->component_index];
for (k = 0; k < inter.num_fields; k++) {
if(s->node1_index == qin->node->list_index){
loss[k] = 1 - s->k21_sqrd_sum[k];
} else if(s->node2_index == qin->node->list_index){
loss[k] = 1 - s->k12_sqrd_sum[k];
}
}
} else if(qin->component_type == DIODE){
diode_t *diode = &inter.diode_list[qin->component_index];
for (k = 0; k < inter.num_fields; k++) {
if(diode->node1_index == qin->node->list_index){
loss[k] = 1 - diode->k21_sqrd_sum[k];
} else if(diode->node2_index == qin->node->list_index){
loss[k] = 2 - diode->k12_sqrd_sum[k] - diode->k32_sqrd_sum[k];
} else if(diode->node3_index == qin->node->list_index){
loss[k] = 1 - diode->k23_sqrd_sum[k];
}
}
} else
bug_error("Not handled\n");
}
int count;
/**
......@@ -881,12 +965,32 @@ void fill_qnoise_input_matrix(){
else
bug_error("not handled");
/* If this component has some loss then we need to fill this in */
for(j=0; j<inter.num_signal_freqs; j++){
for(k=0; k<inter.num_fields; k++){
qin->M_diag_noise_inputs[j*inter.num_fields + k]->re = loss/2.0;
qin->M_diag_noise_inputs[j*inter.num_fields + k]->im = 0;
}
// Here the loss computation differs depending on whether we have
// higher order modes or not.
if (inter.tem_is_set){
// With HOM the coupling matrices *might*
// introduce a per mode losses because power is scattered out
// so vacuum noise must be injected back in again.
double *loss_knm = (double*)calloc(sizeof(double), inter.num_fields);
get_knm_loss(qin, loss_knm);
for(j=0; j<inter.num_signal_freqs; j++){
for(k=0; k<inter.num_fields; k++){
qin->M_diag_noise_inputs[j*inter.num_fields + k]->re = (loss_knm[k] + loss)/2.0;
qin->M_diag_noise_inputs[j*inter.num_fields + k]->im = 0;
}
}
free(loss_knm);
} else {
/* If this component has some loss then we need to fill this in */
for(j=0; j<inter.num_signal_freqs; j++){
for(k=0; k<inter.num_fields; k++){
qin->M_diag_noise_inputs[j*inter.num_fields + k]->re = loss/2.0;
qin->M_diag_noise_inputs[j*inter.num_fields + k]->im = 0;
}
}
}
} else {
bug_error("unhandled value %i", qin->type);
......@@ -1508,8 +1612,7 @@ void compute_squeezing_factor(light_out_t *sd, double *r, double *phi, double *d
*phi = - 0.5 * zphase( b[0][1]) + zphase(b[0][0])/2;
double rr = log(min(eig1, eig2))/2.0;
*r = rr * R2DB;
*r = rr * R2DB;
}
}
......
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