#include #include #include #include #include "mucal.h" #include "panel.h" /* ********************************************** */ vekt normalize( vekt cp, int prt) /* noormaliizes cp */ { double lin; vekt co; lin = sqrt(cp.x * cp.x + cp.y * cp.y + cp.z * cp.z); if (prt) printf("Inside normalize: x: %11.8lf, y: %11.8lf, z: %11.8lf, len: %8.2lf\n", cp.x, cp.y, cp.z, lin); co.x = cp.x / lin; co.y = cp.y / lin; co.z = cp.z / lin; co.len = sqrt(co.x * co.x + co.y * co.y + co.z * co.z); if (prt) printf("At exit from normalize: x: %11.8lf, y: %11.8lf, z: %11.8lf, len: %8.2lf\n\n", co.x, co.y, co.z, co.len); return(co); } vekt crossp( vekt a, vekt b) /*calculates cross product of a and b */ { vekt cp; cp.x = a.y * b.z - b.y * a.z; cp.y = -a.x * b.z + a.z * b.x; cp.z = a.x * b.y - a.y * b.x; cp.len = sqrt(cp.x * cp.x + cp.y * cp.y + cp.z * cp.z); cp.x /= cp.len; cp.y /= cp.len; cp.z /= cp.len; /* printf("Mirror crossp cx %lf ay %lf bz %lf by %lf az %lf\n", cp.x, a.y, b.z, -b.y, a.z); printf("Mirror crossp cy %lf ax %lf bz %lf az %lf bz %lf\n", cp.y, -a.x, b.z, a.y, b.z); printf("Mirror crossp cz %lf ay %lf bz %lf by %lf az %lf\n", cp.z, a.x, b.y, -a.y, b.x); printf("Mirror crossp len %lf cpx %lf cpy %lf by %lf az %lf\n", cp.len, cp.x, cp.y, cp.z); */ return(cp); } /* QQQQ1 Begin of classical (optical) mirror routin */ vekt mirror( vekt Ri, vekt N, int pt) /* calculates mirror vector, Ro, of Ri mirrored in normplane to N. (the mirror-plane) * We split the input vector, Ri, into two composants, Ai and Bi, where Ai is parallel to the mirror plane * and Bi is perpendicular to this plane. The outgoing ray, after reflection, then has the direction: * [vector_Ro] = [vector_Ai] - [vector_Bi]; * */ { vekt C, D, Ai, Bi, Ro; double lenB, lenBi; int prt = 1; if (prt) printf("Useq Mirror \n"); if (prt) printf("Useq Mirror input Ri xyz: %+lf %+lf %+lf, len: %lf\n", Ri.x, Ri.y, Ri.z, Ri.len); if (prt) printf("Useq Mirror input N xyz: %+lf %+lf %+lf, len: %lf\n", N.x, N.y, N.z, N.len); C = crossp(Ri, N); /* [C] is not a unity vector [C] lies in the mirror plane but is perpendicu,ar to both [Ri] and [N] The length is the sine of the angle between [Ri] and [N], that is, to the cosine of the angle between [Ri] and the mirror plane, that is, the length of the projection of [Ri] on the mirror plane */ if (prt) printf("Useq Mirror 1 C xyz: %+lf %+lf %+lf, len: %lf\n", C.x, C.y, C.z, C.len); Ai = crossp(C, N); /* vector [Ai] is the trace of vector [Ri] in the mirror plane */ /* if (prt) printf("Useq Mirror D xyz: %+lf %+lf %+lf, len: %lf\n", D.x, D.y, D.z, D.len); Ai.x = -D.x * C.len; N is a unit vector, Ai is not! Ai.y = -D.y * C.len; Ai.z = -D.z * C.len; the [Ai]-vector is the composant of Ri in the mirror plane Ai.len = sqrt(Ai.x * Ai.x + Ai.y * Ai.y + Ai.z * Ai.z); */ if (prt) printf("Useq Mirror 2 Ai xyz: %+lf %+lf %+lf, len: %lf %lf (Ai is the Ri composant along mirror plane)\n", Ai.x, Ai.y, Ai.z, Ai.len, C.len); Bi.x = Ri.x - N.x * C.len; /* Bi is not a unit vector. */ Bi.y = Ri.y - N.y * C.len; Bi.z = Ri.z - N.z * C.len; lenBi = sqrt(Bi.x * Bi.x + Bi.y * Bi.y + Bi.z * Bi.z); lenB = dotp(Ri, N); /* [Ri] and [N] are unit vectors, lenB is the csine to the angle between them */ if (prt) printf("Useq Mirror Ro_tmp Bi xyz: %+lf %+lf %+lf, lenB=sin(Ri/N): %lf, lenBi: %lf, Ai.len**2+Bi.len**2: %lf\n", Bi.x, Bi.y, Bi.z, lenB, lenBi, sqrt(Ai.len*Ai.len + Bi.len*Bi.len)); Bi.x *= lenB / lenBi; /* Bi is not a unit vector */ Bi.y *= lenB / lenBi; /* Bi is composant of Ri parallel to normplane */ Bi.z *= lenB / lenBi; Bi.len = sqrt(Bi.x * Bi.x + Bi.y * Bi.y + Bi.z * Bi.z); if (prt) printf("Useq Mirror tmp Bi xyz: %+lf %+lf %+lf, len: %lf (Bi is the Ri composant perpendicular to mirror plane.)\n", Bi.x, Bi.y, Bi.z, Bi.len); Ro.x = Ai.x - Bi.x; /* Ro is not a unit vector here */ Ro.y = Ai.y - Bi.y; /* is the direction of the outgoing ray after diffraction */ Ro.z = Ai.z - Bi.z; Ro.len = sqrt(Ro.x * Ro.x + Ro.y * Ro.y + Ro.z * Ro.z); if (prt) printf("Useq Mirror tmp Ro xyz: %+lf %+lf %+lf, len: %lf (Ro = Bi - Ai)\n", Ro.x, Ro.y, Ro.z, Ro.len); Ro.x /= Ro.len; /* Normalize vector components of Ro */ Ro.y /= Ro.len; Ro.z /= Ro.len; Ro.len = sqrt(Ro.x * Ro.x + Ro.y * Ro.y + Ro.z * Ro.z); if (prt) printf("Useq Mirror output Ro xyz: %+lf %+lf %+lf, len: %lf (Ro normalized) \n", Ro.x, Ro.y, Ro.z, Ro.len); return(Ro); } /* QQQQ1 End of classical (optical) mirror routin */ double dotp( vekt a, vekt b) /*calculates dot product of a and b */ { double dot; int prt = 0; if (prt) printf("Useq Dotp tmp A xyz: %+lf %+lf %+lf, len: %lf\n", a.x, a.y, a.z, a.len); if (prt) printf("Useq Dotp tmp B xyz: %+lf %+lf %+lf, len: %lf\n", b.x, b.y, b.z, b.len); dot = a.x * b.x + a.y * b.y + a.z * b.z; if (prt) printf("Useq Dotp output dot. %+lf\n", dot); return(dot); } double angl( vekt a, vekt b) /* calculates angle (radians) between a and b */ { /* assumes that the vectors a and b are normalized */ vekt cp; double angle; cp.x = a.y * b.z - b.y * a.z; cp.y = -a.x * b.z + a.z * b.x; cp.z = a.x * b.y - b.y * a.x; cp.len = sqrt(cp.x * cp.x + cp.y * cp.y + cp.z * cp.z); angle = asin(cp.len); return(angle); } /* QQQQ 2 */ vekt mirBRG( vekt Ri, vekt N, double thetaB, int pt) /* calculates mirror vector, Ro, of Ri mirrored in normplane to N. (the mirror-plane) * We split the input vector, Ri, into two composants, Ai and Bi, where Ai is parallel to the mirror plane * and Bi is perpendicular to this plane. The outgoing ray, after reflection, then has the direction: * [vector_Ro] = [vector_Ai] - [vector_Bi]; * */ { vekt C, D, Ai, Bi, Ro; double lenB, lenBi, lenBrg; int prt = 1; if (prt) printf("Useq Mirror \n"); if (prt) printf("Useq Mirror input Ri xyz: %+lf %+lf %+lf, len: %lf\n", Ri.x, Ri.y, Ri.z, Ri.len); if (prt) printf("Useq Mirror input N xyz: %+lf %+lf %+lf, len: %lf\n", N.x, N.y, N.z, N.len); C = crossp(Ri, N); /* [C] is not a unity vector [C] lies in the mirror plane but is perpendicu,ar to both [Ri] and [N] The length is the sine of the angle between [Ri] and [N], that is, to the cosine of the angle between [Ri] and the mirror plane, that is, the length of the projection of [Ri] on the mirror plane */ if (prt) printf("Useq Mirror 1 C xyz: %+lf %+lf %+lf, len: %lf\n", C.x, C.y, C.z, C.len); Ai = crossp(C, N); /* vector [Ai] is the trace of vector [Ri] in the mirror plane */ /* if (prt) printf("Useq Mirror D xyz: %+lf %+lf %+lf, len: %lf\n", D.x, D.y, D.z, D.len); Ai.x = -D.x * C.len; N is a unit vector, Ai is not! Ai.y = -D.y * C.len; Ai.z = -D.z * C.len; the [Ai]-vector is the composant of Ri in the mirror plane Ai.len = sqrt(Ai.x * Ai.x + Ai.y * Ai.y + Ai.z * Ai.z); */ if (prt) printf("Useq Mirror 2 Ai xyz: %+lf %+lf %+lf, len: %lf %lf (Ai is the Ri composant along mirror plane)\n", Ai.x, Ai.y, Ai.z, Ai.len, C.len); lenBrg = sin(2.0 * thetaB); /* force correct Bragg-angle corresponnding to current energy */ Ro.x = Ri.x + N.x * lenBrg; Ro.y = Ri.y + N.y * lenBrg; Ro.z = Ri.z + N.z * lenBrg; Ro = normalize(Ro, pt); if (prt) printf("Useq Mirror output Ro xyz: %+lf %+lf %+lf, len: %lf (Ro normalized) lenBrg: %+lf \n", Ro.x, Ro.y, Ro.z, Ro.len, lenBrg); return(Ro); lenB = dotp(Ri, N); /* Ri is a unit vector */ Bi.x = Ri.x - N.x * C.len; /* Bi is not a unit vector. */ Bi.y = Ri.y - N.y * C.len; Bi.z = Ri.z - N.z * C.len; lenBi = sqrt(Bi.x * Bi.x + Bi.y * Bi.y + Bi.z * Bi.z); if (prt) printf("Useq Mirror tmp Bi xyz: %+lf %+lf %+lf, sin(Ri/N): %lf, lenBi: %lf, Ai.len**2+Bi.len**2: %lf\n", Bi.x, Bi.y, Bi.z, lenB, lenBi, sqrt(Ai.len*Ai.len + Bi.len*Bi.len)); Bi.x *= lenB / lenBi; /* Bi is not a unit vector */ Bi.y *= lenB / lenBi; /* Bi is composant of Ri parallel to normplane */ Bi.z *= lenB / lenBi; Bi.len = sqrt(Bi.x * Bi.x + Bi.y * Bi.y + Bi.z * Bi.z); if (prt) printf("Useq Mirror tmp Bi xyz: %+lf %+lf %+lf, len: %lf (Bi is the Ri composant perpendicular to mirror plane.)\n", Bi.x, Bi.y, Bi.z, Bi.len); Ro.x = Ai.x - Bi.x; /* Ro is not a unit vector here */ Ro.y = Ai.y - Bi.y; /* is the direction of the outgoing ray after diffraction */ Ro.z = Ai.z - Bi.z; Ro.len = sqrt(Ro.x * Ro.x + Ro.y * Ro.y + Ro.z * Ro.z); if (prt) printf("Useq Mirror tmp Ro xyz: %+lf %+lf %+lf, len: %lf (Ro = Bi - Ai)\n", Ro.x, Ro.y, Ro.z, Ro.len); Ro.x /= Ro.len; /* Normalize vector components of Ro */ Ro.y /= Ro.len; Ro.z /= Ro.len; Ro.len = sqrt(Ro.x * Ro.x + Ro.y * Ro.y + Ro.z * Ro.z); if (prt) printf("Useq Mirror output Ro xyz: %+lf %+lf %+lf, len: %lf (Ro normalized) \n", Ro.x, Ro.y, Ro.z, Ro.len); return(Ro); } /* QQQQ3 */