kat_fortran.c 13.1 KB
Newer Older
1 2
// $Id$

3
/*!
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
 * \file kat_fortran.c
 * \brief Interfaces to legacy fortran routines.
 * 
 * @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.
22
 *  See the GNU General Public License for more details.
23
 *
24
 *  You should have received a copy of the GNU General Public License along with
25 26 27 28
 *  this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place,
 *  Suite 330, Boston, MA 02111-1307 USA
 */

29

30 31
#include "kat.h"
#include "kat_inline.c"
32
#include "kat_fortran.h"
33
#include "kat_optics.h"
34
#include "kat_io.h"
35
#include "kat_aux.h"
36

37
extern local_var_t vlocal;
38

39
extern interferometer_t inter;
40
extern options_t options;
41

42
//! Integration routine ???  Check and expand
43 44 45
/*!
 *
 */
46
extern void dcuhre_(int *, int *, double *, double *, int *, int *,
Daniel Brown's avatar
Daniel Brown committed
47 48 49
        void (*funcs) (int *, double *, int *, double *),
        double *, double *, int *, int *, int *, double *, double *,
        int *, int *, double *);
50

51
//! ???
Daniel Brown's avatar
Daniel Brown committed
52

53 54 55 56 57 58 59 60 61 62 63 64 65
/*!
 * \param n1 mode of input beam
 * \param n2 mode of output beam ???
 * \param m1 mode of input beam
 * \param m2 mode of output beam ???
 * \param qx1 Gaussian beam parameter in x-plane for input beam
 * \param qy1 Gaussian beam parameter in y-plane for input beam
 * \param qx2 Gaussian beam parameter in x-plane for output beam ???
 * \param qy2 Gaussian beam parameter in y-plane for output beam ???
 * \param gamma_x misalignment angle in x-plane
 * \param gamma_y misalignment angle in y-plane
 * \param nr refractive index
 */
66
complex_t faltung_adapt(int n1, int n2, int m1, int m2,
Daniel Brown's avatar
Daniel Brown committed
67 68
        complex_t qx1, complex_t qy1,
        complex_t qx2, complex_t qy2,
Daniel Brown's avatar
Daniel Brown committed
69
        double gamma_x, double gamma_y, double nr, const unsigned int knm_flags) {
Daniel Brown's avatar
Daniel Brown committed
70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143
    double r0x, r0y;
    double w1, w2;
    complex_t z;
    double start[3] = {0};
    double stop[3] = {0};
    double result[3] = {0};
    double error[3] = {0.0};
    complex_t printerr;
    int numofcalls;
    int ifail;
    int dim1, dim2, minop, maxop;
    double relerr, abserr;
    int worklen, key, restar;
    void (*func_ptr) (int *, double *, int *, double *);

    //  error = {0.0};
    numofcalls = 0;
    ifail = 0;

    func_ptr = function_;

    vlocal.lnr = nr;
    vlocal.ln1 = n1;
    vlocal.ln2 = n2;
    vlocal.lm1 = m1;
    vlocal.lm2 = m2;

    vlocal.lqx1 = qx1;
    vlocal.lqx2 = qx2;
    vlocal.lqy1 = qy1;
    vlocal.lqy2 = qy2;

    vlocal.lcx = cos(gamma_x);
    vlocal.lcy = cos(gamma_y);
    vlocal.lsx = sin(gamma_x);
    vlocal.lsy = sin(gamma_y);

    // breite der n-ten mode ist ca. 2sqrt(n) w
    w1 = w_size(qx1, nr);
    w2 = w_size(qx2, nr);
    r0x = 5 * max(sqrt(n1 + 0.5) * w1, sqrt(n2 + 0.5) * w2);

    w1 = w_size(qy1, nr);
    w2 = w_size(qy2, nr);
    r0y = 5 * max(sqrt(m1 + 0.5) * w1, sqrt(m2 + 0.5) * w2);

    start[0] = -r0x;
    start[1] = -r0y;

    stop[0] = -start[0];
    stop[1] = -start[1];

    dim1 = 2;
    dim2 = 2;
    minop = MIN_INT_OP;
    maxop = init.maxintop;
    relerr = init.relerr;
    abserr = init.abserr;
    worklen = mem.worklen;
    key = 0;
    restar = 0;

    numofcalls = ifail = 0;

    dcuhre_(&dim1, &dim2, start, stop, &minop, &maxop, func_ptr,
            &abserr, &relerr, &key, &worklen, &restar, result, error,
            &numofcalls, &ifail, mem.work);

    z.re = result[0];
    z.im = result[1];

    printerr.re = error[0];
    printerr.im = error[1];

Daniel Brown's avatar
Daniel Brown committed
144
    if ((knm_flags & 1) && !options.quiet) {
Daniel Brown's avatar
Daniel Brown committed
145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 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 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
        message("(error %s, calls %d, ifail %d)", complex_form15(printerr), numofcalls, ifail);
    }

    /* 
       SUBROUTINE DCUHRE(NDIM, NUMFUN, A, B, MINPTS, MAXPTS, FUNSUB, EPSABS,
       +                  EPSREL, KEY, NW, RESTAR, RESULT, ABSERR, NEVAL, IFAIL,
       +                  WORK)
       C     NDIM   Integer.
       C            Number of variables. 1 < NDIM <=  15.
       C     NUMFUN Integer.
       C            Number of components of the integral.
       C     A      Real array of dimension NDIM.
       C            Lower limits of integration.
       C     B      Real array of dimension NDIM.
       C            Upper limits of integration.
       C     MINPTS Integer.
       C            Minimum number of function evaluations.
       C     MAXPTS Integer.
       C            Maximum number of function evaluations.
       C            The number of function evaluations over each subregion
       C            is NUM.
       C            If (KEY = 0 or KEY = 1) and (NDIM = 2) Then
       C              NUM = 65
       C            Elseif (KEY = 0 or KEY = 2) and (NDIM = 3) Then
       C              NUM = 127
       C            Elseif (KEY = 0 and NDIM > 3) or (KEY = 3) Then
       C              NUM = 1 + 4*2*NDIM + 2*NDIM*(NDIM-1) + 4*NDIM*(NDIM-1) +
       C                    4*NDIM*(NDIM-1)*(NDIM-2)/3 + 2**NDIM
       C            Elseif (KEY = 4) Then
       C              NUM = 1 + 3*2*NDIM + 2*NDIM*(NDIM-1) + 2**NDIM
       C            MAXPTS >= 3*NUM and MAXPTS >= MINPTS
       C            For 3 < NDIM < 13 the minimum values for MAXPTS are:
       C             NDIM =    4   5   6    7    8    9    10   11    12
       C            KEY = 3:  459 819 1359 2151 3315 5067 7815 12351 20235
       C            KEY = 4:  195 309  483  765 1251 2133 3795  7005 13299
       C     FUNSUB Externally declared subroutine for computing
       C            all components of the integrand at the given
       C            evaluation point.
       C            It must have parameters (NDIM, X, NUMFUN, FUNVLS)
       C            Input parameters:
       C              NDIM   Integer that defines the dimension of the
       C                     integral.
       C              X      Real array of dimension NDIM
       C                     that defines the evaluation point.
       C              NUMFUN Integer that defines the number of
       C                     components of I.
       C            Output parameter:
       C              FUNVLS Real array of dimension NUMFUN
       C                     that defines NUMFUN components of the integrand.
       C
       C     EPSABS Real.
       C            Requested absolute error.
       C     EPSREL Real.
       C            Requested relative error.
       C     KEY    Integer.
       C            Key to selected local integration rule.
       C            KEY = 0 is the default value.
       C                  For NDIM = 2 the degree 13 rule is selected.
       C                  For NDIM = 3 the degree 11 rule is selected.
       C                  For NDIM > 3 the degree  9 rule is selected.
       C            KEY = 1 gives the user the 2 dimensional degree 13
       C                  integration rule that uses 65 evaluation points.
       C            KEY = 2 gives the user the 3 dimensional degree 11
       C                  integration rule that uses 127 evaluation points.
       C            KEY = 3 gives the user the degree 9 integration rule.
       C            KEY = 4 gives the user the degree 7 integration rule.
       C                  This is the recommended rule for problems that
       C                  require great adaptivity.
       C     NW     Integer.
       C            Defines the length of the working array WORK.
       C            Let MAXSUB denote the maximum allowed number of subregions
       C            for the given values of MAXPTS, KEY and NDIM.
       C            MAXSUB = (MAXPTS-NUM)/(2*NUM) + 1
       C            NW should be greater or equal to
       C            MAXSUB*(2*NDIM+2*NUMFUN+2) + 17*NUMFUN + 1
       C            For efficient execution on parallel computers
       C            NW should be chosen greater or equal to
       C            MAXSUB*(2*NDIM+2*NUMFUN+2) + 17*NUMFUN*MDIV + 1
       C            where MDIV is the number of subregions that are divided in
       C            each subdivision step.
       C            MDIV is default set internally in DCUHRE equal to 1.
       C            For efficient execution on parallel computers
       C            with NPROC processors MDIV should be set equal to
       C            the smallest integer such that MOD(2*MDIV, NPROC) = 0.
       C
       C     RESTAR Integer.
       C            If RESTAR = 0, this is the first attempt to compute
       C            the integral.
       C            If RESTAR = 1, then we restart a previous attempt.
       C            In this case the only parameters for DCUHRE that may
       C            be changed (with respect to the previous call of DCUHRE)
       C            are MINPTS, MAXPTS, EPSABS, EPSREL and RESTAR.
       C
       C   ON RETURN
       C
       C     RESULT Real array of dimension NUMFUN.
       C            Approximations to all components of the integral.
       C     ABSERR Real array of dimension NUMFUN.
       C            Estimates of absolute errors.
       C     NEVAL  Integer.
       C            Number of function evaluations used by DCUHRE.
       C     IFAIL  Integer.
       C            IFAIL = 0 for normal exit, when ABSERR(K) <=  EPSABS or
       C              ABSERR(K) <=  ABS(RESULT(K))*EPSREL with MAXPTS or less
       C              function evaluations for all values of K,
       C              1 <= K <= NUMFUN .
       C            IFAIL = 1 if MAXPTS was too small for DCUHRE
       C              to obtain the required accuracy. In this case DCUHRE
       C              returns values of RESULT with estimated absolute
       C              errors ABSERR.
       C            IFAIL = 2 if KEY is less than 0 or KEY greater than 4.
       C            IFAIL = 3 if NDIM is less than 2 or NDIM greater than 15.
       C            IFAIL = 4 if KEY = 1 and NDIM not equal to 2.
       C            IFAIL = 5 if KEY = 2 and NDIM not equal to 3.
       C            IFAIL = 6 if NUMFUN is less than 1.
       C            IFAIL = 7 if volume of region of integration is zero.
       C            IFAIL = 8 if MAXPTS is less than 3*NUM.
       C            IFAIL = 9 if MAXPTS is less than MINPTS.
       C            IFAIL = 10 if EPSABS < 0 and EPSREL < 0.
       C            IFAIL = 11 if NW is too small.
       C            IFAIL = 12 if unlegal RESTAR.
       C     WORK   Real array of dimension NW.
       C            Used as working storage.
       C            WORK(NW) = NSUB, the number of subregions in the data
       C            structure.
       C            Let WRKSUB=(NW-1-17*NUMFUN*MDIV)/(2*NDIM+2*NUMFUN+2)
       C            WORK(1), ..., WORK(NUMFUN*WRKSUB) contain
       C              the estimated components of the integrals over the
       C              subregions.
       C            WORK(NUMFUN*WRKSUB+1), ..., WORK(2*NUMFUN*WRKSUB) contain
       C              the estimated errors over the subregions.
       C            WORK(2*NUMFUN*WRKSUB+1), ..., WORK(2*NUMFUN*WRKSUB+NDIM*
       C              WRKSUB) contain the centers of the subregions.
       C            WORK(2*NUMFUN*WRKSUB+NDIM*WRKSUB+1), ..., WORK((2*NUMFUN+
       C              NDIM)*WRKSUB+NDIM*WRKSUB) contain subregion half widths.
       C            WORK(2*NUMFUN*WRKSUB+2*NDIM*WRKSUB+1), ..., WORK(2*NUMFUN*
       C              WRKSUB+2*NDIM*WRKSUB+WRKSUB) contain the greatest errors
       C              in each subregion.
       C            WORK((2*NUMFUN+2*NDIM+1)*WRKSUB+1), ..., WORK((2*NUMFUN+
       C              2*NDIM+1)*WRKSUB+WRKSUB) contain the direction of
       C              subdivision in each subregion.
       C            WORK(2*(NDIM+NUMFUN+1)*WRKSUB), ..., WORK(2*(NDIM+NUMFUN+1)*
       C              WRKSUB+ 17*MDIV*NUMFUN) is used as temporary
       C              storage in DADHRE.
       C
     */

    return (z);
293 294
}

295 296 297



298
//! Function to be integrated by 'faltung_adapt': u_{n1 m2} times u_{n2 m2}^*
Daniel Brown's avatar
Daniel Brown committed
299

300 301 302 303 304 305
/*!
 * \param din ???
 * \param point ???
 * \param dout ???
 * \param value ???
 */
Daniel Brown's avatar
Daniel Brown committed
306 307 308 309 310 311 312 313
void function_(int *din, double point[3], int *dout, double value[3]) {
    double k;
    double x, y, x1, y1;
    double delta_z;
    complex_t z1, z2, z3;
    double zx, zy;
    complex_t qx, qy;

314 315
    (void) *din;     // silence compiler warnings:
    (void) *dout;
Daniel Brown's avatar
Daniel Brown committed
316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342

    k = TWOPI / init.lambda * vlocal.lnr;

    x = point[0];
    y = point[1];

    qx.im = vlocal.lqx1.im;
    qy.im = vlocal.lqy1.im;
    zx = vlocal.lqx1.re;
    zy = vlocal.lqy1.re;

    x1 = x * vlocal.lcx + y * vlocal.lsy * vlocal.lsx;
    y1 = y * vlocal.lcy;
    // Vorzeichen experimetell bestimmt :(
    delta_z = -x * vlocal.lsx + y * vlocal.lsy * vlocal.lcx;

    qx.re = zx + delta_z;
    qy.re = zy + delta_z;

    z2 = z_by_phr(u_nm(vlocal.ln1, vlocal.lm1, qx, qy, x1, y1, vlocal.lnr),
            -1.0 * k * delta_z);
    z3 = cconj(u_nm(vlocal.ln2, vlocal.lm2, vlocal.lqx2, vlocal.lqy2, x, y,
            vlocal.lnr));
    z1 = z_by_z(z2, z3);

    value[0] = z1.re;
    value[1] = z1.im;
343 344
}

345 346