mirror

m.cpp

#include <iostream>
#include <cmath>
#include <vector>
#include <fstream>
#include <functional>

double fp(double phi, double a, double x, double y) {
    return std::tan(0.5 * (phi + std::atan(y / (x - a))));
}

double fm(double phi, double a, double x, double y) {
return std::tan(0.5 * (phi + std::atan(y / (x - a))));
}

double deg2rad(double deg) {
    const double PI = std::atan(1) * 4.;
    return deg / 180. * PI;
}

int main()
{
    const double PI = std::atan(1) * 4.;
    double lim_th = deg2rad(2.);
    double tag = deg2rad(3.);
    double source = -10.;
    double x, y;
    double k1, k2, k3, k4;
    const double dx = 0.001;
    double theta;
    std::function<double(double, double)> funcp = std::bind(fp, tag, source, std::placeholders::_1, std::placeholders::_2);
    std::function<double(double, double)> funcm = std::bind(fp, tag, source, std::placeholders::_1, std::placeholders::_2);

    std::ofstream file = std::ofstream("mirror.dat",std::ios::out);
     x = 0;
     y = 0;
     while (1) {
         k1 = funcm(x, y);
         k2 = funcm(x - 0.5 * dx, y - 0.5 * dx * k1);
         k3 = funcm(x - 0.5 * dx, y - 0.5 * dx * k2);
         k4 = funcm(x - dx, y - dx * k3);
         theta = (k1 + 2. * k2 + 2. * k3 + k4) / 6.;
         y -= dx * theta;
         x -= dx;
         //file << x << "  " << y << std::endl;

         if (std::atan(theta) < (tag-lim_th) * 0.5) { break; }
     }
     while (1) {
         k1 = funcp(x, y);
         k2 = funcp(x + 0.5 * dx, y + 0.5 * dx * k1);
         k3 = funcp(x + 0.5 * dx, y + 0.5 * dx * k2);
         k4 = funcp(x + dx, y + dx * k3);
         theta = (k1 + 2. * k2 + 2. * k3 + k4) / 6.;
         y += dx * theta;
         x += dx;
         file << x << "  " << y << std::endl;
         
         if (std::atan(theta) > (lim_th + tag) * 0.5) { break; }
     }
}