| test formattazione [messaggio #84379] |
mer, 22 giugno 2011 01:06 |
 Tommaso Russo, Triest
Messaggi: 14
Registrato: luglio 2011 |
Junior Member |
|
|
Ciao to everybody,
I've encountered frequently, in forums and news groups, questions about
a rectangle inscribed in another rectangle, but nowhere a general
discussion on ways to solve related problems (i.e., given three of the
five quantities involved - 2 sides of outer rectangle, 2 sides of inner
rectangle, and one of the two complementary angles between them -
compute remaining two, or state that the problem has no solution).
Probably the problem is felt too trivial to encourage somebody to afford
difficulties - really NOT trivial - that are encountered in some of
these cases.
I did it, and attach below my solutions and some Fortran 95 subroutines
and programs that allow for finding unknowns in each case - or for
stating that the problem, with input data given, has no solution.
Hope this will help somebody in future.
Of course, observations, corrections, and suggestions for improvements,
are welcome.
Code has been written by myself and I release it in the public domain.
june 2011
Tommaso Russo, Trieste(Italy)
Rectangle inscribed in another rectangle: solutions for all cases.
============================================================ ======
Let X and Y be the sides of a rectangle XY exactly inscribed in another
rectangle AB of sides A and B, i.e. each vertex of XY belongs to a side
of AB. Let us call phi (0°<=phi<=90°) one of the angles formed by a side
Y of XY and a side A of AB.
It is graphically obvious that:
(1) A = Y cos(phi) + X sin(phi) ( = a1 + a2)
(2) B = Y sin(phi) + X cos(phi) ( = b1 + b2)
Where we have put
a1 = Y cos(phi)
a2 = X sin(phi) = A - a1
b1 = Y sin(phi)
b2 = X cos(phi) = B - b1
with the obvious meaning of distances between a vertex of AB and one of XY.
We will treat three non trivial cases:
1). Given A,B,phi, find X and Y,
or state that the problem has no solution.
2). Given A,X,Y, find B (and phi),
or state that the problem has no solution.
2.a). Trigonometrical solution
2.b). Analytical solution
3). Given A,B,Y, find X (and phi),
or state that the problem has no solution.
3.a) Graphical solution
3.b) Analytical solution
3.c) Numerical solution
As a byproduct of solutions, useful for a graphical construction of
inner rectangle, also a1,a2,b1 and b2 will be found.
Remaining possible cases can be solved trivially: i.e.,
0.1). Given X,Y,phi, find A and B:
solution is given directly by (1) and (2).
0.2). Given A,Y,phi, find X and B:
a1 = Y cos(phi)
a2 = A-a1 (if negative, there is no solution)
X = a2/sin(phi)
B is given by (2).
Case 1). Given A,B,phi, find X and Y
====================================
This problem is encountered in practice when a picture has been taken,
holding the camera non perfectly right, and one wants to rotate the
photo by an angle phi, so that the horizon becomes horizontal; and then,
cut the photo's oblique sides, saving as much as possible of the image.
This is the simplest case: it is enough to solve the system of equations
(1) and (2) in unknowns X and Y.
X = [B cos(phi)-A sin(phi)]/[cos(phi)^2-sin(phi)^2],
Y = [A cos(phi)-B sin(phi)]/[cos(phi)^2-sin(phi)^2].
If both X and Y turn out non-negative, this is the required solution.
And, of course,
a1 = Y cos(phi)
a2 = X sin(phi)
b1 = Y sin(phi)
b2 = X cos(phi).
In the particular case in which the denominator [cos(phi)^2-sin(phi)^2]
is zero, i.e. phi=45°, solutions exist only if A=B, and in this case X
can assume any value between 0 and A*sqrt(2); Y is determined accordingly:
Y = A*sqrt(2) - X.
Case 2). Given A,X,Y, find B (and phi)
======================================
If BOTH X and Y are greater than A, the problem has no solution.
If the diagonal of XY ( sqrt(X^2 + Y^2) ) is less than A, the problem
has no solution.
OTHERWISE:
2.a). Trigonometrical solution
------------------------------
If Y is less than or equal to A, the problem has at least one solution.
Proof:
Let us draw two parallel straight lines, r and s, at distance A.
Let us insert rectangle XY between r and s, with sides X parallel to
them. This can be done because Y<=A.
Let us rotate XY, until both endpoints of one of its diagonals touch r
and s. This can be done because its diagonal is greater than or equal to A.
Let us draw two other straight lines, orthogonal to both r and s,
passing throught the endpoints of the other diagonal of XY. Their
distance is a value for B that satisfies (1) and (2).
Now:
Let us call theta the acute (or rect) angle between the first diagonal
and r (or s):
theta = arcsin( A / sqrt(X^2 + Y^2) )
Let us call alpha the angle between the same diagonal and one of the X
sides:
alpha = arctan(Y/X)
It turns out
phi = theta - alpha
= arcsin(A/sqrt(X^2 + Y^2)) - arctan(Y/X)
Once found phi, B can be computed from (2):
B = Y sin(phi) + X cos(phi)
and obviously
a1 = Y cos(phi)
a2 = X sin(phi)
b1 = Y sin(phi)
b2 = X cos(phi).
BESIDES:
Also if X is less than or equal to A, the problem has at least one solution.
Proof proceeds as before:
Let us draw two parallel straihgt lines, r and s, at distance A.
Let us insert rectangle XY between r and s, this time with sides Y
parallel to them. This can be done because X<=A.
Let us rotate XY, until both endpoints of one of its diagonals touch r
and s. This can be done because its diagonal is greater than or equal to A.
Let us draw two other straight lines, orthogonal to both r and s,
passing throught the endpoints of the other diagonal of XY. Their
distance is a value for B that satisfies (1) and (2).
In this case, alpha will be the angle between one diagonal of XY and one
of its Y (not X) sides, i.e. the complementary of that seen above; and hence
alpha = arctan(X/Y)
phi = theta - alpha
= arcsin(A/sqrt(X^2 + Y^2)) - arctan(X/Y).
Once found another value for phi, another value for B can be computed
from (2):
B = Y sin(phi) + X cos(phi)
and obviously, four other values for a1,a2,b1 and b2:
a1 = Y cos(phi)
a2 = X sin(phi)
b1 = Y sin(phi)
b2 = X cos(phi).
Thus, if BOTH X and Y are <= A, there are TWO solutions for B (that
coincide if X=Y, i.e. XY is a square).
Of course, on can choose between them the smallest one.
You can find, in the appendix below, a subroutine than implements the
described algorithm and a test program that calls it. Both are written
in Fortran 95 and can be compiled and run on a linux PC using f95.
2.b). Analytical solution
-------------------------
putting cos(phi) = c in (1)
(1) A = Y cos(phi) + X sin(phi)
we obtain
(1') A = Y c + X sqrt(1-c^2)
A - Y c = X sqrt(1 - c^2)
(A - Yc)^2 = X^2 (1 - c^2)
(Y^2+X^2) c^2 - 2AY c +(A^2-C^2) = 0
This is a II degree equation in c that has two solutions:
c_1 = {A Y + sqrt[A^2 Y^2 - (A^2-X^2) (X^2+Y^2)]} / (X^2+Y^2)
c_2 = {A Y - sqrt[A^2 Y^2 - (A^2-X^2) (X^2+Y^2)]} / (X^2+Y^2)
Acceptable solutions must be real (the term under square root must be
non-negative) and within 0 and 1 included (phi=arcos(c) must be acute).
Under these conditions, c1(X) e c2(X) coincide with cos(phi) found in
two subcases of case 2.a) seen above.
From both values of c, a value of phi can be found as arcos(c), and
then B, and then a1,a2,b1,b2, as above.
Caso 3). Given A,B,Y, find X (and phi)
======================================
This is the hardest case. It can be solved graphically, analitycally and
numerically.
3.a) Graphical solution
3.b) Analytical solution
3.c) Numerical solution
3.a). Graphical solution
-----------------------
If Y is greater than the diagonal of AB, Y > sqrt(A^2+B^2), the problem
has no solution.
Otherwise, solutions can be one or more.
To find them, let us consider the solution found in case 2): once found
phi, we obtained
B = Y sin(phi) + X cos(phi)
a1 = Y cos(phi)
b1 = Y sin(phi)
we can also write, straightfully,
X = a2/sin(phi) = (A-a1)/sin(phi)
= (A-Y*cos(phi))/sin(phi)
and so
B = Y sin(phi) + [(A-Y*cos(phi))/sin(phi)] cos(phi).
Let us now consider phi as the independent variable, and plot curves
B(phi), a1(phi), b1(phi), X(phi), and the horizontal straight line
corresponding to the input value of B.
Using gnuplot, we could write (or copy and paste):
A = 1 # replace with your input value
Y = .8 # replace with your input value
Binput = 1 # replace with your input value
phi(x)=x # phi = x, i.e. x-axis *is* angle phi
X(x) = (A-Y*cos(phi(x)))/sin(phi(x))
B(x) = Y*sin(phi(x)) + X(x)*cos(phi(x))
a1(x)= Y*cos(phi(x))
b1(x) = Y*sin(phi(x))
set angles degrees
set xrange [0:90] # phi must be acute
set yrange [0:2*Binput]
set samples 10000 # enough for a good resolution,
# low enough to make zooms in less than 1 second
plot B(x), Binput, X(x), a1(x), b1(x)
Let's now locate intersections of B(x) with the horizontal straight line
y=Binput: between 0° and 90°, there can be 0, 1, 2 or 3 of them (this is
shown below, in 3.b)).
For each intersection, zooming on it, let's locate corresponding angle
(i.e. its x coordinate).
With input data shown in the example given, A = 1, Binput = 1, Y = .8,
we find
phi = 17.1144°, 45°, 72.8856°.
Zooming on vertical straight lines passing throught each intersection,
we can locate corresponding values for X(phi): positive ones can be
accepted; negative ones must be dropped, together with corresponding
value of phi.
In the example given, they turn out to be
X = 0.8, 0.61421, 0.8
that, from a geometrical point of view, is quite obvious: all squares
with sides Y < 1 can be inscribed in a square 1 by 1, but also a
rectangle with two sides of given lenght Y, forming two isosceles
triangles at 45° with two opposite verteces of the square 1 by 1 (second
solution). The third solution is symmetric of the first one, with
respect to both axes of the 1 by 1 square.
Having obtained three values for phi, it is easy to find corresponding
values of a1 and b1, for an easy drawing of three inner rectangles.
3.b). Analytical solution
------------------------
Any manipulation of (1) and (2), aimed at eliminating all unkowns but
one, will lead to a quartic (fourth degree) equation in the remaining
unkown.
This appears somehow discouraging, because, while most people knows how
to solve a quadratic equation with the well-known formula x=[-b +-
sqrt(b^2-4ac)]/(2a), not everybody suspect that similar, thought more
complex, formulas exist to compute exactly complex radices of a cubic
and of a quartic polynomials. They have been found, respectively, by
Girolamo Cardano and by his pupil Lodovico Ferrari, and published
together in "Ars Magna" in 1545.
Here they are:-)
http://planetmath.org/encyclopedia/QuarticFormula.html
Ferrari's algorithm looks actractive, because it is NOT iterative:
however, it is NOT suggested for solutions of quartic equations: it has
been said "beautiful mathematics but a bad algorithm", expensive,
unstable and inaccurate. Use of Jenkins-Traub_algorithm is suggested
instead ( http://en.wikipedia.org/wiki/Jenkins-Traub_algorithm ); or,
create a companion matrix of the polynomial (
http://en.wikipedia.org/wiki/Companion_matrix ) and compute its
eigenvalues with QR method ( http://en.wikipedia.org/wiki/QR_algorithm ).
However, if implemented on a common PC, in this case Ferrari's algorithm
works, following some precautions:
- use complex double precision variables;
- accept as valid solutions "real" numbers with imaginary part
nonzero, negative lenghts, cosines greater than 1, if these differences
vanish in conversion to single precision variables.
Geometrically, rectangular triangles, cut by inner rectangle out of
outer one, are all similar: hence
(3) b1 / a1 = (A-a1) / (B-b1)
and, for Pythagoras' theorem,
(4) a1^2 + b1^2 = Y^2
(5) (A-a1)^2 + (B-b1)^2 = X^2
i.e.
(3') b1(B-b1) = a1(A-a1)
Bb1-b1^2 = Aa1-a1^2
Bb1 = Aa1 + b1^2 - a1^2 elevando al quadrato
(3") B^2b1^2 = A^2a1^2 + b1^4 + a1^4 + 2Aa1b1^2 - 2Aa1^3 -2 b1^2a1^2
(4') b1^2 = Y^2 - a1^2
replacing b1^2 from (4') into (3")
B^2(Y^2-a1^2) = A^2a1^2 + (Y^2-a1^2)^2 + a1^4 + 2Aa1(Y^2-a1^2) -
2Aa1^3 - 2(Y^2-a1^2)a1^2
B^2 Y^2 - B^2 a1^2 = A^2 a1^2 - 4 A a1^3 + 2 A a1 Y^2 + 4 a1^4 -
4 a1^2 Y^2 + Y^4
and, after rearrangement,
4a1^4 - 4Aa1^3 + (A^2+B^2-4Y^2)a1^2 + 2AY^2a1 + (Y^4 - B^2Y^2) = 0
that is a quartic in unkown a1 in canonical form:
ax^4 + bx^3 + cx^2 + dx + e = 0
where
a = 4
b = - 4A
c = A^2+B^2-4Y^2
d = 2AY^2
e = Y^4 - B^2Y^2
An algorithm to compute step by step its solutions using Ferrari's
method is well exposed here:
http://en.wikipedia.org/wiki/Quartic_equation#Solving_a_quar tic_equation
"subroutine Lodovico_Ferrari", that can be found in appendix, follows
strictly this method.
It solves correctly quartic equations built starting from their
solutions, e.g. (x-1)(x-7)(x+2/3)(x+50)=0 or (x-3)(x-15)(X^2+1)=0. It is
written in fortan 95, together with a subroutine that computes a,b,c,d,e
from input data, and a test program. Them all can be compiled with f95
and run on a Linux machine.
3.c). Numerical solution
------------------------
Follows same guidelines of given graphical solution:
X(phi) = (A-Y*cos(phi))/sin(phi)
B(phi) = Y*sin(phi) + X(phi)*cos(phi)
and the problem is reduced to finding zeroes of B(phi)-Binput in phi
interval ]0°,90°]
(0 is excluded because, for phi=0°, B(phi) turns out +infinite,
-infinite, or indeterminate.)
Zeroes of B(phi)-Binput can be found iteratively with several
algorithms, the simplest one being bisection. However, them all require
knowledge of contiguous intervals [phi_1,phi_2], where B(phi)-Binput is
defined and continuous and changes sign only once; i.e. phi_1 and phi_2
"bracket" a zero, and only one zero.
Such intervals must be identified before.
Let us study behaviour of function B(phi) in interval [0°,90°],
considering A and Y parameters.
It's (geometrically) obvious that B(90°) = Y.
Limit of B(phi) for phi->0 has a geometrical meaning only if Y<=A;
otherwise, phi is inferiorly limited by the need of inserting into
rectangle AxB at least a rectangle Yx0 (i.e. a segment of lenght Y): it is
phi_min = arcos(A/Y), e
B(phi_min) = Y sin(phi_min).
case 3.c.1): Y >= A
In this case, if Y>A, B(phi) has always a maximum between phi_min and
90°. If Y=A, B(0°) turns out indetermined: in facts, for any Binput,
choosing X=Binput suffices to obtain a solution (outer rectangle
inscribes itself). However, for phi->0, B(phi) tends to 0, and, in
semi-open interval ]0°,90°], B(phi) and also X(phi) are always positive.
If Y<A instead, lim_phi->0(B(phi))=+infinite, and, for what concerns
behaviour of B(phi), one must distinguish two subcases:
case 3.c.2): Y =< A/sqrt(2)
In this case, between 0° and 90° B(phi) is monotone and decreasing.
In the particular case Y = A/sqrt(2), B(phi) has an horizontal
inflection point for phi=45°, where B(phi)=A; anyway, between 0° e 90°
B(phi) is decreasing.
Hence, if Y =< A/sqrt(2), for any Binput>=Y one and only one solution
exists; for Binput<Y, there are none.
case 3.c.3): A > Y > A/sqrt(2)
In this case, B(phi) has a minimum for phi between 0° and 45°, and a
maximum for phi between 45° and 90°. For any Binput>=Y at least one
solution exists (but there can be 3 of them); for Binput<Y there are 2
solutions, or none.
Hence, the problem can be solved, in each of these cases, determining
maxima ad minima of B(phi) in significant intervals seen, with simple
trisection method (or golden section method); ]0°,90°] is, so, divided
in subintervals in which B(phi) is monotone. In these of them, in which
B(phi)-Binput has at endpoins values (or limits) of opposite signs, the
only point where B(phi)=Binput can be found with simple bisection method
(or equivalent and faster methods, e.g. secants or Newton–Raphson).
Find below, in the appendix, a routine that implements the described
algorithm, and a test program that uses it, written in Fortran 95. Both
can be compiled with f95 and run on a Linux machine.
------------------------------------------------------------
Previous assertions sub 3.c.1), 3.c.2) e 3.c.3) can be proved as follows:
B(phi) can be written in the form
B(phi) = Y*sin(phi) + (A-Y*cos(phi))/sin(phi)*cos(phi)
Without loss of generality, we can put A=1. This is equivalent to use of
A as unit of measure for lenghts.
We obtain
B(x) = y*sin(x) + ((1-y*cos(x))/sin(x))*cos(x)
where we have put x = phi*pi/180, i.e. phi measured in radians, and y =
Y/A, i.e. Y measured in units A.
Deriving B with respect to x,
dB/dx = 2*y*cos(x)+cosec(x)**2*(y*cos(x)-1)
= y*cos(x)*[1/sin(x)**2 + 2] - 1/sin(x)**2
= cosec(x)**2 * [y*cos(x)*(1+2*sin(x)**2) - 1]
in the interval ]0,pi/2] for x, corresponding to ]0°,90°] for phi,
cosec(x)>0; so minima and maxima can fall only where
(a) y*cos(x)*(1+2*sin(x)**2) = 1.
The value of cos(x)*(1+2*sin(x)**2) is 1 for x=0, 0 for x=pi/2
(phi=90°), and its maximum value is sqrt(2) where x=pi/4 (phi=45°).
So, for y = sqrt(1/2), equation (a) can be verified straightforwardly
for x=pi/4 (phi=45°); but B(x) has there an horizontal inflection point.
For any other value of x in the interval, the left member of (a) is <1.
A fortiori, if y < sqrt(1/2), the left member of (a) is always <1, and
so B(x) is monotone.
For y > sqrt(1/2), but < 1, the value of left member of (a) is y (i.e. <
1) for x=0, 0 for x=pi/2 (phi=90°), ad for x=pi/4 (phi=45°) it has a
maximum y*sqrt(2) (that is > 1). So, it equals 1 for two values of x,
one in the interval 0<x<pi/4 (or 0°<phi<45°) and on in pi/4<x<pi/2
(45°<phi<90°).
For y=1, the left member of (a) is 1 for x=0; the second stationary
point of B(x) (that turns out to be a maximum) falls in the interval
pi/4<x<pi/2 (45°<phi<90°).
For y>1, the value of left member of (a) is > 1 for x=0, and (a) holds
only in one point, in the interval pi/4<x<pi/2 (45°<phi<90°). And
finally, to show that this maximum of B falls between phi_min and 90°,
even if phi_min>45°, it's enough to show geometrically that, for
phi=phi_min, B(phi) is increasing: in facts, having inserted exactly a
rectangle Yx0 in AxB as its diagonal, any rotation of one of its Y sides
requires B to increase.
------------------------------------------------------------
Appendix
========
Subroutine and test program for case 2.a). (Fortran 95):
--------------------------------------------------------
! copy from here
!23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
subroutine find_B_sides_by_angles(A,X,Y, containable,B_min_cont,
& inscribableX,B_ins_X,a1X,a2X,b1X,b2X,
& inscribableY,B_ins_Y,a1Y,a2Y,b1Y,b2Y)
! A,X,Y input, others are output
! A is a side of outer rectangle AxB; X,Y are sides of candidate inner
! containable means XxY containable in AxInfinity
! with sides PARALLEL to those of AxB;
! inscribable means XxY in inscribable in AxSomeBSide
! B_min_cont will contain minimum B for which XxY is containable in AxB
! with sides PARALLEL to those of AxB
! (meaningless if containable is .false.).
!
! if X <= A, inscribableX is .true., and
! B_ins_X will contain a B for which XxY is inscribable
! and a1X,a2X,b1X,b2X
! are distances from verteces of A x B and verteces of X x Y
! along sides A and B ( a1X+a2X=A, b1X+b2X=B_ins_X )
! (meaningless if inscribableX is .false.).
!
! if Y <= A, inscribableY is .true.
! B_ins_Y will contain a B for which YxX is inscribable
! and a1Y,a2Y,b1Y,b2Y
! are distances from verteces of A x B and verteces of X x Y
! along sides A and B ( a1X+a2X=A, b1X+b2X=B_ins_Y )
! (meaningless if inscribableY is .false.).
!
! if X<=A and Y <=A,
! both inscribableX and inscribableY are .true.
! but in general B_ins_X =/= B_ins_Y.
logical containable, inscribableX, inscribableY
if(A.lt.0.0.or.X.lt.0.0.or.Y.lt.0.0) then
print *, "find_B_sides - negative input argument"
call abort
endif
if (min(X,Y).gt.A) then
containable=.false.
inscribableX=.false.
inscribableY=.false.
return
endif
! here min(X,Y).le.A, hence X x Y is surely
! containable in A x SomeBSide
containable=.true.
if(X.le.A.and.Y.gt.A) B_min_cont = Y
if(Y.le.A.and.X.gt.A) B_min_cont = X
if(Y.le.A.and.X.le.A) B_min_cont = min(X,Y)
! but, if its diagonal < A, X x Y is NOT inscribable in A x AnySide
if (X**2+Y**2.lt.A**2) then
inscribableX=.false.
inscribableY=.false.
return
endif
if(Y.le.A) then
! if(Y.le.A) XxY IS SURELY inscribable in A x SomeSide.
!
! Proof:
! - insert X x Y between two parallel lines at distance A
! with X side also parallel to them
! (can be done, because Y<=A)
! - rotate X x Y by an angle phi until it touchs both parallel lines
! (can be done, because its diagonal >= A)
! - draw two lines, orthogonal to previous parallel lines,
! passing for free verteces.
! Their distance is B_ins_Y.
!
! Now:
! (1) Y cos(phi) + X sin(phi) = A
! (2) Y sin(phi) + X cos(phi) = B_ins_Y
!
! But:
!
! angle theta formed by diagonal and the parallel lines is
!
! theta = asin (A/diag(XxY))
!
! angle alpha formed by diagonal and XxY axis parallel to Y
!
! alpha = atan((X/2)/(Y/2)) = atan(X/Y)
!
! angle phi is theta - alpha
theta = asin (A/sqrt(X**2+Y**2))
alpha = atan(Y/X)
phi = theta - alpha
B_ins_Y = Y * sin(phi) + X * cos(phi)
inscribableY=.true.
a1Y = Y*cos(phi)
a2Y = X*sin(phi)
b1Y = Y*sin(phi)
b2Y = X*cos(phi)
else
inscribableY=.false.
endif
! Now, if also side X of XxY is <= A, the procedure can
! be repeated giving a different result:
!
if(X.le.A) then
! - insert X x Y between two parallel lines at distance A
! with Y SIDE, NOT X SIDE, also parallel to them
! (can be done, because X<A)
! - rotate X x Y by an angle phi until it touchs both parallel lines
! (can be done, because its diagonal >= A)
! - draw two lines, orthogonal to previous parallel lines,
! passing for free verteces.
! Their distance is ANOTHER B_ins, let call it B_ins_X.
! Again, /mutata mutandis/:
!
! (1) X cos(phi) + Y sin(phi) = A
! (2) X sin(phi) + Y cos(phi) = B_ins_X
!
! as before,angle theta formed by diagonal and the parallel lines is
!
! theta = asin (A/diag(XxY))
!
! but angle alpha formed by diagonal and XxY axis parallel to X is
!
! alpha = atan((Y/2)/(X/2)) = atan(Y/X)
!
! as before, angle phi is theta - alpha.
theta = asin (A/sqrt(X**2+Y**2))
alpha = atan(X/Y)
phi = theta - alpha
B_ins_X = X * sin(phi) + Y * cos(phi)
inscribableX=.true.
a1X = X*sin(phi)
a2X = Y*cos(phi)
b1X = X*cos(phi)
b2X = Y*sin(phi)
else
inscribableX=.false.
endif
return
end
program checkcontainability
logical containable, inscribableX, inscribableY, EOF /.false./
do while (.not.EOF)
print *, "--------------------------------------------------"
print *, "enter A,B (outer rectangle), X,Y (inner rectangle)"
READ (*, *, IOSTAT=IST) A,B,X,Y
EOF = IST.LT.0
IF (.NOT.EOF) THEN
print *, "outer rectangle:",A," x ",B,
& ", inner rectangle: ",X," x ",Y
call find_B_sides_by_angles(A,X,Y, containable,B_min_cont,
& inscribableX,B_ins_X,a1X,a2X,b1X,b2X,
& inscribableY,B_ins_Y,a1Y,a2Y,b1Y,b2Y)
print *,"containable: ",containable," B_min_cont:",B_min_cont
print *,"inscribableX:",inscribableX, " B_ins_X:",B_ins_X
print *,"inscribableY:",inscribableY," B_ins_Y:",B_ins_Y
if(inscribableY) then
print*,"inner rectangle is inscribable in",
& " a rectangle ",A,"x ",B_ins_Y
print *, "distances between verteces of ",
& "inner and of outer rectangles, a1,a2,b1,b2:"
print *, a1Y,a2Y,b1Y,b2Y
endif
if(inscribableX) then
print*,"inner rectangle is inscribable in",
& " a rectangle ",A,"x ",B_ins_X
print *, "distances between verteces of ",
& "inner and of outer rectangles, a1,a2,b1,b2:"
print *, a1X,a2X,b1X,b2X
endif
if (containable) then
print *, "XxY is containable in", A," x ",B_min_cont,
& " with sides PARALLEL to those of AxB"
endif
if (inscribableX .and.
& abs(B-B_ins_X)/B.lt.10.**(-precision(B))) then
print *, "XxY is exactly inscribable in AxB"
endif
if (inscribableY .and.
& abs(B-B_ins_Y)/B.lt.10.**(-precision(B))) then
print *, "XxY is exactly inscribable in AxB"
endif
ENDIF
end do
end
! copy up to here
END Subroutine and test program for case 2.a).
----------------------------------------------
Subroutine and test program for case 3.b). (Fortran 95):
--------------------------------------------------------
! copy from here
!23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
subroutine Lodovico_Ferrari(a,b,c,d,e, x)
! finds analitically solutions of quartic ax^4+bx^2+cx^3+dx+e=0
! from http://en.wikipedia.org/wiki/Quartic_equation
implicit complex(8) (A-H,O-Z)
real(8) a,b,c,d,e
complex(8) x(1:4) ! x is a vector to host four complex solutions
alpha = -(3.*b**2)/(8.*a**2) + c/a
beta = b**3./(8.*a**3)-(b*c)/(2.*a**2)+d/a
gamma = -(3.*b**4)/(256.*a**4) + (c*b**2)/(16.*a**3) -
& (b*d)/(4.*a**2) + e/a
print *, "beta = ",beta
if(beta.eq.0.) then
! sign(t)=+
x(1)=-b/(4*a)+ sqrt( .5*(-alpha + sqrt(alpha**2-4*gamma)) )
x(2)=-b/(4*a)- sqrt( .5*(-alpha + sqrt(alpha**2-4*gamma)) )
! sign(t)=-
x(3)=-b/(4*a)+ sqrt( .5*(-alpha - sqrt(alpha**2-4*gamma)) )
x(4)=-b/(4*a)- sqrt( .5*(-alpha - sqrt(alpha**2-4*gamma)) )
return
else
P = -(alpha**2)/12. - gamma
Q = -(alpha**3)/108.+alpha*gamma/3.-(beta**2)/8.
R = -Q/2.+sqrt((Q**2)/4.+(P**3)/27.)
U = R**(1./3.)
if(U.ne.0.) y= -5.*alpha/6. +U -P/(3.*U)
if(U.eq.0.) y= -5.*alpha/6. -Q**(1./3.)
W = sqrt(alpha+2.*y)
! sign(s)=+
x(1)=-b/(4*a)+.5*(+W+sqrt(-(3*alpha+2.*y+2*beta/W)))
x(2)=-b/(4*a)+.5*(+W-sqrt(-(3*alpha+2.*y+2*beta/W)))
! sign(s)=-
x(3)=-b/(4*a)+.5*(-W+sqrt(-(3*alpha+2.*y-2*beta/W)))
x(4)=-b/(4*a)+.5*(-W-sqrt(-(3*alpha+2.*y-2*beta/W)))
endif
return
end
subroutine find_max_X_side_analytically(sideA,sideB,sideY,
& sideXout,phiout,a1out,a2out,b1out,b2out)
real(8) a,b,c,d,e
complex(8) x(1:4) ! x is a vector for four complex solutions
logical xIsReal
a = 4.
b = - 4.*sideA
c = sideA**2+sideB**2-4.*sideY**2
d = 2.*sideA*sideY**2
e = sideY**4 - sideB**2*sideY**2
print*, "a,b,c,d,e", a,b,c,d,e
call Lodovico_Ferrari(a,b,c,d,e, x(1),x(2),x(3),x(4))
print*,"Ferrari returned: ", (x(i),i=1,4)
! x(1)..x(4) are possible values for a1, i.e. Y cos(phi)
sideXout=-1 ! "if unchanged, no solution found"
do i=1,4
print *, "abs(img(x(",i,")))/sideA:",abs(imagpart(x(i))/sideA)
xIsReal=.true.
if(abs(imagpart(x(i))/sideA).gt.10.**(-precision(sideA)))
& xIsReal=.false.
if(xisReal)then
a1=REALPART(x(i))
! rounding problems can lead to negative a1, avoid them
if(abs(a1)/sideA.lt.10.**(-precision(sideA))) a1=0.
! rounding problems can lead to try arcos(1.+eps), avoid them
if(abs(sideA-a1)/sideA.lt.10.**(-precision(sideA))) a1=sideA
a2=sideA-a1
phi=acos(a1/sideY)
b1=sideY*sin(phi)
b2=sideB-b1
sideX = sqrt(a2**2+b2**2)
print "(a,i1,a,a)", "from x(",i,"): phi in radians,",
& "phi in degrees, sideX, a1,a2,b1,b2"
print*, phi, phi*45./atan(1.), sideX, a1,a2,b1,b2
if(sideX.ge.0.0.and.a1.ge.0.0.and.
& a1.le.sideA.and.sideX.gt.sideXout) then
sideXout=sideX
phiout=phi
a1out=a1
a2out=a2
b1out=b1
b2out=b2
endif
endif
end do
end
program test_find_max_X_side_analytically
logical EOF /.false./
print *, "machine precision(real) is ",precision(sideA)
do while (.not.EOF)
print *, "--------------------------------------------------"
print *, "enter sideA,sideB",
& "(outer rectangle), sideY (inner rectangle)"
READ (*, *, IOSTAT=IST) sideA,sideB,sideY
EOF = IST.LT.0
IF (.NOT.EOF) THEN
print *, "outer rectangle:",sideA," x ",sideB,
& ", inner rectangle: ? x ",sideY
call find_max_X_side_analytically(sideA,sideB,sideY,
& sideX,phi,a1,a2,b1,b2)
if(sideX.lt.0.) then
print *, " No solution"
else
print *, "*** Best solution for second inner side ***"
print *, "max value for sideX =",sideX
print *, "angle between A and Y (radians) =", phi
print *, "angle between A and Y (degrees) =", phi*45/atan(1.)
print *, "distances between verteces of ",
& "inner and of outer rectangles:"
print *, a1,a2,b1,b2
print *, " verify:"
print *, "Y cos(phi) + X sin(phi), i.e.",sideY*cos(phi)+
& sideX*sin(phi), "should be equal to A: ", sideA
print *, "X cos(phi) + Y sin(phi), i.e.",sideX*cos(phi)+
& sideY*sin(phi), "should be equal to B: ", sideB
print *,"atan2(a1,b1)+atan2(a2,b2), i.e.",atan2(a1,b1)
& +atan2(a2,b2), "should be equal to 90°, i.e.", atan(1.)*2.
endif
ENDIF
end do
end
! copy up to here
END Subroutine and test program for case 3.b).
----------------------------------------------
Subroutine and test program for case 3.c). (Fortran 95):
--------------------------------------------------------
! copy from here
!23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
function B(A,Y,phi)
! phi in radians, in range 0° - 90°
if (phi.eq.0.0) then
if(Y.lt.A) then
B = huge(B) ! is nearly +infinite
else if (Y.eq.A) then
B=0.0 ! makes B continuous in 0
else
B= -huge(B) ! is nearly -infinite
endif
else
X = (A-Y*cos(phi))/sin(phi)
B = Y*sin(phi) + X*cos(phi)
endif
return
end
function phi_of_Bmin (A,Y, phiinit,phiend)
! phi1<phi2 and they bracket a local minimum of B
! simple trisection method: golden section would be
! twice faster, but in this case who cares...
phi1 = phiinit
phi2 = phiend
do while ((phi2-phi1)/phi2.gt.10.**(-precision(phi)))
phiu = phi1 + (phi2-phi1)/3.
phiv = phi2 - (phi2-phi1)/3.
Bu = B(A,Y,phiu)
Bv = B(A,Y,phiv)
if (Bu.lt.Bv) then ! searching for a min!
phi2 = phiv
else if (Bu.gt.Bv) then ! searching for a min!
phi1 = phiu
else if (Bu.eq.Bv) then
phi1 = phiu
phi2 = phiv
endif
end do
phi_of_Bmin = (phi1+phi2)/2.
return
end
function phi_of_Bmax (A,Y, phiinit,phiend)
! phi1<phi2 and they bracket a local maximum of B
! simple trisection method: golden section would be
! twice faster, but in this case who cares...
phi1 = phiinit
phi2 = phiend
do while ((phi2-phi1)/phi2.gt.10.**(-precision(phi)))
phiu = phi1 + (phi2-phi1)/3.
phiv = phi2 - (phi2-phi1)/3.
Bu = B(A,Y,phiu)
Bv = B(A,Y,phiv)
if (Bu.gt.Bv) then ! searching for a max!
phi2 = phiv
else if (Bu.lt.Bv) then ! searching for a max!
phi1 = phiu
else if (Bu.eq.Bv) then
phi1 = phiu
phi2 = phiv
endif
end do
phi_of_Bmax = (phi1+phi2)/2.
return
end
function phi_of_B_equal_Binput (A,Y,Binput,phiinit,phiend)
! phiinit<phiend and they must bracket a zero of (B-Binput)
! simple bisection method: secant or Newton–Raphson method
! would be faster, but in this case who cares...
phi1 = phiinit
phi2 = phiend
B1 = B(A,Y,phi1)
if(B1.eq.Binput)then
phi_of_B_equal_Binput = phiinit
return
endif
B2 = B(A,Y,phi2)
if(B2.eq.Binput)then
phi_of_B_equal_Binput = phiend
return
endif
do while ((phi2-phi1)/phi2.gt.10.**(-precision(phi)))
phi = phi1 + (phi2-phi1)/2.
Bmiddle = B(A,Y,phi)
if (Bmiddle.eq.Binput) then
phi_of_B_equal_Binput = phi
return
else if ((B1-Binput)*(Bmiddle-Binput).gt.0.) then
phi1 = phi
B1 = Bmiddle
else
phi2 = phi
B2 = Bmiddle
endif
end do
phi_of_B_equal_Binput = phi
return
end
subroutine find_angles_numerically(A,Binput,Y, phiout,ns)
real phiout(1:4)
! phi is a vector for at most 3 positive solutions:
! but sometimes there will be apparently four
! (problem is malconditioned for A=B and Y=~A/sqrt(2),
! and for A=B, A/sqrt(2)<Y<A, solution at 45° will be found twice.)
ns = 0 ! number of solutions found (here 0)
pi = 4.*atan(1.) ! greek pi
! case 3.c.1
if(Y.ge.A) then
! max between phi_start and 90°,
! 0 or 1 or 2 solutions
phi_start = acos(A/Y)
phi_end = pi/2.
phiBmax = phi_of_Bmax (A,Y, phi_start, phi_end)
print *, "phiBmax=", phiBmax, phiBmax*180./pi
if ((B(A,Y,phi_start).le.Binput).and.
& (B(A,Y,phiBmax).ge.Binput)) then
ns = ns+1
phiout(ns) =
& phi_of_B_equal_Binput (A,Y,Binput, phi_start,phiBmax)
endif
if ((B(A,Y,phi_end).le.Binput).and.
& (B(A,Y,phiBmax).ge.Binput)) then
ns = ns+1
phiout(ns) =
& phi_of_B_equal_Binput (A,Y,Binput, phiBmax,phi_end)
endif
! case 3.c.2
else if (Y.le.A/sqrt(2.)) then
! one solution only between 0 and 90°
ns = ns+1
phiout(ns) = phi_of_B_equal_Binput (A,Y,Binput, 0.,pi/2.)
! case 3.c.3
else if (Y.gt.A/sqrt(2.).and.Y.lt.A) then ! redundant but clearer
! 1 min between 0° and 45°,
! 0 or 1 or 2 solutions
phi_start = 0.
phi_end = pi/4.
phiBmin = phi_of_Bmin (A,Y, phi_start, phi_end)
print *, "phiBmin=", phiBmin, phiBmin*180./pi
if ((B(A,Y,phi_start).ge.Binput).and.
& (B(A,Y,phiBmin).le.Binput)) then
ns = ns+1
phiout(ns) =
& phi_of_B_equal_Binput (A,Y,Binput, phi_start,phiBmin)
endif
if ((B(A,Y,phi_end).ge.Binput).and.
& (B(A,Y,phiBmin).le.Binput)) then
ns = ns+1
phiout(ns) =
& phi_of_B_equal_Binput (A,Y,Binput, phiBmin,phi_end)
endif
! 1 max between 45° and 90°,
! 0 or 1 or 2 solutions
phi_start = pi/4.
phi_end = pi/2.
phiBmax = phi_of_Bmax (A,Y, phi_start, phi_end)
print *, "phiBmax=", phiBmax, phiBmax*180./pi
if ((B(A,Y,phi_start).le.Binput).and.
& (B(A,Y,phiBmax).ge.Binput)) then
ns = ns+1
phiout(ns) =
& phi_of_B_equal_Binput (A,Y,Binput, phi_start,phiBmax)
endif
if ((B(A,Y,phi_end).le.Binput).and.
& (B(A,Y,phiBmax).ge.Binput)) then
ns = ns+1
phiout(ns) =
& phi_of_B_equal_Binput (A,Y,Binput, phiBmax,phi_end)
endif
endif
end
program test_find_X_S_numerically
logical EOF /.false./
real phi(1:4)
! phi is a vector for at most 3 positive solutions
! negative if not valid
print *, "machine precision(real) is ",precision(A)
do while (.not.EOF)
print *, "--------------------------------------------------"
print *, "enter A,B",
& "(outer rectangle), Y (inner rectangle)"
READ (*, *, IOSTAT=IST) A,Binput,Y
EOF = IST.LT.0
IF (.NOT.EOF) THEN
print *, "outer rectangle:",A," x ",Binput,
& ", inner rectangle: ? x ",Y
call find_angles_numerically(A,Binput,Y, phi,ns)
if(ns.eq.0.) then
print *, " No solution"
else
if(ns.eq.4.and.phi(2).eq.phi(3)) then
! may occurr if one solution is exactly 45°
! an is found twice, in two contiguous intervals
! (occurrs if A = B)
phi(3) = phi(4)
ns = 3
endif
do i=1,ns
!! read *,dummy
X = (A-Y*cos(phi(i)))/sin(phi(i))
Bobtained = Y*sin(phi(i)) + X*cos(phi(i))
a1 = Y*cos(phi(i))
b1 = Y*sin(phi(i))
a2 = A-a1
b2 = Bobtained - b1
print *, "*** solution for inner side number",i
print *, "value for X =",X
print *, "B obtained = ", Bobtained
print *, "angle between A and Y (radians) =", phi(i),
& " (degrees) =", phi(i)*45/atan(1.)
print *, "distances between verteces of ",
& "inner and of outer rect:"
print *, a1,a2,b1,b2
print *, " verify:"
print *, "Y cos(phi) + X sin(phi), i.e.",Y*cos(phi(i))+
& X*sin(phi(i)), "should be equal to A: ", A
print *, "X cos(phi) + Y sin(phi), i.e.",X*cos(phi(i))+
& Y*sin(phi(i)), "should be equal to B: ", Binput
print *,"atan2(a1,b1)+atan2(a2,b2), i.e.",atan2(a1,b1)
& +atan2(a2,b2), "should be equal to 90°, i.e.", atan(1.)*2.
end do
endif
ENDIF
end do
end
! copy up to here
END Subroutine and test program for case 3.c).
----------------------------------------------
--
TRu-TS
Buon vento e cieli sereni
|
|
|
My forum, my way! Il forum dei newsgroup
Utenti
F.A.Q.
Registrati
Login
Home
E-mail ad un amico 
