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

Adding additional option for three port isolator to use for injecting squeezed...

Adding additional option for three port isolator to use for injecting squeezed light whilst not diverting the return signal field. Fixing isolator beam tracing hopefully, this needs checking though. I added an additional pass condition to trace_beam as it was complaining that the beam wasn't traced to all nodes but all nodes had a parameter set. Surface motion gouy phase still needs work
parent cd86bb6c
......@@ -1826,6 +1826,8 @@ typedef struct diode {
bitflag knm_flags;
int nnz_count;
int option;
} diode_t;
//! information describing a lens
......
......@@ -1396,7 +1396,12 @@ int set_k_mirror(int mirror_index) {
}
// 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]));
// warn("s1: %i%i->%i%i = %s\n",n1,m1,n2,m2, complex_form(mirror->knm_surf_x_a_1[k][n][m]));
// warn("s2: %i%i->%i%i = %s\n",n1,m1,n2,m2, complex_form(mirror->knm_surf_x_a_2[k][n][m]));
// warn("1i: %i%i->%i%i = %s\n",n1,m1,n2,m2, complex_form(mirror->knm_surf_motion_1i[k][n][m]));
// warn("1o: %i%i->%i%i = %s\n",n1,m1,n2,m2, complex_form(mirror->knm_surf_motion_1o[k][n][m]));
// warn("2i: %i%i->%i%i = %s\n",n1,m1,n2,m2, complex_form(mirror->knm_surf_motion_2i[k][n][m]));
// warn("2o: %i%i->%i%i = %s\n",n1,m1,n2,m2, complex_form(mirror->knm_surf_motion_2o[k][n][m]));
// }
}
}
......@@ -2780,17 +2785,17 @@ int set_k_diode(int diode_index) {
//node2 = inter.node_list[node2_index];
nr3 = *node3.n;
if (node3.component_index == component) {
qx3 = node3.qx;
qy3 = node3.qy;
} else {
// reverse q so that beam goes _from_ the diode to node2
// reverse q so that beam goes _from_ the diode to node3
qx3 = cminus(cconj(node3.qx));
qy3 = cminus(cconj(node3.qy));
} else {
qx3 = node3.qx;
qy3 = node3.qy;
}
}
if ((nr1 != nr2 && !node1.gnd_node && !node2.gnd_node) || (nr3 != nr2 && !node3.gnd_node && !node2.gnd_node)) {
gerror("index of refraction not consistent at diode '%s'\n", diode->name);
gerror("index of refraction not consistent at diode '%s' (%g, %g, %g)\n", diode->name, nr1, nr2, nr3);
}
if (inter.debug & 32 || inter.trace & 64) {
......@@ -4356,10 +4361,9 @@ int trace_cavity(ABCD_t *Mlocal_x, ABCD_t *Mlocal_y,
message("trace_cavity(): new node %s (index %d)\n",
get_node_name(node_indices[i]), node_indices[i]);
}
err_x = component_matrix(&transform_x, component, node1_index,
node_indices[i], TANGENTIAL);
err_y = component_matrix(&transform_y, component, node1_index,
node_indices[i], SAGITTAL);
err_x = component_matrix(&transform_x, component, node1_index, node_indices[i], TANGENTIAL);
err_y = component_matrix(&transform_y, component, node1_index, node_indices[i], SAGITTAL);
if (NOT err_x && NOT err_y) {
temp_x = multiply_ABCD(*Mlocal_x, inv_ABCD(transform_x));
......@@ -4383,7 +4387,7 @@ int trace_cavity(ABCD_t *Mlocal_x, ABCD_t *Mlocal_y,
*Mlocal_y = temp_y;
}
} else {
bug_error("errxy3");
bug_error("errxy3 %i %i", err_x, err_y);
}
}
}
......@@ -4419,6 +4423,9 @@ void trace_beam(int startnode) {
which_components(startnode, &component1_index, &component2_index);
//warn("%s\n", get_component_name(component1_index));
//warn("%s\n", get_component_name(component2_index));
if (component1_index != NOT_FOUND) {
set_qbase(component1_index);
}
......@@ -4446,13 +4453,16 @@ void trace_beam(int startnode) {
// nodes, if not, then barf.
if ((n_index + inter.num_dump_nodes) != inter.num_nodes) {
fprintf(stderr,
"Interferometer may be broken;"
" no parameters were set for the\nfollowing nodes:\n");
"Interferometer may be broken (%i/%i);"
" no parameters were set for the\nfollowing nodes:\n", n_index + inter.num_dump_nodes, inter.num_nodes);
int i;
bool missing = false;
for (i = 0; i < inter.num_nodes; i++) {
if ((NOT inter.node_list[i].q_is_set) && (NOT inter.node_list[i].gnd_node)) {
fprintf(stderr, "%s ", inter.node_list[i].name);
missing = true;
j++;
if (j > 6) {
j = 0;
......@@ -4460,8 +4470,11 @@ void trace_beam(int startnode) {
}
}
}
fprintf(stderr, "\n");
gerror("Could not completely trace beam.\n");
if (missing)
gerror("Could not completely trace beam.\n");
}
}
......@@ -4493,6 +4506,10 @@ void tracing(int startnode, int start_comp) {
assert(start_comp < inter.num_components);
which_components(startnode, &component1_index, &component2_index);
//warn("%s\n", get_component_name(component1_index));
//warn("%s\n", get_component_name(component2_index));
if (component2_index == start_comp) {
component_index = component1_index;
} else {
......@@ -4770,10 +4787,11 @@ int component_matrix(ABCD_t *matrix, int component_index,
return (1);
}
// there is no connection between node 2 and node 3
if((node1_index == diode->node2_index && node2_index == diode->node3_index)
|| (node1_index == diode->node3_index && node2_index == diode->node2_index))
// there is no connection between node 1 and node 3
if( (node1_index == diode->node1_index && node2_index == diode->node3_index)
|| (node1_index == diode->node3_index && node2_index == diode->node1_index))
return (1);
*matrix = (diode->qq);
break;
......@@ -4932,8 +4950,7 @@ void set_qbase(int component_index) {
for (j = 0; j < no; j++) {
node_j = inter.node_list[node_indices[j]];
// set q if q is not yet set
if (NOT node_j.q_is_set &&
NOT node_j.gnd_node) { // || get_component_type(component)==SPACE)
if (NOT node_j.q_is_set && NOT node_j.gnd_node) { // || get_component_type(component)==SPACE)
err = component_matrix(&transform, component_index,
node_indices[i], node_indices[j], TANGENTIAL);
if (!err) {
......@@ -5065,184 +5082,190 @@ void which_nodes(int component_index, int current_node,
bug_error("component index less than zero");
}
mirror_t mirror = inter.mirror_list[component_index];
node_t node1 = inter.node_list[mirror.node1_index];
node_t node2 = inter.node_list[mirror.node2_index];
mirror_t *mirror = &inter.mirror_list[component_index];
node_t *node1 = &inter.node_list[mirror->node1_index];
node_t *node2 = &inter.node_list[mirror->node2_index];
if (NOT node1.gnd_node) {
node_indices[no++] = mirror.node1_index;
if (NOT node1->gnd_node) {
node_indices[no++] = mirror->node1_index;
}
if (NOT node2.gnd_node) {
node_indices[no++] = mirror.node2_index;
if (NOT node2->gnd_node) {
node_indices[no++] = mirror->node2_index;
}
}// experimental: BS returns only the two available nodes
// depending on the node of the incident field
else if ((component_index -= inter.num_mirrors) < inter.num_beamsplitters) {
beamsplitter_t beamsplitter = inter.bs_list[component_index];
beamsplitter_t *beamsplitter = &inter.bs_list[component_index];
node_t node1 = inter.node_list[beamsplitter.node1_index];
node_t node2 = inter.node_list[beamsplitter.node2_index];
node_t node3 = inter.node_list[beamsplitter.node3_index];
node_t node4 = inter.node_list[beamsplitter.node4_index];
node_t *node1 = &inter.node_list[beamsplitter->node1_index];
node_t *node2 = &inter.node_list[beamsplitter->node2_index];
node_t *node3 = &inter.node_list[beamsplitter->node3_index];
node_t *node4 = &inter.node_list[beamsplitter->node4_index];
if (current_node == GND_NODE ||
inter.node_list[current_node].gnd_node) {
if (NOT node1.gnd_node) {
node_indices[no++] = beamsplitter.node1_index;
if (NOT node1->gnd_node) {
node_indices[no++] = beamsplitter->node1_index;
}
if (NOT node2.gnd_node) {
node_indices[no++] = beamsplitter.node2_index;
if (NOT node2->gnd_node) {
node_indices[no++] = beamsplitter->node2_index;
}
if (NOT node3.gnd_node) {
node_indices[no++] = beamsplitter.node3_index;
if (NOT node3->gnd_node) {
node_indices[no++] = beamsplitter->node3_index;
}
if (NOT node4.gnd_node) {
node_indices[no++] = beamsplitter.node4_index;
if (NOT node4->gnd_node) {
node_indices[no++] = beamsplitter->node4_index;
}
} else {
if (current_node == beamsplitter.node1_index) {
if (NOT node2.gnd_node) {
node_indices[no++] = beamsplitter.node2_index;
if (current_node == beamsplitter->node1_index) {
if (NOT node2->gnd_node) {
node_indices[no++] = beamsplitter->node2_index;
}
if (NOT node3.gnd_node) {
node_indices[no++] = beamsplitter.node3_index;
if (NOT node3->gnd_node) {
node_indices[no++] = beamsplitter->node3_index;
}
}
if (current_node == beamsplitter.node2_index) {
if (NOT node1.gnd_node) {
node_indices[no++] = beamsplitter.node1_index;
if (current_node == beamsplitter->node2_index) {
if (NOT node1->gnd_node) {
node_indices[no++] = beamsplitter->node1_index;
}
if (NOT node4.gnd_node) {
node_indices[no++] = beamsplitter.node4_index;
if (NOT node4->gnd_node) {
node_indices[no++] = beamsplitter->node4_index;
}
}
if (current_node == beamsplitter.node3_index) {
if (NOT node4.gnd_node) {
node_indices[no++] = beamsplitter.node4_index;
if (current_node == beamsplitter->node3_index) {
if (NOT node4->gnd_node) {
node_indices[no++] = beamsplitter->node4_index;
}
if (NOT node1.gnd_node) {
node_indices[no++] = beamsplitter.node1_index;
if (NOT node1->gnd_node) {
node_indices[no++] = beamsplitter->node1_index;
}
}
if (current_node == beamsplitter.node4_index) {
if (NOT node3.gnd_node) {
node_indices[no++] = beamsplitter.node3_index;
if (current_node == beamsplitter->node4_index) {
if (NOT node3->gnd_node) {
node_indices[no++] = beamsplitter->node3_index;
}
if (NOT node2.gnd_node) {
node_indices[no++] = beamsplitter.node2_index;
if (NOT node2->gnd_node) {
node_indices[no++] = beamsplitter->node2_index;
}
}
}
} else if ((component_index -= inter.num_beamsplitters) < inter.num_gratings) {
grating_t grating = inter.grating_list[component_index];
grating_t *grating = &inter.grating_list[component_index];
node_t node1 = inter.node_list[grating.node1_index];
node_t node2 = inter.node_list[grating.node2_index];
node_t node3 = inter.node_list[grating.node3_index];
node_t node4 = inter.node_list[grating.node4_index];
node_t *node1 = &inter.node_list[grating->node1_index];
node_t *node2 = &inter.node_list[grating->node2_index];
node_t *node3 = &inter.node_list[grating->node3_index];
node_t *node4 = &inter.node_list[grating->node4_index];
if (NOT node1.gnd_node) {
node_indices[no++] = grating.node1_index;
if (NOT node1->gnd_node) {
node_indices[no++] = grating->node1_index;
}
if (NOT node2.gnd_node) {
node_indices[no++] = grating.node2_index;
if (NOT node2->gnd_node) {
node_indices[no++] = grating->node2_index;
}
switch (grating.num_of_ports) {
switch (grating->num_of_ports) {
case 3:
if (NOT node3.gnd_node) {
node_indices[no++] = grating.node3_index;
if (NOT node3->gnd_node) {
node_indices[no++] = grating->node3_index;
}
case 4:
if (NOT node4.gnd_node) {
node_indices[no++] = grating.node4_index;
if (NOT node4->gnd_node) {
node_indices[no++] = grating->node4_index;
}
}
} else if ((component_index -= inter.num_gratings) < inter.num_spaces) {
space_t space = inter.space_list[component_index];
space_t *space = &inter.space_list[component_index];
node_t node1 = inter.node_list[space.node1_index];
node_t node2 = inter.node_list[space.node2_index];
node_t *node1 = &inter.node_list[space->node1_index];
node_t *node2 = &inter.node_list[space->node2_index];
if (NOT node1.gnd_node) {
node_indices[no++] = space.node1_index;
if (NOT node1->gnd_node) {
node_indices[no++] = space->node1_index;
}
if (NOT node2.gnd_node) {
node_indices[no++] = space.node2_index;
if (NOT node2->gnd_node) {
node_indices[no++] = space->node2_index;
}
} else if ((component_index -= inter.num_spaces) < inter.num_lenses) {
lens_t lens = inter.lens_list[component_index];
lens_t *lens = &inter.lens_list[component_index];
node_t node1 = inter.node_list[lens.node1_index];
node_t node2 = inter.node_list[lens.node2_index];
node_t *node1 = &inter.node_list[lens->node1_index];
node_t *node2 = &inter.node_list[lens->node2_index];
if (NOT node1.gnd_node) {
node_indices[no++] = lens.node1_index;
if (NOT node1->gnd_node) {
node_indices[no++] = lens->node1_index;
}
if (NOT node2.gnd_node) {
node_indices[no++] = lens.node2_index;
if (NOT node2->gnd_node) {
node_indices[no++] = lens->node2_index;
}
} else if ((component_index -= inter.num_lenses) < inter.num_diodes) {
diode_t diode = inter.diode_list[component_index];
diode_t *diode = &inter.diode_list[component_index];
node_t node1 = inter.node_list[diode.node1_index];
node_t node2 = inter.node_list[diode.node2_index];
node_t node3 = inter.node_list[diode.node3_index];
node_t *node1 = &inter.node_list[diode->node1_index];
node_t *node2 = &inter.node_list[diode->node2_index];
node_t *node3 = &inter.node_list[diode->node3_index];
if (NOT node1.gnd_node) {
node_indices[no++] = diode.node1_index;
}
if (NOT node2.gnd_node) {
node_indices[no++] = diode.node2_index;
}
if (NOT node3.gnd_node) {
node_indices[no++] = diode.node3_index;
if (current_node == GND_NODE || inter.node_list[current_node].gnd_node) {
if (NOT node1->gnd_node) node_indices[no++] = diode->node1_index;
if (NOT node2->gnd_node) node_indices[no++] = diode->node2_index;
if (NOT node3->gnd_node) node_indices[no++] = diode->node3_index;
} else {
// otherwise we need to return only the nodes that the current node
// can reach
if (current_node == diode->node1_index){
if (NOT node2->gnd_node) node_indices[no++] = diode->node2_index;
} else if (current_node == diode->node2_index){
if (NOT node1->gnd_node) node_indices[no++] = diode->node1_index;
if (NOT node3->gnd_node) node_indices[no++] = diode->node3_index;
} else if (current_node == diode->node3_index){
if (NOT node2->gnd_node) node_indices[no++] = diode->node2_index;
} else
bug_error("unhandled (%i)\n", current_node);
}
} else if ((component_index -= inter.num_diodes) < inter.num_light_inputs) {
light_in_t light_in = inter.light_in_list[component_index];
light_in_t *light_in = &inter.light_in_list[component_index];
node_t light_in_node = inter.node_list[light_in.node_index];
node_t light_in_node = inter.node_list[light_in->node_index];
if (NOT light_in_node.gnd_node) {
node_indices[no++] = light_in.node_index;
node_indices[no++] = light_in->node_index;
}
} else if ((component_index -= inter.num_light_inputs) < inter.num_modulators) {
modulator_t modulator = inter.modulator_list[component_index];
modulator_t *modulator = &inter.modulator_list[component_index];
node_t node1 = inter.node_list[modulator.node1_index];
node_t node2 = inter.node_list[modulator.node2_index];
node_t *node1 = &inter.node_list[modulator->node1_index];
node_t *node2 = &inter.node_list[modulator->node2_index];
if (NOT node1.gnd_node) {
node_indices[no++] = modulator.node1_index;
if (NOT node1->gnd_node) {
node_indices[no++] = modulator->node1_index;
}
if (NOT node2.gnd_node) {
node_indices[no++] = modulator.node2_index;
if (NOT node2->gnd_node) {
node_indices[no++] = modulator->node2_index;
}
} else if ((component_index -= inter.num_modulators) < inter.num_signals) {
bug_error("which_node 3");
no = 2;
} else if ((component_index -= inter.num_signals) < inter.num_outputs) {
output_data_t output_data = inter.output_data_list[component_index];
output_data_t *output_data = &inter.output_data_list[component_index];
node_t output_data_node = inter.node_list[output_data.node_index];
node_t output_data_node = inter.node_list[output_data->node_index];
if (NOT output_data_node.gnd_node) {
node_indices[no++] = output_data.node_index;
node_indices[no++] = output_data->node_index;
}
} else {
bug_error("component index outside correct range");
......
......@@ -943,6 +943,10 @@ void calc_mirror_knm_surf_motions_map(mirror_t *m, double nr1, double nr2, compl
#else
integrate_riemann_sum_surf(&value, &p);
#endif
// surface integration is just the real kernal of the hermite
// polynomials, have to add in the gouy phase here.
// Even though we reverse gouy it later it needs to be in
// here for merging with any static distortions
complex_t knm = z_by_xphr(complex_1, value, (p.n1 - p.n2) * gx11 + (p.m1 - p.m2) * gy11);
//warn("%i%i->%i%i %s\n", p.n1, p.m1, p.n2, p.m2, complex_form15(knm));
......@@ -954,7 +958,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] = knm; //cconj(z_by_x(knm, msign(p.n1+p.n2)));
m->knm_surf_motion_1i[k][in][out] = 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));
......@@ -1001,16 +1005,7 @@ void calc_mirror_knm_surf_motions_map(mirror_t *m, double nr1, double nr2, compl
m->knm_surf_motion_2o[k][out][in] = cconj(m->knm_surf_motion_2o[k][in][out]);
m->knm_surf_motion_2i[k][out][in] = cconj(m->knm_surf_motion_2i[k][in][out]);
}
} else {
m->knm_surf_motion_2o[k][in][out] = z_by_x(complex_1, value* msign(p.n1+p.n2));
m->knm_surf_motion_2i[k][in][out] = cconj(z_by_x(complex_1, value));
if(out!=in) {
m->knm_surf_motion_2o[k][out][in] = cconj(m->knm_surf_motion_2o[k][in][out]);
m->knm_surf_motion_2i[k][out][in] = cconj(m->knm_surf_motion_2i[k][in][out]);
}
}
}
if(out!=in)
current_coeff+=2; // add two because we also compute the transposed elements too
......
......@@ -446,6 +446,7 @@ void pre_scan(FILE *fp) {
mem.num_node_names++;
mem.num_nodes++;
} else if (string_matches(s, "isol ") ||
string_matches(s, "isol* ") ||
string_matches(s, "diode ")) {
mem.num_diodes++;
mem.num_node_names++;
......@@ -2481,6 +2482,8 @@ int allocate_memory_for_diode_list(long int *bytes) {
d->k12_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
d->k21_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
d->k23_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
d->k32_sqrd_sum = (double*) malloc(mem.num_fields * sizeof(double));
err = alloc_field_freq_ptrs(&d->a12f, bytes, 24, 24, 25);
if(err) return err;
......
......@@ -492,9 +492,9 @@ size_t count_quantum_components(bool count){
} else if(type == DIODE){
diode_t *m = &inter.diode_list[j];
// TODO Handle loss at diode properly
__check_quant_noise_for_node(count, m->node1_index, &num, m->name, j, type, 0, false);
__check_quant_noise_for_node(count, m->node2_index, &num, m->name, j, type, 0, false);
__check_quant_noise_for_node(count, m->node3_index, &num, m->name, j, type, 0, false);
//__check_quant_noise_for_node(count, m->node1_index, &num, m->name, j, type, 0, false);
//__check_quant_noise_for_node(count, m->node2_index, &num, m->name, j, type, 0, false);
//__check_quant_noise_for_node(count, m->node3_index, &num, m->name, j, type, 0, false);
} else if(type == LENS){
lens_t *m = &inter.lens_list[j];
......@@ -962,8 +962,11 @@ void fill_qnoise_input_matrix(){
loss = 1 - inter.mirror_list[qin->component_index].R - inter.mirror_list[qin->component_index].T;
else if(qin->component_type == BEAMSPLITTER)
loss = 1 - inter.bs_list[qin->component_index].R - inter.bs_list[qin->component_index].T;
else
bug_error("not handled");
else if(qin->component_type == DIODE)
loss = 1 - inter.diode_list[qin->component_index].loss;
else {
bug_error("not handled (%i)", qin->component_type);
}
// Here the loss computation differs depending on whether we have
// higher order modes or not.
......
......@@ -323,6 +323,7 @@ void check_all_commands(FILE *fp){
(strncasecmp(s, "laser ", 6) == 0) ||
(strncasecmp(s, "light ", 6) == 0) ||
(strncasecmp(s, "isol ", 5) == 0) ||
(strncasecmp(s, "isol* ", 6) == 0) ||
(strncasecmp(s, "diode ", 6) == 0) ||
(strncasecmp(s, "variable ", 9) == 0) ||
(strncasecmp(s, "var ", 4) == 0) ||
......@@ -554,7 +555,7 @@ void read_file(FILE *fp) {
read_verbose_grating(s);
} else if (strncasecmp(s, "gr", 2) == 0) {
read_grating(s);
} else if (strncasecmp(s, "isol ", 5) == 0) {
} else if (strncasecmp(s, "isol ", 5) == 0 || strncasecmp(s, "isol* ", 6) == 0) {
read_diode(s);
} else if (strncasecmp(s, "diode ", 6) == 0) {
read_diode(s);
......@@ -12353,6 +12354,12 @@ void read_diode(const char *command_string) {
diode->node1_reduced_index = get_reduced_index_for(diode->node1_index);
diode->node2_reduced_index = get_reduced_index_for(diode->node2_index);
if(command_name[strlen(command_name)-1] == '*'){
diode->option = 1;
} else {
diode->option = 0;
}
++inter.num_diodes;
++inter.num_components;
check_diode(diode);
......
......@@ -1025,7 +1025,8 @@ void fill_surf_motion_coupling(int motion_idx, int surface_idx, mirror_t *m, com
// we don't have to worry about the memory alignment, as we are using CCS
// matrix.
cblas_zgemv(CblasRowMajor, CblasNoTrans, inter.num_fields, inter.num_fields,
&fac, &(m->knm_surf_x_a_1[surface_idx][0][0]), inter.num_fields, ac_1i, 1, &complex_0, x_a1u[0], 1);
&fac, &(m->knm_surf_x_a_1[surface_idx][0][0]), inter.num_fields,
ac_1i, 1, &complex_0, x_a1u[0], 1);
// lower sideband creation is simply the conjugate of the upper
// plus we must also propagate the new sidebands away from the mirror
......@@ -2924,6 +2925,7 @@ void fill_matrix_lens_elements(bool conjugate) {
void fill_matrix_diode_elements(bool conjugate) {
// fill matrix for diodes
int diode_index;
for (diode_index = 0; diode_index < inter.num_diodes; diode_index++) {
diode_t *diode;
diode = &(inter.diode_list[diode_index]);
......@@ -2936,38 +2938,72 @@ void fill_matrix_diode_elements(bool conjugate) {
int j, k;
for (j = 0; j < inter.num_fields; j++) {
for (k = 0; k < inter.num_fields; k++) {
switch(diode->option) {
case 0:
if (diode->a12[j][k] != NULL) {
(diode->a12[j][k])->re = -1.0;
(diode->a12[j][k])->im = 0.0;
*(diode->a12[j][k]) = z_by_x(z_by_z(*(diode->a12[j][k]), diode->k12[j][k]), loss);
if(conjugate) diode->a12[j][k]->im *= -1;
}
if (diode->a12[j][k] != NULL) {
(diode->a12[j][k])->re = -1.0;
(diode->a12[j][k])->im = 0.0;
if (diode->a21[j][k] != NULL) {
(diode->a21[j][k])->re = -leak;
(diode->a21[j][k])->im = 0.0;
*(diode->a12[j][k]) = z_by_x(z_by_z(*(diode->a12[j][k]), diode->k12[j][k]), loss);
if(conjugate) diode->a12[j][k]->im *= -1;
}
*(diode->a21[j][k]) = z_by_x(z_by_z(*(diode->a21[j][k]), diode->k21[j][k]), loss);
if(conjugate) diode->a21[j][k]->im *= -1;
}
if (diode->a21[j][k] != NULL) {
(diode->a21[j][k])->re = -leak;
(diode->a21[j][k])->im = 0.0;
if (diode->a23[j][k] != NULL) {
// output from node 2 to node 3 is whatever is not suppressed
(diode->a23[j][k])->re = -out;
(diode->a23[j][k])->im = 0