tut05_interface.cpp

/// @file tut05_interface.cpp
/// @brief Tutorial 05 — the L2 surface-physics layer (GSInterface and friends).
///
/// @page tutorial_05 Tutorial 05 — GSInterface: surface physics
///
/// This tutorial is a guided, runnable walkthrough of the L2 physics layer: the
/// abstract @ref GSInterface and its three concrete implementations
/// @ref GSInterfaceMirror, @ref GSInterfaceRefractive, and
/// @ref GSInterfaceAbsorber. An "interface" in GOST is not a C++ interface in the
/// Java sense — it is the physics of what happens to a ray when it strikes a
/// surface: does it reflect, bend, or stop? The layer is deliberately decoupled
/// from geometry (L1); an interface is handed only an incident direction and a
/// surface normal and returns the outgoing direction.
///
/// @section tut05_contract The one contract
/// Every method here obeys @c GetNextDirection(GSVector d, GSVector n): both @p d
/// (incident direction) and @p n (surface normal) are UNIT vectors in the WORLD
/// frame, and the returned outgoing direction is likewise a WORLD-frame unit
/// vector. There is no local-frame math at this layer — the frame bridge lives in
/// GSPlacement (see @ref tutorial_06).
///
/// @section tut05_topics What it demonstrates
/// - @b Mirror: specular reflection @f$ d' = d - 2(d\cdot n)n @f$, verified against
///   the analytic formula and the law of reflection (angle in == angle out).
/// - @b Refraction and the @b mu @b convention: @f$ \mu = n_{incident}/n_{transmitted} @f$
///   with @f$ \sin\theta_t = \mu\,\sin\theta_i @f$, checked numerically against Snell.
/// - The @b back-face @b auto-flip: when @f$ d\cdot n < 0 @f$ the routine sets
///   @f$ n\to -n,\ \mu\to 1/\mu @f$, making one stored surface reversible; shown by
///   firing the same object from both sides and proving reversibility.
/// - @b Total @b internal @b reflection past the critical angle (mu > 1), where
///   refraction falls back to specular reflection.
/// - The @b absorber and its out-of-band @c IsAbsorbing() termination signal.
///
/// @section tut05_build Build and run
/// Built as @c tut05_interface at the build root and linked against the GOST core
/// only (no ROOT). It runs standalone and doubles as a smoke test: every claim is
/// asserted and @c main() returns 0 on success, non-zero on any failure.
/// @code
/// g++ -std=c++11 -Iinclude tutorials/tut05_interface.cpp src/*.cpp -o tut05
/// @endcode

// tut05_interface.cpp — GOST tutorial 05: the L2 surface-physics layer.
//
// This is a guided, runnable walkthrough of GSInterface and its three concrete
// implementations. An "interface" in GOST is NOT a C++ interface in the Java
// sense — it is the physics of what happens to a ray when it strikes a surface:
// does it reflect (mirror), bend (refract), or stop (absorb)?
//
// The layer is deliberately decoupled from geometry (GSGeometry, L1). A
// GSInterface knows nothing about where the surface is; it is handed only an
// incident direction and the surface normal, and returns the outgoing direction.
// Geometry and interface are later glued together by GSElement (L3, see tut06).
//
// -------------------------------------------------------------------------
// THE ONE CONTRACT YOU MUST REMEMBER (DESIGN.md §4):
//   GSVector GetNextDirection(GSVector d, GSVector n) const;
//     * d — incident direction. UNIT vector, WORLD frame.
//     * n — surface normal at the hit point. UNIT vector, WORLD frame.
//     * returns — outgoing direction. UNIT vector, WORLD frame.
// Everything here operates on UNIT vectors expressed in the WORLD frame. There
// is no local-frame math at this layer; the frame bridge lives in GSPlacement.
// (GSVector's accessors are not const-correct, hence directions are passed BY
//  VALUE throughout the codebase — that is why the signature is not const-ref.)
// -------------------------------------------------------------------------
//
// This program also doubles as a smoke test: every claim below is checked with
// an Assert helper (like demo_procedure.cpp), and main() returns non-zero if any
// check fails. Build against the GOST core only (no ROOT):
//   g++ -std=c++11 -Iinclude tutorials/tut05_interface.cpp src/*.cpp -o tut05

#include "GSConstants.hpp"
#include "GSVector.hpp"
#include "GSInterface.hpp"
#include "GSInterfaceMirror.hpp"
#include "GSInterfaceRefractive.hpp"
#include "GSInterfaceAbsorber.hpp"

#include "cmath"
#include "iostream"
#include "iomanip"
#include "string"

// ------------------------------------------------------------------ utilities
static int gFail = 0;
static void Assert(bool cond, const std::string& what){
    std::cout << (cond ? "  [OK]   " : "  [FAIL] ") << what << std::endl;
    if(!cond) ++gFail;
}

// Component-wise vector equality (by value: GSVector accessors are non-const).
static bool VecEq(GSVector a, GSVector b, double tol = 1e-9){
    return std::fabs(a.X()-b.X()) < tol &&
           std::fabs(a.Y()-b.Y()) < tol &&
           std::fabs(a.Z()-b.Z()) < tol;
}

// sin of the angle between two vectors, via |a x b| / (|a||b|). With unit inputs
// the denominator is 1, so this is just |a x b| = sin(angle).
static double SinBetween(GSVector a, GSVector b){
    return a.Cross(b).Norm() / (a.Norm() * b.Norm());
}

int main(){
    std::cout << "======================================================\n";
    std::cout << " GOST tutorial 05 — the L2 surface-physics interfaces\n";
    std::cout << "======================================================\n";
    std::cout << std::fixed << std::setprecision(4);

    const double s2 = std::sqrt(2.0);
    const double deg = GSConstant::Pi() / 180.0;

    // ==== 1. GSInterfaceMirror — specular reflection ====================
    //
    // A perfect mirror reflects the incident ray about the surface normal. The
    // physics is the classic vector reflection formula:
    //
    //     d' = d - 2 (d . n) n
    //
    // This formula is symmetric in the sign of n (flipping n leaves d' unchanged),
    // so a mirror needs no notion of "front" vs "back" face. We verify (a) the
    // output equals the analytic formula, and (b) the law of reflection: the angle
    // of incidence equals the angle of reflection (measured from the normal).
    {
        std::cout << "\n---- 1. Mirror: specular reflection ----\n";
        GSInterfaceMirror mirror;

        GSVector n(0, 1, 0);                       // surface normal, +y
        GSVector d(1.0/s2, -1.0/s2, 0);            // 45 deg incidence, heading down
        GSVector out = mirror.GetNextDirection(d, n);

        // (a) Analytic reflection formula, computed independently.
        double k = d.Dot(n);
        GSVector analytic = d - (2.0 * k) * n;     // d - 2(d.n)n
        std::cout << "    in  " << d.Print(true) << "  n " << n.Print(true)
                  << "  ->  out " << out.Print(true) << "\n";
        Assert(VecEq(out, analytic), "mirror output equals d - 2(d.n)n");
        Assert(std::fabs(out.Norm() - 1.0) < 1e-9, "mirror output is a unit vector");

        // (b) Law of reflection: angle in == angle out (from the normal). The
        //     incident ray travels toward the surface, so its angle to n is taken
        //     against -d; the reflected ray leaves along +out.
        double angIn  = SinBetween((-1.0) * d, n); // sin of incidence angle
        double angOut = SinBetween(out, n);        // sin of reflection angle
        std::cout << "    sin(incidence)=" << angIn << "  sin(reflection)=" << angOut << "\n";
        Assert(std::fabs(angIn - angOut) < 1e-9, "angle of incidence == angle of reflection");

        // Normal incidence is the degenerate case: the ray bounces straight back.
        GSVector back = mirror.GetNextDirection(GSVector(0,-1,0), n);
        Assert(VecEq(back, GSVector(0,1,0)), "normal-incidence ray reflects straight back");

        Assert(!mirror.IsAbsorbing(), "a mirror is NOT absorbing (ray continues)");
    }

    // ==== 2. GSInterfaceRefractive — Snell's law & the mu convention =====
    //
    // A refracting (dielectric) surface bends the ray according to Snell's law.
    // GOST parameterises it by a SINGLE relative index mu, with the convention:
    //
    //     mu = n_incident_side / n_transmitted_side
    //
    // and, with the normal oriented along the propagation direction, the outgoing
    // ray satisfies:  sin(theta_t) = mu * sin(theta_i).
    //
    //   * mu < 1  => entering a DENSER medium (air->glass, mu = 1/1.5 ~ 0.667);
    //                the ray bends TOWARD the normal; TIR is impossible.
    //   * mu > 1  => entering a RARER medium (glass->air, mu = 1.5); the ray bends
    //                AWAY from the normal; TIR occurs past the critical angle.
    //
    // Here we send a ray into a front face (d.n > 0, no flip) with air->glass and
    // check the sine ratio numerically against Snell.
    {
        std::cout << "\n---- 2. Refraction: Snell's law (air -> glass) ----\n";
        const double nGlass = 1.5;
        const double mu = 1.0 / nGlass;            // air -> glass, mu ~ 0.667
        GSInterfaceRefractive glass(mu);
        Assert(std::fabs(glass.GetMu() - mu) < 1e-12, "GetMu() returns the stored mu");

        GSVector n(0, 1, 0);
        double thI = 30.0 * deg;                   // 30 deg incidence
        // Front-face incidence: d points generally ALONG n (d.n > 0), so no flip.
        GSVector d(std::sin(thI), std::cos(thI), 0);
        GSVector out = glass.GetNextDirection(d, n);

        double sinI = SinBetween(d, n);
        double sinT = SinBetween(out, n);
        std::cout << "    mu=" << mu << "  in " << d.Print(true)
                  << "  ->  out " << out.Print(true) << "\n";
        std::cout << "    sin(theta_i)=" << sinI << "  sin(theta_t)=" << sinT
                  << "  ratio sinT/sinI=" << sinT/sinI << "  (expect mu=" << mu << ")\n";
        Assert(std::fabs(out.Norm() - 1.0) < 1e-9, "refracted output is a unit vector");
        Assert(std::fabs(sinT - mu * sinI) < 1e-9, "sin(theta_t) = mu * sin(theta_i)");
        Assert(std::fabs(sinT - sinI / nGlass) < 1e-9, "matches Snell for n_glass = 1.5");
        // Denser medium => ray bends toward the normal => theta_t < theta_i.
        Assert(sinT < sinI, "ray bends TOWARD the normal entering a denser medium");
        Assert(!glass.IsAbsorbing(), "a refractive surface is NOT absorbing");
    }

    // ==== 3. The back-face auto-flip — one surface, both sides ==========
    //
    // A real surface is crossed from either side, but we store only ONE mu. The
    // routine makes the surface REVERSIBLE by auto-flipping when the ray arrives
    // on the back face: when d.n < 0 it sets  n -> -n  and  mu -> 1/mu, so the
    // stored mu always describes the front->back crossing. We demonstrate this by
    // firing at the SAME interface object from both sides.
    {
        std::cout << "\n---- 3. Back-face auto-flip (reversibility) ----\n";
        const double muFront = 1.0 / 1.5;          // front crossing: air -> glass
        GSInterfaceRefractive iface(muFront);
        GSVector n(0, 1, 0);
        double thI = 30.0 * deg;

        // Front-face hit: d.n > 0. Effective index = mu -> ray bends toward normal.
        GSVector dFront(std::sin(thI), std::cos(thI), 0);
        GSVector oFront = iface.GetNextDirection(dFront, n);
        double sinFrontT = SinBetween(oFront, n);
        Assert(std::fabs(sinFrontT - muFront * std::sin(thI)) < 1e-9,
               "front face: sin_t = mu * sin_i (bends toward normal)");

        // Back-face hit: d.n < 0. The routine flips to effective mu = 1/mu = 1.5,
        // so the SAME object now behaves as glass -> air and bends AWAY from normal.
        GSVector dBack(std::sin(thI), -std::cos(thI), 0);   // d.n = -cos(30) < 0
        GSVector oBack = iface.GetNextDirection(dBack, n);
        double sinBackT = SinBetween(oBack, n);
        std::cout << "    stored mu=" << muFront
                  << "   front sin_t=" << sinFrontT << " (bends in)"
                  << "   back sin_t=" << sinBackT << " (bends out)\n";
        Assert(std::fabs(sinBackT - (1.0/muFront) * std::sin(thI)) < 1e-9,
               "back face: effective index inverts to 1/mu (bends away from normal)");

        // The strongest proof of reversibility: refract front->back, send the
        // outgoing ray straight back through the surface, and we must recover the
        // exact reversed incident ray. Optical path is time-reversal symmetric.
        GSVector reversed = (-1.0) * oFront;               // fire oFront backwards
        GSVector recovered = iface.GetNextDirection(reversed, n);
        Assert(VecEq(recovered, (-1.0) * dFront),
               "refraction is reversible: reversing the exit ray recovers the entry ray");
    }

    // ==== 4. Total internal reflection (TIR) ============================
    //
    // Going from a denser to a rarer medium (mu > 1), beyond the CRITICAL angle
    // there is no real transmitted ray: sin(theta_t) = mu*sin(theta_i) would exceed
    // 1. Physics resolves this by reflecting the ray entirely. The GOST routine
    // detects this and falls back to the specular reflection d - 2(d.n)n.
    //
    // For glass->air (mu = 1.5) the critical angle has sin(theta_c) = 1/mu ~ 0.667
    // (theta_c ~ 41.8 deg). We pick 45 deg (sin = 0.707 > 0.667) to force TIR.
    {
        std::cout << "\n---- 4. Total internal reflection (glass -> air) ----\n";
        const double mu = 1.5;                     // glass -> air; TIR possible
        GSInterfaceRefractive glass(mu);
        double sinCrit = 1.0 / mu;
        std::cout << "    mu=" << mu << "  critical sin(theta_c)=" << sinCrit
                  << "  (theta_c~" << std::asin(sinCrit)/deg << " deg)\n";

        GSVector n(0, 1, 0);
        GSVector d(1.0/s2, 1.0/s2, 0);             // 45 deg, sin=0.707 > 0.667 -> TIR
        GSVector out = glass.GetNextDirection(d, n);

        // Past the critical angle the outgoing ray is the pure specular reflection.
        double k = d.Dot(n);
        GSVector reflected = d - (2.0 * k) * n;
        std::cout << "    in " << d.Print(true) << "  ->  out " << out.Print(true)
                  << "   (expected specular " << reflected.Print(true) << ")\n";
        Assert(std::fabs(out.Norm() - 1.0) < 1e-9, "TIR output is a unit vector");
        Assert(VecEq(out, reflected), "beyond critical angle refraction becomes specular reflection");

        // Cross-check: a standalone mirror gives the identical direction.
        GSInterfaceMirror mirror;
        Assert(VecEq(out, mirror.GetNextDirection(d, n)),
               "TIR direction equals a mirror's reflection");
    }

    // ==== 5. GSInterfaceAbsorber — the ray stops here ===================
    //
    // Detector faces and beam dumps ABSORB the ray: there is no outgoing ray at
    // all. A single "outgoing direction" cannot express "no ray", so absorption is
    // signalled out-of-band by IsAbsorbing() returning true. The propagator MUST
    // check IsAbsorbing() after a hit and, if true, terminate the ray instead of
    // trusting GetNextDirection() (which merely returns d unchanged as a formality).
    {
        std::cout << "\n---- 5. Absorber: ray termination ----\n";
        GSInterfaceAbsorber absorber;
        Assert(absorber.IsAbsorbing(), "absorber: IsAbsorbing() is true (propagator must stop the ray)");

        GSVector d(0.3, -0.7, 0.1);
        GSVector n(0, 1, 0);
        GSVector out = absorber.GetNextDirection(d, n);
        Assert(VecEq(out, d), "absorber: GetNextDirection returns d unchanged (a formality; do not propagate it)");
    }

    // ==== 6. Polymorphism through the base pointer ======================
    //
    // All three share the abstract GSInterface base, so the propagator can hold a
    // GSInterface* and let virtual dispatch pick the physics — including the
    // IsAbsorbing() decision. This is how GSElement stores "surface physics"
    // without knowing which kind it is (see tut06).
    {
        std::cout << "\n---- 6. Base-pointer dispatch ----\n";
        GSInterfaceMirror     mirror;
        GSInterfaceAbsorber   absorber;
        GSInterface* pm = &mirror;
        GSInterface* pa = &absorber;
        Assert(!pm->IsAbsorbing() && pa->IsAbsorbing(),
               "IsAbsorbing() dispatches correctly through GSInterface*");
    }

    std::cout << "\n------------------------------------------------------\n";
    if(gFail == 0){
        std::cout << "tut05_interface: all checks passed.\n";
        return 0;
    }
    std::cout << "tut05_interface: " << gFail << " check(s) FAILED.\n";
    return 1;
}