kat_aa.c 239 KB
Newer Older
1
// $Id: kat_aa.c 3374 2011-09-12 16:29:43Z dzb610 $
2

cochrane's avatar
cochrane committed
3
/*!  
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
 * \file kat_aa.c 
 * \brief Routines for the automatic alignment of the interferometer
 * 
 * @section Copyright notice
 *
 *  This file is part of the interferometer simulation Finesse
 *  http://www.gwoptics.org/finesse
 *
 *  Copyright (C) 1999 onwards Andreas Freise
 *  with parts of the code written by Daniel Brown, Paul Cochrane
 *  and Gerhard Heinzel.
 *
 *  This program is free software; you can redistribute it and/or modify it under
 *  the terms of the GNU General Public License version 3 as published
 *  by the Free Software Foundation.
 *
 *  This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
 *  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
Daniel Brown's avatar
Daniel Brown committed
22
 *  See the GNU General Public License for more details.
23
 *
Daniel Brown's avatar
Daniel Brown committed
24
 *  You should have received a copy of the GNU General Public License along with
25 26
 *  this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place,
 *  Suite 330, Boston, MA 02111-1307 USA
27
 */
cochrane's avatar
cochrane committed
28

29 30
#include "kat.h"
#include "kat_inline.c"
31
#include "kat_io.h"
32
#include "kat_fortran.h"
33
#include "kat_optics.h"
34
#include "kat_aux.h"
35
#include "kat_dump.h"
36
#include "kat_aa.h"
37
#include "kat_check.h"
Daniel Brown's avatar
Daniel Brown committed
38
#include "kat_knm_mirror.h"
39
#include "kat_knm_bs.h"
40
#include "kat_knm_bayer.h"
41
#include "kat_knm_aperture.h"
42
#include "kat_calc.h"
43
#include "kat_knm_int.h"
44
#include <gsl/gsl_cblas.h>
45
#include <stdlib.h>
46

47
//#include <omp.h>
48

49
extern init_variables_t init;
50
extern interferometer_t inter;
51
extern options_t options;
52

53
extern local_var_t vlocal;
54

55
extern FILE *fp_log;
56

57 58
extern const complex_t complex_i; //!<  sqrt(-1) or 0 + i
extern const complex_t complex_1; //!<  1 but in complex space:  1 + 0i
59
extern const complex_t complex_m1; //!< -1 but in complex space: -1 + 0i
60
extern const complex_t complex_0; //!<  0 but in complex space:  0 + 0i
61

Daniel Brown's avatar
Daniel Brown committed
62 63
// define some temporary variable for computations
// these are extern definined in other files
Daniel Brown's avatar
Daniel Brown committed
64
bs_knm_t bstmp;
65

66
FILE * ipfile, * idfile;
67

68
const ABCD_t Unity = {1.0, 0.0, 0.0, 1.0}; //!< The identity matrix
69

70
int n_index; //!< node index; number of nodes in beam trace
71
bool defaultKNMChangeWarning = false;
72

73 74
zmatrix tmp_knm = NULL;

Daniel Brown's avatar
Daniel Brown committed
75 76 77 78 79
double mismatch(complex_t q1, complex_t q2) {
    /**
     * Given two beam parameters this computes a figure of merit for
     * how much they overlap, or the amount of mismatch between them.
     * 
Daniel Brown's avatar
Daniel Brown committed
80 81
     *    1 being no matching
     *    0 being perfect matching
Daniel Brown's avatar
Daniel Brown committed
82 83 84
     * 
     * Will return nan if the beam parameter imaginary part is zero or negative.
     * 
Daniel Brown's avatar
Daniel Brown committed
85
     * See ddb thesis equation 1.63 (this uses 1 - Eq.1.63 though)
Daniel Brown's avatar
Daniel Brown committed
86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
     * 
     * @param q1
     * @param q2
     * @return Mismatch figure of merit or nan if not calculable
     */
    double result;
    
    if(q1.im <= 0 || q2.im <= 0) {
        result = NAN;
    } else {
        complex_t t1 = z_m_z(cconj(q1), q2);
        double abs_t1 = zmod(t1);
        result = 4.0 * fabs(q1.im * q2.im) / (abs_t1 * abs_t1);
    }
    
Daniel Brown's avatar
Daniel Brown committed
101
    assert((result >= 0 && result <=(1+1e-15)) || isnan(result));
Daniel Brown's avatar
Daniel Brown committed
102
    
Daniel Brown's avatar
Daniel Brown committed
103
    return 1- result;
Daniel Brown's avatar
Daniel Brown committed
104 105
}

106 107 108
bool __print_mismatch(const char *dir, double mx, double my, char *output) {
    const char *fmt1 = "%s\t%s: %18.10g (x) %18.10g (y)\n";
    const char *fmt2 = "%s\t%s: %18.10g\n";
Daniel Brown's avatar
Daniel Brown committed
109 110
    
    if(!(isnan(mx) && isnan(my))){
111 112 113 114
        if(inter.mismatches_options & AVERAGE_Q){
            double am = (mx+my)/2;
            
            if(am >= inter.mismatches_lower){
115 116
                sprintf(output, fmt2, output, dir, am);
                return true;
117 118 119
            }
        } else {
            if(mx >= inter.mismatches_lower || my >= inter.mismatches_lower){
120 121
                sprintf(output, fmt1, output, dir, mx, my);
                return true;
122
            }
Daniel Brown's avatar
Daniel Brown committed
123 124
        }
    }
125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
    
    return false;
}

void __print_two_port_mismatches(const char *component_name, char *name, q_io_t *qm12, q_io_t *qm21, char *output){

    bool displayed = false;

    displayed |= __print_mismatch("1->2", mismatch(qm12->qxi, qm12->qxo),
                                          mismatch(qm12->qyi, qm12->qyo),
                                          output);

    displayed |= __print_mismatch("2->1", mismatch(qm21->qxi, qm21->qxo),
                                          mismatch(qm21->qyi, qm21->qyo),
                                          output);

    if(displayed) {
        message("%s: %s\n", component_name, name);
        message(output);
    }
Daniel Brown's avatar
Daniel Brown committed
145
}
146

Daniel Brown's avatar
Daniel Brown committed
147 148 149 150
/**
 * When called this will output the current mismatch status of the interferometer.
 */
void output_mismatch() {
Daniel Brown's avatar
Daniel Brown committed
151 152 153
    if(!inter.mismatches)
        return;
    
Daniel Brown's avatar
Daniel Brown committed
154 155 156 157
    // Idea here is to loop through each component and output the mode matching
    // for each beam direction within the component.
    
    int i = 0;
Daniel Brown's avatar
Daniel Brown committed
158
   
Daniel Brown's avatar
Daniel Brown committed
159
    message("\n---- Start interferometer mismatches ----\n");
Daniel Brown's avatar
Daniel Brown committed
160
    
161
    char output[1000] = {0};
Daniel Brown's avatar
Daniel Brown committed
162
    
Daniel Brown's avatar
Daniel Brown committed
163 164
    for(i=0; i<inter.num_mirrors; i++) {
        mirror_t *m = &inter.mirror_list[i];
165
        output[0] = '\0';
Daniel Brown's avatar
Daniel Brown committed
166
        
167
        bool print_refl  = !((inter.mismatches_options & IGNORE_LOW_R_MIRROR_BS) && m->R < inter.mismatches_lower);
Daniel Brown's avatar
Daniel Brown committed
168
        
169
        bool print_trans = !((inter.mismatches_options & IGNORE_LOW_T_MIRROR_BS) && m->T < inter.mismatches_lower);
Daniel Brown's avatar
Daniel Brown committed
170
        
171 172
        bool displayed = false;
        
173
        if(print_refl) {
174 175 176
            displayed |= __print_mismatch("1->1", mismatch(m->qm.qxt1_11, m->qm.qxt2_11),
                                                  mismatch(m->qm.qyt1_11, m->qm.qyt2_11),
                                                  output);
177 178 179
        }
        
        if(print_trans){
180 181 182
            displayed |= __print_mismatch("1->2", mismatch(m->qm.qxt1_12, m->qm.qxt2_12),
                                                  mismatch(m->qm.qyt1_12, m->qm.qyt2_12),
                                                  output);
183

184 185 186
            displayed |= __print_mismatch("2->1", mismatch(m->qm.qxt1_21, m->qm.qxt2_21),
                                                  mismatch(m->qm.qyt1_21, m->qm.qyt2_21),
                                                  output);
187
        }
Daniel Brown's avatar
Daniel Brown committed
188
        
189
        if(print_refl) {
190 191 192 193 194 195 196 197
            displayed |= __print_mismatch("2->2", mismatch(m->qm.qxt1_22, m->qm.qxt2_22),
                                                  mismatch(m->qm.qyt1_22, m->qm.qyt2_22),
                                                  output);
        }
        
        if(displayed) {
            message("mirror: %s\n", m->name);
            message(output);
198
        }
Daniel Brown's avatar
Daniel Brown committed
199 200 201 202
    }
    
    for(i=0; i<inter.num_beamsplitters; i++) {
        beamsplitter_t *bs = &inter.bs_list[i];
203
        output[0] = '\0';
Daniel Brown's avatar
Daniel Brown committed
204
        
205 206 207 208
        bool print_refl  = !((inter.mismatches_options & IGNORE_LOW_R_MIRROR_BS) && bs->R < inter.mismatches_lower);
        
        bool print_trans = !((inter.mismatches_options & IGNORE_LOW_T_MIRROR_BS) && bs->T < inter.mismatches_lower);
        
209
        bool displayed = false;
Daniel Brown's avatar
Daniel Brown committed
210
       
211 212 213 214 215 216 217 218
        if(print_refl) displayed |= __print_mismatch("1->2", mismatch(bs->qm.qxt1_12, bs->qm.qxt2_12),
                               mismatch(bs->qm.qyt1_12, bs->qm.qyt2_12), output);
        
        if(print_trans) displayed |= __print_mismatch("1->3", mismatch(bs->qm.qxt1_13, bs->qm.qxt2_13),
                               mismatch(bs->qm.qyt1_13, bs->qm.qyt2_13), output);
        
        if(print_refl) displayed |= __print_mismatch("2->1", mismatch(bs->qm.qxt1_21, bs->qm.qxt2_21),
                               mismatch(bs->qm.qyt1_21, bs->qm.qyt2_21), output);
Daniel Brown's avatar
Daniel Brown committed
219
        
220 221
        if(print_trans) displayed |= __print_mismatch("2->4", mismatch(bs->qm.qxt1_24, bs->qm.qxt2_24),
                               mismatch(bs->qm.qyt1_24, bs->qm.qyt2_24), output);
Daniel Brown's avatar
Daniel Brown committed
222
        
223 224
        if(print_trans) displayed |= __print_mismatch("3->1", mismatch(bs->qm.qxt1_31, bs->qm.qxt2_31),
                               mismatch(bs->qm.qyt1_31, bs->qm.qyt2_31), output);
Daniel Brown's avatar
Daniel Brown committed
225
        
226 227
        if(print_refl) displayed |= __print_mismatch("3->4", mismatch(bs->qm.qxt1_34, bs->qm.qxt2_34),
                               mismatch(bs->qm.qyt1_34, bs->qm.qyt2_34), output);
Daniel Brown's avatar
Daniel Brown committed
228
        
229 230
        if(print_trans) displayed |= __print_mismatch("4->2", mismatch(bs->qm.qxt1_42, bs->qm.qxt2_42),
                               mismatch(bs->qm.qyt1_42, bs->qm.qyt2_42), output);
Daniel Brown's avatar
Daniel Brown committed
231
        
232 233
        if(print_refl) displayed |= __print_mismatch("4->3", mismatch(bs->qm.qxt1_43, bs->qm.qxt2_43),
                               mismatch(bs->qm.qyt1_43, bs->qm.qyt2_43), output);
Daniel Brown's avatar
Daniel Brown committed
234
        
235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262
        if(displayed) {
            message("bs: %s\n", bs->name);
            message(output);
        }
    }
    
    for(i=0; i<inter.num_modulators; i++) {
        modulator_t *m = &inter.modulator_list[i];
        __print_two_port_mismatches("modulator",
                                    m->name,
                                    &m->qm12,
                                    &m->qm21,
                                    output);
    }
    
    for(i=0; i<inter.num_dbss; i++) {
        dbs_t *dbs = &inter.dbs_list[i];
        output[0] = '\0';
        
        bool displayed = false;

        displayed |= __print_mismatch("1->3", mismatch(dbs->qm13.qxi, dbs->qm13.qxo),
                                              mismatch(dbs->qm13.qyi, dbs->qm13.qyo),
                                              output);
        
        displayed |= __print_mismatch("2->1", mismatch(dbs->qm21.qxi, dbs->qm21.qxo),
                                              mismatch(dbs->qm21.qyi, dbs->qm21.qyo),
                                              output);
Daniel Brown's avatar
Daniel Brown committed
263
        
264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293
        displayed |= __print_mismatch("3->4", mismatch(dbs->qm34.qxi, dbs->qm34.qxo),
                                              mismatch(dbs->qm34.qyi, dbs->qm34.qyo),
                                              output);
        
        displayed |= __print_mismatch("4->2", mismatch(dbs->qm42.qxi, dbs->qm42.qxo),
                                              mismatch(dbs->qm42.qyi, dbs->qm42.qyo),
                                              output);
        
        if(displayed) {
            message("dbs: %s\n", dbs->name);
            message(output);
        }
    }
    
    for(i=0; i<inter.num_lenses; i++) {
        lens_t *m = &inter.lens_list[i];
        __print_two_port_mismatches("lens",
                                    m->name,
                                    &m->qm12,
                                    &m->qm21,
                                    output);
    }
    
    for(i=0; i<inter.num_diodes; i++) {
        diode_t *m = &inter.diode_list[i];
        __print_two_port_mismatches("diode",
                                    m->name,
                                    &m->qm12,
                                    &m->qm21,
                                    output);
Daniel Brown's avatar
Daniel Brown committed
294 295 296 297 298 299
    }
    
    message("\n---- End interferometer mismatches ----\n");
}


cochrane's avatar
cochrane committed
300
//! Compute coupling coefficients between two beams
301

302 303 304 305 306 307 308 309 310 311 312 313 314 315
/**
 * Computes matrix product:
 *  C = A*B
 * 
 * @param A
 * @param B
 * @param C
 */
void knm_matrix_mult(zmatrix A, zmatrix B, zmatrix C) {
    assert(A != NULL);
    assert(B != NULL);
    assert(C != NULL);

    zmatrix tmp = NULL;
316

317 318 319 320 321 322 323
    // if result != A | B then no need to use a temporary matrix to store result 
    // just bung it in the result matrix
    if ((C == A) || (C == B))
        tmp = tmp_knm;
    else
        tmp = C;

324 325
    size_t n, m, l;

326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342
    for (n = 0; n < mem.num_fields; n++) {
        for (m = 0; m < mem.num_fields; m++) {
            // initialise to 0 here just incase some other data is there
            tmp[n][m] = complex_0;

            for (l = 0; l < mem.num_fields; l++) {
                tmp[n][m] = z_pl_z(tmp[n][m], z_by_z(A[n][l], B[l][m]));
            }
        }
    }

    // if either A or B is the target matrix to store the result
    // we need to copy over the result
    if ((C == A) || (C == B))
        memcpy(C[0], tmp[0], mem.num_fields * mem.num_fields * sizeof (complex_t));
}

cochrane's avatar
cochrane committed
343 344 345 346 347 348 349 350 351 352 353
/*!
 * Each coefficient describes the coupling of the TEM_n1m1 mode into the
 * TEM_n2m2 mode as a function of the misalignments (gamma_x, gamma_y). In
 * TEM_nm, n gives the mode number in the horizontal plane and m the number
 * for the vertical modes. qx gives the beam parameter in the horizontal
 * plane, qy for the vertical modes. Please note that the notation of
 * gamma_x and gamma_y is not intuitive: gamma_x describes the rotation
 * around the vertical axis and gamma_y the rotation around a horizontal
 * axis through the optical component.

 * Utility functions called by the routines below:
354
 * - cconj(z) : return complex conjugate of z
355
 * - ceq(z1, z2) : check whether two complex numbers z1, z2 are equal
cochrane's avatar
cochrane committed
356 357 358 359 360 361
 * - z_by_z(z1, z2) : return multiplication of two complex numbers
 * - z_pl_z(z1, z2) : return sum of two complex numbers
 * - z_by_x(z, x) : multiply complex number by a real number
 * - inv_complex(z) : return the inverse of a complex number
 * - zsqrt(z) : complex square root of z
 * - zr(q) : compute Rayleigh range from Gaussian beam parameter q
cochrane's avatar
cochrane committed
362
 * - w0_size(q, nr) : compute beam waist size from beam parameter q and 
cochrane's avatar
cochrane committed
363
 *                   refrective index nr
cochrane's avatar
cochrane committed
364
 * - pow_complex(z, n1, n2) : return z**(n1/n2) with n1, n2 as integers
cochrane's avatar
cochrane committed
365 366
 * - msign(n) : return -1.0**n, with n as an integer
 * - fac(n) : return factorial of n
367 368 369 370 371 372 373 374 375
 *
 * \param n1 mode number of beam 1
 * \param m1 mode number of beam 1
 * \param n2 mode number of beam 2
 * \param m2 mode number of beam 2
 * \param qx1 Gaussian beam q parameter in x-plane for beam 1
 * \param qy1 Gaussian beam q parameter in y-plane for beam 1
 * \param qx2 Gaussian beam q parameter in x-plane for beam 2
 * \param qy2 Gaussian beam q parameter in y-plane for beam 2
adf's avatar
 
adf committed
376
 * \param gamma_x misalignment in x-plane 
377 378
 * \param gamma_y misalignment in y-plane
 * \param nr refractive index
379
 *
380
 * \see Test_k_mnmn()
cochrane's avatar
cochrane committed
381
 */
Daniel Brown's avatar
Daniel Brown committed
382
complex_t k_nmnm(int n1, int m1, int n2, int m2, complex_t qx1, complex_t qy1,
Daniel Brown's avatar
Daniel Brown committed
383
        complex_t qx2, complex_t qy2, double gamma_x, double gamma_y,
Daniel Brown's avatar
Daniel Brown committed
384
        double nr, const unsigned int knm_flag) {
Daniel Brown's avatar
Daniel Brown committed
385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409
    complex_t z, z1, z2, ztmp;

    // sanity checks on inputs
    // the mode numbers shouldn't be less than zero
    assert(n1 >= 0);
    assert(n2 >= 0);
    assert(m1 >= 0);
    assert(m2 >= 0);
    // the refractive index shouldn't be zero or less
    assert(nr > 0.0);

    // if the misalignment is zero and ...
    if ((gamma_x == 0) && (gamma_y == 0)) {
        // if the Gaussian beam parameters are equal 
        // (i.e. there is no mode mismatch) ...
        if (ceq(qx1, qx2) && ceq(qy1, qy2)) {
            // then return 1 for coupling into the same mode number and zero 
            // for all other couplings
            if (n1 == n2 && m1 == m2) {
                return (complex_1);
            } else {
                return (complex_0);
            }
        }
    }
410

Daniel Brown's avatar
Daniel Brown committed
411 412 413 414
    // check for inconsistent mode numbers
    if (n1 < 0 || n2 < 0 || m1 < 0 || m2 < 0) {
        bug_error("negative mode numbers");
    }
415

Daniel Brown's avatar
Daniel Brown committed
416
    // inter.knm is a global setting for chosing the type of algorithm that 
Daniel Brown's avatar
Daniel Brown committed
417 418 419 420 421 422 423
    // is used for computing the coupling coefficients (mostly for debugging)

    // The sign of gamma_y has to be different between the numeric algorithm
    // and the Bayer-Helms formula. This is probaly due to a wrong definition from
    // me when implementing angles.

    /*
Daniel Brown's avatar
Daniel Brown committed
424 425
        if ((inter.knm & 2 && !(gamma_x == 0.0) && !(gamma_y == 0.0)) ||
                (inter.knm & 4 && !(gamma_x == 0.0 && gamma_y == 0.0)) || inter.knm & 8) {
Daniel Brown's avatar
Daniel Brown committed
426
     */
Daniel Brown's avatar
Daniel Brown committed
427 428 429
    // added new knm debug for old routine
    if (knm_flag & INT_BAYER_HELMS_OLD) {
        if ((knm_flag & 1) && !options.quiet) { // TODO should be replaced by ENUM too?
Daniel Brown's avatar
Daniel Brown committed
430 431 432 433 434 435 436 437 438
            z1 = k_mm(n1, n2, qx1, qx2, gamma_x, nr);
            z2 = k_mm(m1, m2, qy1, qy2, -gamma_y, nr); // gamma + oder - ? 251004
            ztmp = z_by_z(z1, z2);

            message("%d %d %d %d | qx1 %s qy1 %s qx2 %s qy2 %s | "
                    "gammax %g, gammay %g) ",
                    n1, m1, n2, m2, complex_form(qx1), complex_form(qy1),
                    complex_form(qx2), complex_form(qy2), gamma_x, -gamma_y);
        }
439

Daniel Brown's avatar
Daniel Brown committed
440 441
        // do a numerical integration on a grid
        // with adaptive integration routine `dcuhre':
Daniel Brown's avatar
Daniel Brown committed
442
        z = faltung_adapt(n1, n2, m1, m2, qx1, qy1, qx2, qy2, gamma_x, gamma_y, nr, knm_flag);
Daniel Brown's avatar
Daniel Brown committed
443 444

        // compare Bayer-Helms with numericly derived coefficients
Daniel Brown's avatar
Daniel Brown committed
445
        if ((knm_flag & 1) && !options.quiet) {
Daniel Brown's avatar
Daniel Brown committed
446 447 448 449 450 451 452 453 454 455 456 457 458 459 460
            message("\nNumeric:     %s\nBayer-Helms: %s\n", complex_form(z),
                    complex_form(ztmp));
            if (capprox(z, ztmp)) {
                message("--> OK   ------\n");
            } else {
                message("--> DIFF ++++++\n");
            }
        }
    } else {
        // use coupling coefficients as in F. Bayer-Helms (see below)
        if (n1 == vlocal.kn1 && n2 == vlocal.kn2 && ceq(vlocal.kqx1, qx1)
                && ceq(vlocal.kqx2, qx2) && (gamma_x == vlocal.kgammax)
                && (nr == vlocal.knr)) {
            z1 = vlocal.zk1;
        } else {
461
            //omp_set_lock(&knmLock);
Daniel Brown's avatar
Daniel Brown committed
462 463 464 465 466 467 468 469
            z1 = k_mm(n1, n2, qx1, qx2, gamma_x, nr);
            vlocal.zk1 = z1;
            vlocal.kn1 = n1;
            vlocal.kn2 = n2;
            vlocal.kqx1 = qx1;
            vlocal.kqx2 = qx2;
            vlocal.kgammax = gamma_x;
            vlocal.knr = nr;
470
            //omp_unset_lock(&knmLock);
Daniel Brown's avatar
Daniel Brown committed
471
        }
472

Daniel Brown's avatar
Daniel Brown committed
473 474 475 476 477
        // -gammy_y, probably because of my definition of rotation as right-handed 
        // (around y axis)
        z2 = k_mm(m1, m2, qy1, qy2, -gamma_y, nr);
        z = z_by_z(z1, z2);

Daniel Brown's avatar
Daniel Brown committed
478
        if ((knm_flag & 1) && !options.quiet) {
Daniel Brown's avatar
Daniel Brown committed
479 480 481 482 483 484 485 486
            message("%d %d %d %d | qx1 %s qy1 %s qx2 %s qy2 %s | "
                    "gammax %g, gammay %g) ",
                    n1, m1, n2, m2,
                    complex_form(qx1), complex_form(qy1),
                    complex_form(qx2), complex_form(qy2), gamma_x, -gamma_y);
            message("\nBayer-Helms: %s\n", complex_form(z));
        }
    }
487

Daniel Brown's avatar
Daniel Brown committed
488 489 490 491 492 493 494 495 496 497 498 499
    // `reverse Gouy phase': I found out that on changing from one TEM 
    // base system to another I have to turn back the Gouy phase with 
    // respect to the old beam parameter and add the Gouy phase with 
    // respect to the new beam parameter. This is a required because of 
    // the way the coupling coefficients are defined here and the way 
    // FINESSE handles the Gouy phase: the Gouy phase is added explicitly 
    // to the complex amplitude coefficiants in every `space' whereas the 
    // coupling coefficients are derived using a formula in which the Gouy 
    // phase is implicitly given by spacial distribution.
    z = rev_gouy(z, n1, m1, n2, m2, qx1, qx2, qy1, qy2);

    return (z);
500 501
}

502
// same as above but without rev_gouy, to be used with new set_k_mirror function
503

Daniel Brown's avatar
Daniel Brown committed
504
complex_t k_nmnm_new(int n1, int m1, int n2, int m2, complex_t qx1, complex_t qy1,
505
        complex_t qx2, complex_t qy2, double gamma_x, double gamma_y, double nr, const unsigned int knm_flag) {
Daniel Brown's avatar
Daniel Brown committed
506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522
    complex_t z, z1, z2, ztmp;

    // sanity checks on inputs
    // the mode numbers shouldn't be less than zero
    assert(n1 >= 0);
    assert(n2 >= 0);
    assert(m1 >= 0);
    assert(m2 >= 0);
    // the refractive index shouldn't be zero or less
    assert(nr > 0.0);

    // if the misalignment is zero and ...
    if ((gamma_x == 0) && (gamma_y == 0)) {
        // if the Gaussian beam parameters are equal 
        // (i.e. there is no mode mismatch) ...
        if (ceq(qx1, qx2) && ceq(qy1, qy2)) {
            // then return 1 for coupling into the same mode number and zero 
523
            // for all other couplingsf
Daniel Brown's avatar
Daniel Brown committed
524 525 526 527 528 529 530
            if (n1 == n2 && m1 == m2) {
                return (complex_1);
            } else {
                return (complex_0);
            }
        }
    }
531

Daniel Brown's avatar
Daniel Brown committed
532 533 534 535
    // check for inconsistent mode numbers
    if (n1 < 0 || n2 < 0 || m1 < 0 || m2 < 0) {
        bug_error("negative mode numbers");
    }
536

Daniel Brown's avatar
Daniel Brown committed
537
    // inter.knm is a global setting for chosing the type of algorithm that 
Daniel Brown's avatar
Daniel Brown committed
538 539 540 541 542 543 544
    // is used for computing the coupling coefficients (mostly for debugging)

    // The sign of gamma_y has to be different between the numeric algorithm
    // and the Bayer-Helms formula. This is probaly due to a wrong definition from
    // me when implementing angles.

    /*
Daniel Brown's avatar
Daniel Brown committed
545 546
        if ((inter.knm & 2 && !(gamma_x == 0.0) && !(gamma_y == 0.0)) ||
                (inter.knm & 4 && !(gamma_x == 0.0 && gamma_y == 0.0)) || inter.knm & 8) {
Daniel Brown's avatar
Daniel Brown committed
547
     */
Daniel Brown's avatar
Daniel Brown committed
548
    // added new knm debug for old routine
549

Daniel Brown's avatar
Daniel Brown committed
550 551
    if (knm_flag & INT_BAYER_HELMS_OLD) {
        if ((knm_flag & 1) && !options.quiet) {
Daniel Brown's avatar
Daniel Brown committed
552 553 554 555 556 557 558 559 560
            z1 = k_mm(n1, n2, qx1, qx2, gamma_x, nr);
            z2 = k_mm(m1, m2, qy1, qy2, -gamma_y, nr); // gamma + oder - ? 251004
            ztmp = z_by_z(z1, z2);

            message("%d %d %d %d | qx1 %s qy1 %s qx2 %s qy2 %s | "
                    "gammax %g, gammay %g) ",
                    n1, m1, n2, m2, complex_form(qx1), complex_form(qy1),
                    complex_form(qx2), complex_form(qy2), gamma_x, -gamma_y);
        }
561

Daniel Brown's avatar
Daniel Brown committed
562 563
        // do a numerical integration on a grid
        // with adaptive integration routine `dcuhre':
Daniel Brown's avatar
Daniel Brown committed
564
        z = faltung_adapt(n1, n2, m1, m2, qx1, qy1, qx2, qy2, gamma_x, gamma_y, nr, knm_flag);
Daniel Brown's avatar
Daniel Brown committed
565 566

        // compare Bayer-Helms with numericly derived coefficients
Daniel Brown's avatar
Daniel Brown committed
567
        if ((knm_flag & 1) && !options.quiet) {
Daniel Brown's avatar
Daniel Brown committed
568 569 570 571 572 573 574 575 576
            message("\nNumeric:     %s\nBayer-Helms: %s\n", complex_form(z),
                    complex_form(ztmp));
            if (capprox(z, ztmp)) {
                message("--> OK   ------\n");
            } else {
                message("--> DIFF ++++++\n");
            }
        }
    } else {
577

Daniel Brown's avatar
Daniel Brown committed
578
        // use coupling coefficients as in F. Bayer-Helms (see below)
579
        // cache previously computed value for z1 if nothing has changed
580

Daniel Brown's avatar
Daniel Brown committed
581 582 583 584 585
        if (n1 == vlocal.kn1 && n2 == vlocal.kn2 && ceq(vlocal.kqx1, qx1)
                && ceq(vlocal.kqx2, qx2) && (gamma_x == vlocal.kgammax)
                && (nr == vlocal.knr)) {
            z1 = vlocal.zk1;
        } else {
586

Daniel Brown's avatar
Daniel Brown committed
587 588 589 590 591 592 593 594
            z1 = k_mm(n1, n2, qx1, qx2, gamma_x, nr);
            vlocal.zk1 = z1;
            vlocal.kn1 = n1;
            vlocal.kn2 = n2;
            vlocal.kqx1 = qx1;
            vlocal.kqx2 = qx2;
            vlocal.kgammax = gamma_x;
            vlocal.knr = nr;
595

Daniel Brown's avatar
Daniel Brown committed
596
        }
597

Daniel Brown's avatar
Daniel Brown committed
598 599 600 601 602
        // -gammy_y, probably because of my definition of rotation as right-handed 
        // (around y axis)
        z2 = k_mm(m1, m2, qy1, qy2, -gamma_y, nr);
        z = z_by_z(z1, z2);

Daniel Brown's avatar
Daniel Brown committed
603
        if ((knm_flag & 1) && !options.quiet) {
Daniel Brown's avatar
Daniel Brown committed
604 605 606 607 608 609 610 611
            message("%d %d %d %d | qx1 %s qy1 %s qx2 %s qy2 %s | "
                    "gammax %g, gammay %g) ",
                    n1, m1, n2, m2,
                    complex_form(qx1), complex_form(qy1),
                    complex_form(qx2), complex_form(qy2), gamma_x, -gamma_y);
            message("\nBayer-Helms: %s\n", complex_form(z));
        }
    }
612

Daniel Brown's avatar
Daniel Brown committed
613 614 615 616 617 618 619 620 621
    // `reverse Gouy phase': I found out that on changing from one TEM 
    // base system to another I have to turn back the Gouy phase with 
    // respect to the old beam parameter and add the Gouy phase with 
    // respect to the new beam parameter. This is a required because of 
    // the way the coupling coefficients are defined here and the way 
    // FINESSE handles the Gouy phase: the Gouy phase is added explicitly 
    // to the complex amplitude coefficiants in every `space' whereas the 
    // coupling coefficients are derived using a formula in which the Gouy 
    // phase is implicitly given by spacial distribution.
Daniel Brown's avatar
Daniel Brown committed
622
    // 22/02/2012 This will be moved to after the knm merging - Daniel
Daniel Brown's avatar
Daniel Brown committed
623 624 625
    //z = rev_gouy(z, n1, m1, n2, m2, qx1, qx2, qy1, qy2);

    return (z);
626 627
}

adf's avatar
 
adf committed
628

629

cochrane's avatar
cochrane committed
630
//! Coupling coefficients for TEM modes
631

cochrane's avatar
cochrane committed
632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650
/*!
 * F.  Bayer-Helms `Coupling coefficients of an incident wave and modes
 * of an spherical optical resonator in the case of mismatching and
 * misalignment', Appl. Optics Vol23 No9 1369-1380 (1984).
 * There is also a report (in German) which gives a more detailed
 * derivation: `Ableitung der Kopplunskoeffizienten fuer die Eigenfunktionen
 * eines sphaerischen optischen Resonators bei Fehlanpassung und 
 * Fehljustierung der einfallenden Welle' PTB Bericht PTB-Me-43 (1983)
 * 
 * k_mm ist for one degree of freedom (x or y) only. m1 and m2 are the 
 * mode numbers of the incoming and outgoing beam, q1 and q2 the 
 * respective Gaussian beam parameters, gamma is the misalignment angle 
 * between the beams(!) and nr the index of refraction of the medium in 
 * which the coupling happens. The latter is needed to compute the diffraction 
 * angle from the Gaussian beam parameter.
 * 
 * The function calls S_u and S_g (see below). The names of parameters and
 * variables are chosen such that they are similar to the notation used in the
 * paper.
651 652 653 654 655 656 657
 *
 * \param m1 mode number for incoming beam
 * \param m2 mode number for outgoing beam
 * \param q1 Gaussian beam parameter for incoming beam
 * \param q2 Gaussian beam parameter for outgoing beam
 * \param gamma misalignment angle between the beams
 * \param nr refractive index
cochrane's avatar
cochrane committed
658
 */
659
complex_t k_mm(int m1, int m2, complex_t q1, complex_t q2,
Daniel Brown's avatar
Daniel Brown committed
660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692
        double gamma, double nr) {
    double Theta, K0;
    double gg0, z2z0;
    double zr;
    double out1;
    complex_t K;
    complex_t XS, X, FS, F;
    complex_t ex;
    complex_t k1, x1, x2;
    complex_t out2, out3;

    // sanity checks on inputs
    // mode numbers bounds; they should not be less than zero
    assert(m1 >= 0);
    assert(m2 >= 0);

    // refractive index should be > 0 
    assert(nr > 0.0);

    if (ceq(q1, q2)) {
        K = complex_0;
        K0 = 0.0;
        k1 = complex_1;
        FS = complex_0;
        F = complex_0;
    } else {
        zr = 1.0 / q2.im;
        K.re = 0.5 * (q1.im - q2.im) * zr;
        K.im = 0.5 * (q1.re - q2.re) * zr;
        K0 = 2 * K.re;
        k1 = inv_complex(zsqrt(z_pl_z(complex_1, cconj(K))));
        FS = z_by_x(K, 1.0 / (2 * (1.0 + K0)));
        F = z_by_x(cconj(K), 0.5);
693
    }
694

Daniel Brown's avatar
Daniel Brown committed
695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722
    // I think this theta variable is the denominator that comes from bayer-helms equation 12
    Theta = w0_size(q2, nr) / z_r(q2);

    //printf("\nw0=%g\n",w0_size(q2, nr) );
    //printf("zr=%g\n",z_r(q2));

    if (gamma == 0) {
        x1 = x2 = XS = X = ex = complex_0;
        gg0 = 0.0;
    } else {
        //    gg0 = -1.0 * gamma / Theta; // sin (gamma) ???

        //See equation 12, gamma/gamma^hat_0
        gg0 = gamma / Theta; // changed sign of gamma 010708

        //See equation 19 bayer-helms, z^hat_2 / z^hat_0
        z2z0 = z_q(q2) / z_r(q2);

        x1 = z_m_z(co(z2z0, 0), complex_i); // inner bracket of X^hat eq 19
        x2 = z_pl_z(co(z2z0, 0), ci(z_pl_z(complex_1, z_by_x(cconj(K), 2)))); //inner bracket of X eq 19

        // multiply by k1 which is (1+K*)^-0.5
        // also multiply inner brackets calculated above by gamma/gamma^hat_0 factor
        // TODO it appears as though we are missing minus signs according to equation 19 bayer-helms?
        XS = z_by_z(k1, z_by_x(x1, gg0)); //X^hat eq 19
        X = z_by_z(k1, z_by_x(x2, gg0)); //X eq 19
        ex = z_by_x(z_by_z(XS, X), 0.5);
    }
723

Daniel Brown's avatar
Daniel Brown committed
724 725 726
    out1 = msign(m2);
    out2 = z_by_x(pow_complex(z_pl_z(complex_1, cconj(K)), -m1 - m2 - 1, 2),
            sqrt(fac(m1) * fac(m2) * pow(1.0 + K0, m1 + 0.5)));
727

Daniel Brown's avatar
Daniel Brown committed
728 729
    complex_t sg = S_g(X, XS, F, FS, m1, m2);
    complex_t su = S_u(X, XS, F, FS, m1, m2);
730

Daniel Brown's avatar
Daniel Brown committed
731
    out3 = z_by_xphr(z_by_z(out2, z_m_z(sg, su)), out1 * exp(-ex.re), -ex.im);
732

Daniel Brown's avatar
Daniel Brown committed
733
    return (out3);
734 735
}

adf's avatar
 
adf committed
736
//! 
737

738
/*!
adf's avatar
 
adf committed
739
 * 
740 741 742 743 744 745
 * \param X ???
 * \param XS ???
 * \param F ???
 * \param FS ???
 * \param m1 ???
 * \param m2 ???
746 747
 *
 * \see Test_S_u()
748
 */
749
complex_t S_u(complex_t X, complex_t XS, complex_t F, complex_t FS,
Daniel Brown's avatar
Daniel Brown committed
750 751 752 753 754 755 756 757 758
        int m1, int m2) {
    int s3, s4;
    complex_t nom1, nom2;
    double denom1, denom2;
    complex_t z, z1, z2;

    if (m1 == 0 || m2 == 0) {
        return (complex_0);
    }
759

Daniel Brown's avatar
Daniel Brown committed
760 761 762 763 764 765 766 767 768
    if (m1 == 1 || m1 == 2) {
        s3 = 0;
    } else {
        if (odd(m1 - 1)) {
            s3 = (int) ((m1 - 2) / 2);
        } else {
            s3 = (int) ((m1 - 1) / 2);
        }
    }
769

Daniel Brown's avatar
Daniel Brown committed
770 771 772 773 774 775 776 777 778
    if (m2 == 1 || m2 == 2) {
        s4 = 0;
    } else {
        if (odd(m2 - 1)) {
            s4 = (int) ((m2 - 2) / 2);
        } else {
            s4 = (int) ((m2 - 1) / 2);
        }
    }
779

Daniel Brown's avatar
Daniel Brown committed
780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797
    z = complex_0;
    int i, j;
    for (i = 0; i <= s3; i++) {
        for (j = 0; j <= s4; j++) {
            nom1 = z_by_z(pow_complex(XS, m1 - 2 * i - 1, 1),
                    pow_complex(X, m2 - 2 * j - 1, 1));
            denom1 = fac(m1 - 2 * i - 1) * fac(m2 - 2 * j - 1);
            z1 = z_by_x(nom1, msign(i) / denom1);
            z2 = complex_0;
            int k;
            for (k = 0; k <= nmin(i, j); k++) {
                nom2 = z_by_z(pow_complex(FS, i - k, 1), pow_complex(F, j - k, 1));
                denom2 = fac(2 * k + 1) * fac(i - k) * fac(j - k);
                z2 = z_pl_z(z2, z_by_x(nom2, msign(k) / denom2));
            }
            z = z_pl_z(z, z_by_z(z1, z2));
        }
    }
798

Daniel Brown's avatar
Daniel Brown committed
799
    return (z);
800 801
}

adf's avatar
 
adf committed
802
//! 
803

804
/*! 
adf's avatar
 
adf committed
805
 * 
806 807 808 809 810 811
 * \param X ???
 * \param XS ???
 * \param F ???
 * \param FS ???
 * \param m1 ???
 * \param m2 ???
812 813
 *
 * \see Test_S_g()
814
 */
815
complex_t S_g(complex_t X, complex_t XS, complex_t F, complex_t FS,
Daniel Brown's avatar
Daniel Brown committed
816 817 818 819 820 821 822 823 824 825 826 827 828 829 830
        int m1, int m2) {
    int s1, s2;
    complex_t z, z1, z2;
    complex_t nom1, nom2;
    double denom1, denom2;

    if (m1 == 0 || m1 == 1) {
        s1 = 0;
    } else {
        if (odd(m1)) {
            s1 = (int) ((m1 - 1) / 2);
        } else {
            s1 = (int) (m1 / 2);
        }
    }
831

Daniel Brown's avatar
Daniel Brown committed
832 833 834 835 836 837 838 839 840
    if (m2 == 0 || m2 == 1) {
        s2 = 0;
    } else {
        if (odd(m2)) {
            s2 = (int) ((m2 - 1) / 2);
        } else {
            s2 = (int) (m2 / 2);
        }
    }
841

Daniel Brown's avatar
Daniel Brown committed
842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859
    z = complex_0;
    int i, j;
    for (i = 0; i <= s1; i++) {
        for (j = 0; j <= s2; j++) {
            nom1 =
                    z_by_z(pow_complex(XS, m1 - 2 * i, 1), pow_complex(X, m2 - 2 * j, 1));
            denom1 = fac(m1 - 2 * i) * fac(m2 - 2 * j);
            z1 = z_by_x(nom1, msign(i) / denom1);
            z2 = complex_0;
            int k;
            for (k = 0; k <= nmin(i, j); k++) {
                nom2 = z_by_z(pow_complex(FS, i - k, 1), pow_complex(F, j - k, 1));
                denom2 = fac(2 * k) * fac(i - k) * fac(j - k);
                z2 = z_pl_z(z2, z_by_x(nom2, msign(k) / denom2));
            }
            z = z_pl_z(z, z_by_z(z1, z2));
        }
    }
860

Daniel Brown's avatar
Daniel Brown committed
861
    return (z);
862 863
}

864 865
void __print_coupling_info_line(const char* name, coupling_info_t *c) {
    if (inter.debug)
866
        message("%10.10s %5.5s  max=%i \n", name, c->coupling_off ? "off" : "on", c->max_coupling_order);
867 868
}

869
void __set_no_reduce_coupling_info(coupling_info_t *c) {
870 871 872 873 874 875 876 877
    c->coupling_off = false;
    c->has_xeven = true;
    c->has_yeven = true;
    c->has_xodd = true;
    c->has_yodd = true;
    c->max_coupling_order = inter.tem;
}

878 879 880 881 882 883 884 885
/**
 * Iterates through each objects and sets up the
 * coupling_info_t objects. Must be called before
 * matrices are built as this information is used to
 * determine which elements to allocate.
 */
void set_coupling_info() {
    int i;
886

887
    int max = inter.tem;
888 889

    if (inter.debug && options.use_coupling_reduction) {
890 891 892 893
        message("-----------------------------------\n");
        message("Reduced order coupling information\n");
        message("-----------------------------------\n");
    }
894 895

    for (i = 0; i < inter.num_mirrors; i++) {
896
        mirror_t *mirror = &inter.mirror_list[i];
897 898 899 900

        coupling_info_t * v[4] = {&mirror->a_cplng_11, &mirror->a_cplng_22,
            &mirror->a_cplng_12, &mirror->a_cplng_21};

901
        int j;
902 903 904

        for (j = 0; j < 4; j++) {
            if (!options.use_coupling_reduction) {
905 906 907 908
                __set_no_reduce_coupling_info(v[j]);
            } else {
                // switch all coupling off to begin with
                v[j]->coupling_off = true;
Daniel Brown's avatar
Daniel Brown committed
909 910 911 912
                v[j]->has_xeven = true;
                v[j]->has_yeven = true;
                v[j]->has_xodd = true;
                v[j]->has_yodd = true;
913
                v[j]->max_coupling_order = max;
914

915
                // if there is any mismatch turn coupling on
916
                if (mirror->mismatching & pow_two(j)) {
917 918 919
                    v[j]->coupling_off = false;
                }

920
                if (mirror->beta_x != 0.0 || is_var_tuned(&mirror->beta_x)) {
921 922 923
                    v[j]->coupling_off = false;
                }

924
                if (mirror->beta_y != 0.0 || is_var_tuned(&mirror->beta_y)) {
925 926 927 928 929 930 931
                    v[j]->coupling_off = false;
                }

                // if there are maps then ensure coupling is on
                // could probably be more picky about which type of 
                // map triggers this, e.g. a transmission map doesn't
                // affect reflection.
932
                if (mirror->num_maps > 0) {
933 934
                    v[j]->coupling_off = false;
                }
935 936

                if (mirror->map_rom != NULL) {
937 938
                    v[j]->coupling_off = false;
                }
939 940 941

                __print_coupling_info_line(mirror->name, v[j]);
            }
942 943
        }
    }
944 945

    for (i = 0; i < inter.num_beamsplitters; i++) {
946
        beamsplitter_t *bs = &inter.bs_list[i];
947 948 949 950 951 952 953 954

        coupling_info_t * v[12] = {&bs->a_cplng_12, &bs->a_cplng_21,
            &bs->a_cplng_34, &bs->a_cplng_43,
            &bs->a_cplng_13, &bs->a_cplng_31,
            &bs->a_cplng_24, &bs->a_cplng_42,
            &bs->a_cplng_11, &bs->a_cplng_22,
            &bs->a_cplng_33, &bs->a_cplng_44};

955
        int j;
956 957 958

        for (j = 0; j < 12; j++) {
            if (!options.use_coupling_reduction) {
959 960 961 962
                __set_no_reduce_coupling_info(v[j]);
            } else {
                // switch all coupling off to begin with
                v[j]->coupling_off = true;
Daniel Brown's avatar
Daniel Brown committed
963 964
                v[j]->has_xeven = true;
                v[j]->has_yeven = true;
965 966
                v[j]->has_xodd = true;
                v[j]->has_yodd = true;
967 968 969
                v[j]->max_coupling_order = max;

                // if there is any mismatch turn coupling on
970
                if (bs->mismatching & pow_two(j)) {
971 972 973
                    v[j]->coupling_off = false;
                }

974
                if (bs->beta_x != 0.0 || is_var_tuned(&bs->beta_x)) {
975 976 977
                    v[j]->coupling_off = false;
                }

978
                if (bs->beta_y != 0.0 || is_var_tuned(&bs->beta_y)) {
979 980 981 982 983 984 985
                    v[j]->coupling_off = false;
                }

                // if there are maps then ensure coupling is on
                // could probably be more picky about which type of 
                // map triggers this, e.g. a transmission map doesn't
                // affect reflection.
986
                if (bs->num_maps > 0) {
987 988 989 990 991
                    v[j]->coupling_off = false;
                }

                __print_coupling_info_line(bs->name, v[j]);
            }
992 993
        }
    }
994 995

    for (i = 0; i < inter.num_spaces; i++) {
996
        space_t *space = &inter.space_list[i];
997

998 999 1000
        // These should be order so that they
        // represent the same mode mismatch ordering
        // in the bitflag space->mismatching
1001 1002
        coupling_info_t * v[2] = {&space->a_cplng_12, &space->a_cplng_21};

1003
        int j;
1004

1005 1006
        // check if each direction has any mode mismatch
        // as that is the only coupling possible for a space
1007 1008
        for (j = 0; j < 2; j++) {
            if (!options.use_coupling_reduction) {
1009 1010 1011 1012 1013
                __set_no_reduce_coupling_info(v[j]);
            } else {
                v[j]->coupling_off = !(space->mismatching & pow_two(j));
                v[j]->has_xeven = true;
                v[j]->has_yeven = true;
1014 1015
                v[j]->has_xodd = true;
                v[j]->has_yodd = true;