Commit 3a63690a authored by Daniel Brown's avatar Daniel Brown
Browse files

Wide range of updates to get ROMs working for surface motion computations

parent 0e058a71
......@@ -904,6 +904,7 @@ typedef struct mirror {
UT_array *motion_links;
/** Indices to surface motions applied to this mirror, see inter.surface_motion_map_list*/
bool surface_motions_isROM[MAX_MAPS];
int surface_motions[MAX_MAPS];
transfer_func_t* surface_motions_tf[MAX_MAPS];
......
......@@ -1266,7 +1266,8 @@ int set_k_mirror(int mirror_index) {
if(mirror->num_surface_motions) {
// as we are computing this coupling matrix first we must use
// the outgoing beam parameters.
mirror_calc_knm_surf_motions(mirror, nr1, nr2, knm_q.qxt2_11, knm_q.qyt2_11, knm_q.qxt2_22, knm_q.qyt2_22);
calc_mirror_knm_surf_motions_map(mirror, nr1, nr2, knm_q.qxt2_11, knm_q.qyt2_11, knm_q.qxt2_22, knm_q.qyt2_22);
calc_mirror_knm_surf_motions_rom(mirror, nr1, nr2, knm_q.qxt2_11, knm_q.qyt2_11, knm_q.qxt2_22, knm_q.qyt2_22, astigmatism);
for(k=0; k<mirror->num_surface_motions; k++){
complex_t **A = mirror->knm_surf_motion_1o[k];
......@@ -1274,18 +1275,19 @@ int set_k_mirror(int mirror_index) {
complex_t **C = mirror->knm_surf_x_a_1[k];
// zero the output matrix
memset(&(C[0][0]),0,inter.num_fields*inter.num_fields*sizeof(complex_t));
memset(&(C[0][0]), 0, inter.num_fields * inter.num_fields * sizeof(complex_t));
int n,m,l;
double gx = gouy(knm_q.qxt2_11);
double gy = gouy(knm_q.qyt2_11);
int n, m, l, _n1, _m1;
for (n = 0; n < inter.num_fields; n++) {
get_tem_modes_from_field_index(&_n1, &_m1, n);
for (m = 0; m < inter.num_fields; m++) {
get_tem_modes_from_field_index(&n1, &m1, m);
for (l = 0; l < inter.num_fields; l++) {
get_tem_modes_from_field_index(&n2, &m2, l);
z_inc_z(&C[n][m], z_by_z(B[n][l], z_by_phr(A[l][m], (n1-n2)*gx + (m1-m2)*gy)));
z_inc_z(&C[n][m], z_by_z(B[n][l], A[l][m]));
}
}
}
......@@ -1296,26 +1298,18 @@ int set_k_mirror(int mirror_index) {
memset(&(C[0][0]),0,inter.num_fields*inter.num_fields*sizeof(complex_t));
gx = gouy(knm_q.qxt2_22);
gy = gouy(knm_q.qyt2_22);
for (n = 0; n < inter.num_fields; n++) {
for (m = 0; m < inter.num_fields; m++) {
get_tem_modes_from_field_index(&n1, &m1, m);
for (l = 0; l < inter.num_fields; l++) {
get_tem_modes_from_field_index(&n2, &m2, l);
z_inc_z(&C[n][m], z_by_z(B[n][l], z_by_phr(A[l][m], (n1-n2)*gx + (m1-m2)*gy)));
z_inc_z(&C[n][m], z_by_z(B[n][l], A[l][m]));
}
}
}
}
}
// if(strcmp(mirror->name, "ETMY") == 0 && !mirror->map_rom){
// mirror->knm_no_rgouy.k11[0][0].re = 0.9814582625035128 + 1e-8;
// mirror->knm_no_rgouy.k11[0][0].im = 0;
// }
// Here we adjust the odd reflected horizontal modes. We do this because on
// reflection, just like in a mirror, the horizontal direction is flipped
// left to right. We also reverse the extra gouy phase present and rescale the
......@@ -1336,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][0][0] = rev_gouy(mirror->knm_surf_x_a_1[k][0][0], 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][0][0] = rev_gouy(mirror->knm_surf_x_a_2[k][0][0], 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);
}
}
}
......@@ -1345,7 +1339,7 @@ int set_k_mirror(int mirror_index) {
// Here we store the phases of the TEM00 so that we can 'zero' the phase, as
// typically we do not care for the exact phase here just the relative ones
// between the modes
double p11=0, p12=0, p21=0, p22=0, psm1[MAX_MAPS], psm2[MAX_MAPS];
double p11=0, p12=0, p21=0, p22=0, psm1[MAX_MAPS] = {0}, psm2[MAX_MAPS] = {0};
if(inter.set_tem_phase_zero & 1) {
p11 = zphase(mirror->knm.k11[0][0]);
......@@ -1354,8 +1348,22 @@ int set_k_mirror(int mirror_index) {
p22 = zphase(mirror->knm.k22[0][0]);
for(k=0; k<mirror->num_surface_motions; k++){
psm1[k] = zphase(mirror->knm_surf_x_a_1[k][0][0]);
psm2[k] = zphase(mirror->knm_surf_x_a_2[k][0][0]);
// At times the surface motion will have zero 00->00 coupling
// for example with tilt motions. Thus if it is zero then we don't
// compute a phase, otherwise we just get some phase of noise which
// gives weird results.
if(zmod(mirror->knm_surf_x_a_1[k][0][0]) > ZERO){
psm1[k] = zphase(mirror->knm_surf_x_a_1[k][0][0]);
} else {
psm1[k] = 0.0;
}
if(zmod(mirror->knm_surf_x_a_2[k][0][0]) > ZERO){
psm2[k] = zphase(mirror->knm_surf_x_a_2[k][0][0]);
} else {
psm2[k] = 0.0;
}
}
}
......@@ -4193,7 +4201,7 @@ void compute_cavity_params(void) {
cavity->pole = 0.5 * init.clight / cavity_length / finesse(1.0 - cavity_power);
if (inter.trace & 2) {
message(" finesse : %g, round-trip power loss: %g\n",
message(" finesse : %g, round-trip power loss: %.8g [\%/100]\n",
cavity->finesse, cavity->loss);
message(" opt. length: %sm, FSR: %sHz\n",
xdouble_form(cavity->length), xdouble_form(cavity->FSR));
......
......@@ -1117,9 +1117,7 @@ void fill_mirror_signal_rhs(signal_t *signal, ifo_matrix_vars_t* M_car, ifo_matr
break;
case SURFACE:
warn("%i %i\n", i, mirror->num_motions);
if(i == mirror->num_motions + signal->surface_mode_index){
if(i == (mirror->num_motions-1) + signal->surface_mode_index) {
if(signal->type == SIG_SURF) {
assert(mirror->surface_motions_tf != NULL);
assert(mirror->surface_motions_tf[signal->surface_mode_index] != NULL);
......
......@@ -68,7 +68,7 @@
#define MAX_TOKEN_LEN 256
#define REST_STRING_LEN 82 //!< Length of "rest" of string when reading cmds
#define SMALL 1e-9 //!< The definition of a small number
#define ZERO 1e-17 //!< The definition of zero
#define ZERO 1e-15 //!< The definition of zero
#define BIG 1e100 //!< The definition of a big number
#define ALLOC_LEN 1024 //!< Allocation length
#define FORM_LEN 64 //!< Length of home-made num2str output
......
......@@ -1069,7 +1069,7 @@ void check_refractive_index_at_beamsplitters(void) {
* \todo Need to split this routine up into subroutines
*/
void setup_hermite_gauss_extension(void) {
int tn, tm, i;
int tn, tm, i, j;
// ###########################################################
// Hermite Gauss extension :
......@@ -1173,6 +1173,12 @@ void setup_hermite_gauss_extension(void) {
if(inter.mirror_list[i].map_rom) {
allocate_romhom_workspace(inter.mirror_list[i].map_rom, MIRROR);
}
for(j=0; j < inter.mirror_list[i].num_surface_motions; j++) {
if(inter.mirror_list[i].surface_motions_isROM[j]) {
allocate_romhom_workspace(&inter.rom_maps[inter.mirror_list[i].surface_motions[j]], MIRROR);
}
}
if(inter.mirror_list[i].map_merged.integration_method == NEWTON_COTES) {
allocate_newton_cotes_workspace(&inter.mirror_list[i].map_merged.knm_ws,
......
......@@ -620,7 +620,7 @@ double eval_surface_map(surface_map_t *map, double x, double y){
f /= ((x2-x1)*(y2-y1));
return f * map->scaling;
return f;
}
#if INCLUDE_CUBA == 1
......@@ -688,7 +688,7 @@ void integrate_riemann_sum_surf(double *result, knm_surf_cuba_params_t *p) {
double Hm1 = hermite(p->m1, y1);
double Hm2 = (p->m1==p->m2) ? Hm1 : hermite(p->m2, y1);
double func = p->map->data[j][i] * p->map->scaling;
double func = p->map->data[j][i];
*result += Hn1*Hn2 * Hm1*Hm2* exp(- x1*x1 - y1*y1) * func;
}
}
......@@ -698,7 +698,153 @@ void integrate_riemann_sum_surf(double *result, knm_surf_cuba_params_t *p) {
*result *= da * p->constant;
}
void mirror_calc_knm_surf_motions(mirror_t *m, double nr1, double nr2, complex_t qx1, complex_t qy1, complex_t qx2, complex_t qy2) {
void calc_mirror_knm_surf_motions_rom(mirror_t *mirror, double nr1, double nr2,
complex_t qx11, complex_t qy11,
complex_t qx22, complex_t qy22,
bitflag astigmatism) {
assert(mirror != NULL);
int timer = startTimer("ROMHOM");
if(mirror->x_off != 0.0 || mirror->y_off != 0.0)
warn("ROMHOM can't take into account transverse offsets in mirror position, use normal maps for this or include offset in map when making ROM.", mirror->name);
// boolean value that states whether we should calculate knm using knm without
// having to do the whole integral and just by applying some phase factor
int calc_knm_transpose = (init.calc_knm_transpose && (astigmatism == 0));
int max_m, max_n;
get_tem_modes_from_field_index(&max_n, &max_m, inter.num_fields - 1);
max_m = max_n; // max m is no the m is also the max_n. e.g. we get n=1 m=2 for maxtem 2
// K11
if (CALC_MR_KNM(mirror,11)) {
u_nm_accel_get(acc_11_nr1_1, max_n, max_m, qx11, qy11, nr1);
u_nm_accel_get(acc_11_nr1_2, max_n, max_m, qx11, qy11, nr1);
}
// K22
if (CALC_MR_KNM(mirror,22)) {
u_nm_accel_get(acc_22_nr2_1, max_n, max_m, qx22, qy22, nr2);
u_nm_accel_get(acc_22_nr2_2, max_n, max_m, qx22, qy22, nr2);
}
int num_coeffs = inter.num_fields * inter.num_fields;
int current_coeff = 1;
time_t starttime = time(NULL);
int n, m, n1, m1, n2, m2, l;
set_progress_action_text("Calculating ROMHOM surface knm for %s", mirror->name);
for(l=0; l < mirror->num_surface_motions; l++){
if(!mirror->surface_motions_isROM[l])
continue;
rom_map_t *rom = &inter.rom_maps[mirror->surface_motions[l]];
knm_workspace_t *ws11 = &rom->roq11.knm_ws;
knm_workspace_t *ws22 = &rom->roq22.knm_ws;
if (CALC_MR_KNM(mirror, 11) && rom->roq11.enabled){
double wx11 = w0_size(qx11, nr1);
double wy11 = w0_size(qy11, nr1);
double zx11 = fabs(qx11.re);
double zy11 = fabs(qy11.re);
assert(wx11 <= rom->roq11.w0max && wx11 >= rom->roq11.w0min);
assert(zx11 <= rom->roq11.zmax && zx11 >= rom->roq11.zmin);
assert(wy11 <= rom->roq11.w0max && wy11 >= rom->roq11.w0min);
assert(zy11 <= rom->roq11.zmax && zy11 >= rom->roq11.zmin);
fill_unn_cache(&(ws11->ux_cache_11), acc_11_nr1_1->acc_n, acc_11_nr1_2->acc_n, true);
fill_unn_cache(&(ws11->uy_cache_11), acc_11_nr1_1->acc_m, acc_11_nr1_2->acc_m, true);
}
if (CALC_MR_KNM(mirror, 22) && rom->roq22.enabled){
double wx22 = w0_size(qx22, nr2);
double wy22 = w0_size(qy22, nr2);
double zx22 = fabs(qx22.re);
double zy22 = fabs(qy22.re);
assert(wx22 <= rom->roq22.w0max && wx22 >= rom->roq22.w0min);
assert(zx22 <= rom->roq22.zmax && zx22 >= rom->roq22.zmin);
assert(wy22 <= rom->roq22.w0max && wy22 >= rom->roq22.w0min);
assert(zy22 <= rom->roq22.zmax && zy22 >= rom->roq22.zmin);
fill_unn_cache(&(ws22->ux_cache_22), acc_22_nr2_1->acc_n, acc_22_nr2_2->acc_n, true);
fill_unn_cache(&(ws22->uy_cache_22), acc_22_nr2_1->acc_m, acc_22_nr2_2->acc_m, true);
}
// Iterate over each mode to calculate k_n1,m1,n2,m2
for (n = 0; n < inter.num_fields; n++) {
for (m = 0; m < inter.num_fields; m++) {
// If we are using the knm to calc knm then we do not need to bother
// with doing all the integrating for any of the lower triangle.
if (!calc_knm_transpose || (n >= m)) {
//Transform linear mode index system into actual TEM_NM mode
get_tem_modes_from_field_index(&n1, &m1, n);
get_tem_modes_from_field_index(&n2, &m2, m);
// compute some constants that are n1,n2,m1,m2 dependant
if (CALC_MR_KNM(mirror, 11) && rom->roq11.enabled){
complex_t znm1c1 = z_by_zc(z_by_z(acc_11_nr1_1->acc_n->prefac[n1], acc_11_nr1_1->acc_m->prefac[m1]),
z_by_z(acc_11_nr1_2->acc_n->prefac[n2], acc_11_nr1_2->acc_m->prefac[m2]));
mirror->knm_surf_motion_1o[l][n][m] = z_by_z(do_romhom_real_int(&rom->roq11, ws11->d_u_xy, &ws11->ux_cache_11, &ws11->uy_cache_11, n1, m1, n2, m2), znm1c1);
// the only difference between the incoming and outgoing computation
// is that q_1 = cminus(cconj(q_2)) which means that the gouy phase is
// opposite sign, or just the conjugate. As the knm matrix is hermitian
// we can compute everything at once from one integral computation.
// the msign(n1+n2) is the difference between the incoming beam
// coordinate system and that of the maps, the x-axis is inverted.
mirror->knm_surf_motion_1i[l][n][m] = cconj(z_by_x(mirror->knm_surf_motion_1o[l][n][m], msign(n1+n2)));
} else {
mirror->knm_surf_motion_1o[l][n][m] = (n1==n2 && m1==m2) ? complex_1 : complex_0;
mirror->knm_surf_motion_1i[l][n][m] = (n1==n2 && m1==m2) ? complex_1 : complex_0;
}
if (CALC_MR_KNM(mirror, 22) && rom->roq22.enabled) {
complex_t znm2c2 = z_by_zc(z_by_z(acc_22_nr2_1->acc_n->prefac[n1], acc_22_nr2_1->acc_m->prefac[m1]),
z_by_z(acc_22_nr2_2->acc_n->prefac[n2], acc_22_nr2_2->acc_m->prefac[m2]));
mirror->knm_surf_motion_2o[l][n][m] = z_by_z(do_romhom_real_int(&rom->roq22, ws22->d_u_xy, &ws22->ux_cache_22, &ws22->uy_cache_22, n1, m1, n2, m2), (znm2c2));
// the only difference between the incoming and outgoing computation
// is that q_1 = cminus(cconj(q_2)) which means that the gouy phase is
// opposite sign, or just the conjugate. As the knm matrix is hermitian
// we can compute everything at once from one integral computation.
// the msign(n1+n2) is the difference between the incoming beam
// coordinate system and that of the maps, the x-axis is inverted.
mirror->knm_surf_motion_2i[l][n][m] = cconj(z_by_x(mirror->knm_surf_motion_2o[l][n][m], msign(n1+n2)));
} else {
mirror->knm_surf_motion_2o[l][n][m] = (n1==n2 && m1==m2) ? complex_1 : complex_0;
mirror->knm_surf_motion_2i[l][n][m] = (n1==n2 && m1==m2) ? complex_1 : complex_0;
}
if (calc_knm_transpose && !(n1 == n2 && m1 == m2)) {
mirror->knm_surf_motion_1i[l][m][n] = cconj(mirror->knm_surf_motion_1i[l][n][m]);
mirror->knm_surf_motion_1o[l][m][n] = cconj(mirror->knm_surf_motion_1o[l][n][m]);
mirror->knm_surf_motion_2i[l][m][n] = cconj(mirror->knm_surf_motion_2i[l][n][m]);
mirror->knm_surf_motion_2o[l][m][n] = cconj(mirror->knm_surf_motion_2o[l][n][m]);
}
}
current_coeff++;
print_progress_and_time(num_coeffs, current_coeff, starttime);
}
}
}
endTimer(timer);
}
void calc_mirror_knm_surf_motions_map(mirror_t *m, double nr1, double nr2, complex_t qx1, complex_t qy1, complex_t qx2, complex_t qy2) {
assert(m!=NULL);
assert(m->num_surface_motions > 0);
assert(m->surface_motions != NULL);
......@@ -749,6 +895,9 @@ void mirror_calc_knm_surf_motions(mirror_t *m, double nr1, double nr2, complex_t
int current_coeff = 0;
for(k=0; k<m->num_surface_motions; k++){
if(m->surface_motions_isROM[k])
continue;
p.map = &inter.surface_motion_map_list[m->surface_motions[k]];
current_coeff = 1;
......@@ -762,15 +911,16 @@ void mirror_calc_knm_surf_motions(mirror_t *m, double nr1, double nr2, complex_t
for(out=in; out<inter.num_fields; out++){
get_tem_modes_from_field_index(&p.n2, &p.m2, out);
// Note that we do not apply the gouy phase here, which means that
// the result is already reversed gouy'd. For merging with other
// knm matrices for motion to field calculations later we need to
// apply the gouy
double const_fac = 1.0/(sqrt(pow(2, p.n1+p.n2+p.m1+p.m2-2) * fac(p.n1) * fac(p.n2) * fac(p.m1) * fac(p.m2)) * PI);
double const_fac = 1.0 / (sqrt(pow(2, p.n1+p.n2+p.m1+p.m2-2) * fac(p.n1) * fac(p.n2) * fac(p.m1) * fac(p.m2)) * PI);
double gx11 = gouy(qx1);
double gy11 = gouy(qy1);
double gx22 = gouy(qx2);
double gy22 = gouy(qy2);
double value = 0.0;
if(compute1){
if(compute1) {
double wx = w_size(qx1, nr1);
double wy = w_size(qy1, nr1);
p.constant = const_fac / (wx*wy);
......@@ -793,17 +943,16 @@ void mirror_calc_knm_surf_motions(mirror_t *m, double nr1, double nr2, complex_t
#else
integrate_riemann_sum_surf(&value, &p);
#endif
complex_t knm = z_by_xphr(complex_1, value, (p.n1 - p.n2)*gx11 + (p.m1 - p.m2)*gy11);
// the only difference between the incoming and outgoing computation
// is that q_1 = cminus(cconj(q_2)) which means that the gouy phase is
// opposite sign, or just the conjugate. As the knm matrix is hermitian
// we can compute everything at once from one integral computation.
// 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] = z_by_x(complex_1, value);
m->knm_surf_motion_1i[k][in][out] = cconj(z_by_x(complex_1, value* msign(p.n1+p.n2)));
//printf("%i%i %i%i %s\n",p.n1,p.m1,p.n2,p.m2, complex_form15(m->knm_surf_motion_1o[k][in][out]));
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)));
if(out!=in) {
m->knm_surf_motion_1o[k][out][in] = cconj(m->knm_surf_motion_1o[k][in][out]);
......@@ -839,11 +988,13 @@ void mirror_calc_knm_surf_motions(mirror_t *m, double nr1, double nr2, complex_t
// as the motion is opposite on n2 side from n1
value *= -1;
complex_t knm = z_by_xphr(complex_1, value, (p.n1 - p.n2)*gx22 + (p.m1 - p.m2)*gy22);
// use a minus sign (complex_m1) here because the surface motion will appear
// opposite from node 2 side, we also have the msign term here because the coord
// system of the outgoing beam has a negative x compared to the map coord system
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));
m->knm_surf_motion_2o[k][in][out] = z_by_x(knm, msign(p.n1+p.n2));
m->knm_surf_motion_2i[k][in][out] = cconj(knm);
if(out!=in) {
m->knm_surf_motion_2o[k][out][in] = cconj(m->knm_surf_motion_2o[k][in][out]);
......
......@@ -62,7 +62,8 @@ typedef struct knm_mirror_map_int_params {
void alloc_knm_accel_mirror_mem(long *bytes);
void mirror_calc_knm_surf_motions(mirror_t *m, double nr1, double nr2, complex_t qx1, complex_t qy1, complex_t qx2, complex_t qy2);
void calc_mirror_knm_surf_motions_map(mirror_t *m, double nr1, double nr2, complex_t qx1, complex_t qy1, complex_t qx2, complex_t qy2);
void calc_mirror_knm_surf_motions_rom(mirror_t *mirror, double nr1, double nr2, complex_t qx11, complex_t qy11, complex_t qx22, complex_t qy22, bitflag astigmatism);
void mirror_knm_matrix_copy(mirror_knm_t* src, mirror_knm_t* dest);
void mirror_knm_matrix_mult(mirror_knm_t* A, mirror_knm_t* B, mirror_knm_t *result);
......
......@@ -855,8 +855,8 @@ int allocate_memory(long int *bytes, long int *bytes1, long int *bytes2) {
return error_code;
}
if(mem.num_maps) {
inter.rom_maps = (rom_map_t*) calloc(mem.num_maps, sizeof(rom_map_t));
if(mem.num_maps || mem.num_surface_motion_maps) {
inter.rom_maps = (rom_map_t*) calloc(mem.num_maps + mem.num_surface_motion_maps, sizeof(rom_map_t));
}
if (mem.num_dof){
......
......@@ -6991,7 +6991,7 @@ void read_conf(const char *command_string) {
GET_INT
if (ival < 1 || ival > 2)
gerror("knm_change_q must be either 1 (Bayer-Helms) or 2 (Merged map)\n");
gerror("knm_change_q must be either 1 (Merged map) or 2 (Bayer-Helms)\n");
mirror->knm_change_q = ival;
} else if (strcmp(config_name, "knm_force_saved") == 0) {
......@@ -7068,15 +7068,35 @@ void read_surface_motion_map(const char *command_string) {
warn("Line `%s':\ntext '%s' ignored\n", command_string, rest_string);
}
surface_map_t *map = &inter.surface_motion_map_list[inter.num_surface_motion_maps];
const char *dot = strrchr(map_filename, '.');
bool isROMMap = false;
strcpy(map->filename, map_filename);
read_surface_map(map);
if(dot != NULL && dot != map_filename){
if(strncasecmp(dot, ".rom", 4) == 0) {
isROMMap = true;
}
}
if(map->type != PHASE_MAP)
gerror("Line `%s':\nSurface motion map file '%s' was not a phase map\n", command_string, map_filename);
if(isROMMap) {
rom_map_t *rom = &inter.rom_maps[inter.num_ROM_maps];
strcpy(rom->filename, map_filename);
read_rom_map(rom);
if(!rom->roq11.enabled || !rom->roq22.enabled)
gerror("Line `%s':\n"
"ROM file '%s' doesn't have either the 11 or 22 reflection weights.\n", command_string, map_filename);
} else {
surface_map_t *map = &inter.surface_motion_map_list[inter.num_surface_motion_maps];
strcpy(map->filename, map_filename);
read_surface_map(map);
if(map->type != PHASE_MAP)
gerror("Line `%s':\nSurface motion map file '%s' was not a phase map\n", command_string, map_filename);
}
int component_index = get_component_index_from_name(component_name);
if (component_index == NOT_FOUND) {
gerror("Line `%s':\nNo such component '%s'\n",
command_string, component_name);
......@@ -7091,8 +7111,14 @@ void read_surface_motion_map(const char *command_string) {
mirror_t *m = &inter.mirror_list[component_index];
if(m->num_surface_motions >= MAX_MAPS) gerror("Line `%s':\nMaximum number of surface motions reached for '%s'\n", command_string, component_name);
m->surface_motions[m->num_surface_motions] = inter.num_surface_motion_maps;
m->surface_motions_isROM[m->num_surface_motions] = isROMMap;
if(isROMMap){
m->surface_motions[m->num_surface_motions] = inter.num_ROM_maps;
} else {
m->surface_motions[m->num_surface_motions] = inter.num_surface_motion_maps;
}
if(strlen(transfer_func) > 0){
int i;
......@@ -7108,7 +7134,10 @@ void read_surface_motion_map(const char *command_string) {
} else
bug_error("unhandled");
inter.num_surface_motion_maps++;
if(isROMMap)
inter.num_ROM_maps++;
else
inter.num_surface_motion_maps++;
}
/*!
......
......@@ -1462,17 +1462,10 @@ void fill_optic_motion_coupling(){
fill_mirror_field_to_rot_motion(ROTY, m, z_by_x(G_Omega_ry, -factor_a_x), fcar, a1i_y, a1o_y, a2i_y, a2o_y, ac_1i, ac_1o, ac_2i, ac_2o);
} else if(m->motion_type[k] == SURFACE) {
// The coupling coefficients are generated without taking into account
// the map scaling as we want to know the coupling from a unit distance motion
// to various modes. Here though when we actually take into account how much
// the surface is actually moving we need to rescale the coupling factors.
double map_scale = inter.surface_motion_map_list[m->surface_motions[k]].scaling;
// There are 2 indicies effectively, the motion index(X,Rx, Ry, S0, S1,...) and surface motion index (S0, S1, S2,...))
// There are 2 indices effectively, the motion index(X, Rx, Ry, S0, S1,...) and surface motion index (S0, S1, S2,...))
size_t surf_idx = k - (m->num_motions - m->num_surface_motions);
complex_t H = z_by_x(calc_transfer_func(m->surface_motions_tf[surf_idx], s_fsig), factor_a_x * map_scale);
complex_t H = z_by_x(calc_transfer_func(m->surface_motions_tf[surf_idx], s_fsig), factor_a_x);
fill_surf_motion_coupling(surf_idx, m, H, factor_x_a, fcar, ac_1i, ac_1o, ac_2i, ac_2o, tuning1u, tuning1l, tuning2u, tuning2l);
}
......
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