fork download
copy 
  1. /* 
  2.  * File:   main.cpp
  3.  * Author: Captain Giraffe
  4.  *
  5.  
  6. In laser physics, a "white cell" is a mirror system that acts as a delay line for the laser beam. The beam enters the cell, bounces around on the mirrors, and eventually works its way back out.
  7.  
  8. The specific white cell we will be considering is an ellipse with the equation 4x2 + y2 = 100
  9.  
  10. The section corresponding to −0.01 ≤ x ≤ +0.01 at the top is missing, allowing the light to enter and exit through the hole.
  11.  
  12. The light beam in this problem starts at the point (0.0,10.1) just outside the white cell, and the beam first impacts the mirror at (1.4,-9.6).
  13.  
  14. Each time the laser beam hits the surface of the ellipse, it follows the usual law of reflection "angle of incidence equals angle of reflection." That is, both the incident and reflected beams make the same angle with the normal line at the point of incidence.
  15.  
  16. In the figure on the left, the red line shows the first two points of contact between the laser beam and the wall of the white cell; the blue line shows the line tangent to the ellipse at the point of incidence of the first bounce.
  17.  
  18. The slope m of the tangent line at any point (x,y) of the given ellipse is: m = −4x/y
  19.  
  20. The normal line is perpendicular to this tangent line at the point of incidence.
  21.  
  22. The animation on the right shows the first 10 reflections of the beam.
  23.  
  24. How many times does the beam hit the internal surface of the white cell before exiting?
  25.  
  26.  * Created on July 8, 2012, 3:19 PM
  27.  */
  28.  
  29. #include <cstdlib>
  30. #include <utility>
  31. #include <cmath>
  32. #include <iostream>
  33. #include <cassert>
  34.  
  35. using std::cout;
  36. using std::endl;
  37.  
  38. struct point : std::pair<double, double> {
  39.  
  40. 	point() : x(first), y(second) {
  41. 	}
  42.  
  43. 	point(double x, double y) : std::pair<double, double>(x, y), x(first), y(second) {
  44. 	}
  45.  
  46. 	point(const point & p) : std::pair<double, double>(p.first, p.second), x(first), y(second) {
  47. 	}
  48.  
  49. 	point& operator=(const point & p) {
  50. 		this->first = p.first;
  51. 		this->second = p.second;
  52. 		return *this;
  53. 	}
  54. 	double &x;
  55. 	double &y;
  56. };
  57.  
  58. point operator +(point p1, point p2) {
  59. 	return point(p1.x + p2.x, p1.y + p2.y);
  60. }
  61.  
  62. point operator -(point p1, point p2) {
  63. 	p2.x = -p2.x;
  64. 	p2.y = -p2.y;
  65. 	return p1 + p2;
  66. }
  67.  
  68. double angle(point p1, point p2) {
  69. 	point p = p2 - p1;
  70. 	return atan2(p.y, p.x);
  71. }
  72.  
  73. /**
  74.  * Solves
  75.  * a*x*x + b*x + c = 0 
  76.  * @param a
  77.  * @param b
  78.  * @param c
  79.  * @return 
  80.  */
  81. point quadraticSolve(double a, double b, double c) {
  82. 	double discriminant = b * b - 4 * a*c;
  83. 	assert(discriminant >= 0); // is should be two good roots always.
  84.  
  85. 	point solution;
  86. 	solution.first = (-b + sqrt(discriminant)) / (2 * a);
  87. 	solution.second = (-b - sqrt(discriminant)) / (2 * a);
  88. 	return solution;
  89. }
  90.  
  91. /**
  92.  * Calculates the slope for the line containing these points
  93.  * @param p1 from 
  94.  * @param p2 to
  95.  * @return the slope k
  96.  */
  97. double calcK(point p1, point p2) {
  98. 	return (p2.y - p1.y) / (p2.x - p1.x);
  99. }
  100.  
  101. double calcM(point p, double slope) {
  102. 	double m = -slope * p.x + p.y;
  103. 	return m;
  104.  
  105. }
  106.  
  107. const double EPSILON = 1E-12;
  108.  
  109. /*
  110.  * Direct model
  111.  */
  112. int main(int argc, char** argv) {
  113.  
  114. 	point prev(0.0, 10.1);
  115. 	point next(1.4, -9.6);
  116. 	double k, m;
  117.  
  118. 	// Need to calculate step 2 outside the loop since, 1.0, 10.1 is not on the ellipsis
  119. 	cout << (k = calcK(prev, next)) << endl;
  120. 	cout << (m = calcM(prev, k)) << endl;
  121. 	double a = (4 + k * k);
  122. 	double b = (2 * k * m);
  123. 	double c = m * m - 100;
  124. 	cout << "a: " << a << ", b: " << b << ", c:" << c << endl;
  125. 	point solved = quadraticSolve(a, b, c);
  126.  
  127. 	// get the root not relating to this point
  128. 	double nextX = next.x;
  129. 	double nextY = k * nextX + m;
  130. 	//prev = next;
  131. 	next = point(nextX, nextY);
  132. 	double slopeLine = calcK(prev, next);
  133. 	double slopeEllipsis = -4 * next.x / next.y;
  134. 	point diff = next - prev;
  135. 	double incidentAngle = atan2(diff.y, diff.x);
  136. 	double reversedIncident = (incidentAngle > 0) ? incidentAngle - M_PI : incidentAngle + M_PI;
  137. 	double ellipseAngle = atan(slopeEllipsis);
  138. 	double ellipseNormal = M_PI / 2 + ellipseAngle;
  139. 	double angleDiff = reversedIncident - ellipseNormal;
  140. 	double angleReflected = reversedIncident - 2 * angleDiff;
  141. 	double slopeReflection = tan(angleReflected);
  142. 	k = slopeReflection;
  143. 	m = calcM(next, slopeReflection);
  144. 	a = (4 + k * k);
  145. 	b = (2 * k * m);
  146. 	c = m * m - 100;
  147. 	solved = quadraticSolve(a, b, c);
  148.  
  149. 	// get the root not relating to this point
  150. 	nextX = (abs(solved.first - next.x) < EPSILON) ? solved.second : solved.first;
  151. 	nextY = k * nextX + m;
  152. 	prev = next;
  153. 	next = point(nextX, nextY);
  154. 	bool foundExit = false;
  155. 	while ( foundExit == false) {
  156.  
  157. 		slopeLine = calcK(prev, next);
  158. 		slopeEllipsis = -4 * next.x / next.y;
  159. 		point diff = next - prev;
  160. 		incidentAngle = atan2(diff.y, diff.x);
  161. 		reversedIncident = (incidentAngle > 0) ? incidentAngle - M_PI : incidentAngle + M_PI;
  162. 		ellipseAngle = atan(slopeEllipsis);
  163. 		ellipseNormal = M_PI / 2 + ellipseAngle;
  164. 		angleDiff = reversedIncident - ellipseNormal;
  165. 		angleReflected = reversedIncident - 2 * angleDiff;
  166. 		slopeReflection = tan(angleReflected);
  167. 		k = slopeReflection;
  168. 		m = calcM(next, slopeReflection);
  169. 		a = (4 + k * k);
  170. 		b = (2 * k * m);
  171. 		c = m * m - 100;
  172. 		solved = quadraticSolve(a, b, c);
  173.  
  174. 		// get the root not relating to this point
  175. 		nextX = (abs(solved.first - next.x) < EPSILON) ? solved.second : solved.first;
  176. 		nextY = k * nextX + m;
  177. 		prev = next;
  178. 		next = point(nextX, nextY);
  179. 		cout << next.x << endl;
  180. 		cout << next.y << endl;
  181. 		if ( abs(next.x) < 0.01){
  182. 			if ( next.y > 0){
  183. 				foundExit = true;
  184. 			}
  185. 		}
  186. 	}
  187. 	return 0;
  188. }
  189.  
  190.  
Success #stdin #stdout 0.01s 2684KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
-14.0714
10.1
a: 202.005, b: -284.243, c:2.01
0.324583
9.97891
https://ideone.com/L7KMu
language:
C++ (gcc 8.3)
created:
11 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!