#include <stdio.h>
#include <stdlib.h>
#include <math.h>
// #include <omp.h>

#define phillu p2l*d*cos(alpha)
#define phrefl p2l*d*cos(beta)

#define pxlsize 7.56e-6
#define lambda 560e-9
#define p2l 2*M_PI/lambda

int main(int argc, char* argv[]) {
    int nbeta   = atoi(argv[1]);
    int nmirror = atoi(argv[2]);
    double alpha = atoi(argv[3]) * M_PI/180.;
    printf("#alpha = %.12E\n", alpha);
    printf("#b\tbeta\tAmplitude\tRealteil\tImaginärteil\n");
    for(int b = 0; b<nbeta; b++) {  //Winkel durchlaufen
        double R=0, I=0;            //Realteil und Imaginärteil
        double beta =  b * M_PI / nbeta;
        for (int i=-2; i<3; i+=2) { //Pixel duchlaufen
            for (float j=0.5; j<nmirror; j++) { //Elementarwellen eines Pixel durchlaufen
                double d = (i + j/nmirror)*pxlsize;
                double phase = phillu + phrefl;
                R+=cos(phase);
                I+=sin(phase);
            }
        }
        printf("%d\t%E\t%.12E\t%.12E\t%.12E\n",b,beta,sqrt(R*R + I*I),R,I);
    }
}

/* Short algorithm description:
 *
 *  1) if |x|==0:    sin(x)=x,
 *                   cos(x)=1.
 *  2) if |x|<2^-27: sin(x)=x-x*DP_SMALL, raising underflow only when needed,
 *                   cos(x)=1-|x|.
 *  3) if |x|<2^-5 : sin(x)=x+x*x^2*DP_SIN2_0+x^5*DP_SIN2_1,
 *                   cos(x)=1+1*x^2*DP_COS2_0+x^5*DP_COS2_1
 *  4) if |x|< Pi/4: sin(x)=x+x*x^2*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))),
 *                   cos(x)=1+1*x^2*(C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4)))).
 *  5) if |x| < 9*Pi/4:
 *      5.1) Range reduction:
 *          k=trunc(|x|/(Pi/4)), j=(k+1)&0x0e, n=k+1, t=|x|-j*Pi/4.
 *      5.2) Reconstruction:
 *          sign_sin = sign(x) * (-1.0)^(( n   >>2)&1)
 *          sign_cos =           (-1.0)^(((n+2)>>2)&1)
 *          poly_sin = ((((S4*t^2 + S3)*t^2 + S2)*t^2 + S1)*t^2 + S0)*t^2*t+t
 *          poly_cos = ((((C4*t^2 + C3)*t^2 + C2)*t^2 + C1)*t^2 + C0)*t^2*s+s
 *          if(n&2 != 0) {
 *              using cos(t) and sin(t) polynomials for |t|<Pi/4, results are
 *              cos(x) = poly_sin * sign_cos
 *              sin(x) = poly_cos * sign_sin
 *          } else {
 *              sin(x) = poly_sin * sign_sin
 *              cos(x) = poly_cos * sign_cos
 *          }
 *  6) if |x| < 2^23, large args:
 *      6.1) Range reduction:
 *          k=trunc(|x|/(Pi/4)), j=(k+1)&0xfffffffe, n=k+1, t=|x|-j*Pi/4
 *      6.2) Reconstruction same as (5.2).
 *  7) if |x| >= 2^23, very large args:
 *      7.1) Range reduction:
 *          k=trunc(|x|/(Pi/4)), j=(k+1)&0xfffffffe, n=k+1, t=|x|-j*Pi/4.
 *      7.2) Reconstruction same as (5.2).
 *  8) if x is Inf, return x-x, and set errno=EDOM.
 *  9) if x is NaN, return x-x.
 *
 * Special cases:
 *  sin/cos(+-0) = +-0/1 not raising inexact/underflow,
 *  sin/cos(subnormal) raises inexact/underflow,
 *  sin/cos(min_normalized) raises inexact/underflow,
 *  sin/cos(normalized) raises inexact,
 *  sin/cos(Inf) = NaN, raises invalid, sets errno to EDOM,
 *  sin/cos(NaN) = NaN.
 */

-    mx.expandMacros(exec);
+    mx.expandMacrosShellQuote(exec);