Sicily 1719. Ellipse Intersection

1719. Ellipse Intersection

Constraints

Time Limit: 1 secs, Memory Limit: 32 MB

Description

Recently the astronomers have discovered a peculiar pair of planets, named A and B. As we know, a planet usually moves in an ellipse orbit, so do A and B. But their orbits are quite special:
(1) Their orbits are in the same plane and share the same center;
(2) The segments formed by joining their own focuses respectively are perpendicular to each other.
If we denote the center as point O, and the focuses of orbit A as F1, F2, then we can build Cartesian coordinates, with point O being the origin, and the line through F1,F2 being the x-axis.
Here is a sample as follows:

As the astronomers would like to know more about the planets, they decide to calculate the intersection area first. But unfortunately calculating the intersection area is so complex and beyond their ability, so they have to turn to you, a talent of programming. Now your task is: given the description of two ellipse orbits satisfying the conditions above, calculate their intersection area.

Input

Input may contain multiple test cases. The first line is a positive integer n (n<=100), denoting the number of test cases below. Each case is composed of two lines, the first one is the description of orbit A, and the second one the description of orbit B. each description contains two positive integers, a, b, (a, b<=100) denoting the ellipse equation x²/a²+y²/b²=1. It’s guaranteed that the focuses of orbit A are on x-axis, and the focuses of orbit B are on y-axis.

Output

For each test case, output one line containing a single real number, the area of the intersection, accurate to three decimals after the decimal point.

Sample Input

1
2 1
1 2

Sample Output

3.709

// Problem#: 1719
// Submission#: 3585206
// The source code is licensed under Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License
// URI: http://creativecommons.org/licenses/by-nc-sa/3.0/
// All Copyright reserved by Informatic Lab of Sun Yat-sen University
#include <stdio.h>
#include <math.h>

const double pi = acos(-1.0);
double a1, b1, a2, b2;

double solve() {
    if (a1 <= a2) return pi * a1 * b1;
    if (b1 >= b2) return pi * a2 * b2;
    double ans;
    double x, y;
    x = sqrt((b1 * b1 - b2 * b2) * a1 * a1 * a2 * a2 / (b1 * b1 * a2 * a2 - b2 * b2 * a1 * a1));
    y = sqrt((a1 * a1 - a2 * a2) * b1 * b1 * b2 * b2 / (a1 * a1 * b2 * b2 - b1 * b1 * a2 * a2));
    double theta1 = asin(y / b1), theta2 = asin(y / b2);
    ans = a1 * b1 * ((pi - 2 * theta1) + sin(2 * theta1)) + a2 * b2 * (2 * theta2 + sin(2 * theta2)) - 4 * x * y;
    return ans;
}

int main() {
    int n;
    scanf("%d", &n);
    while (n--) {
        scanf("%lf%lf%lf%lf", &a1, &b1, &a2, &b2);
        printf("%.3lf\n", solve());
    }
    return 0;
}