/*

kissing ellipses



These animations are constructed by shrinking and rotating a sequence of concentric and similar ellipses,
so that each ellipse lies inside the previous ellipse and is tangent to it.

https://benice-equation.blogspot.com/2019/01/nested-ellipses.html

==================================================
https://math.stackexchange.com/questions/3773593/given-ellipse-of-axes-a-and-b-find-axes-of-tangential-and-concentric-ellips

tangential concentric ellipse and insribed ellipses

Let’s say I have an ellipse with horizontal axis $a$ and vertical axis $b$, centered at $(0,0)$. 
I want to compute $a’$ and $b’$ of a smaller ellipse centered at $(0,0)$, 
with the axes rotated by some angle $t$, tangent to the bigger ellipse and $\frac{a’}{b’}=\frac{a}{b}$.




---------------------

The standard parametric equation is:

(x,y)->(a cos(t),b sin(t))


---------------------------

Rotation counterclockwise about the origin through an angle α carries 

(x, y) to (x cos α − ysin α, ycos α+x sin α) 

https://www.maa.org/external_archive/joma/Volume8/Kalman/General.html

=====================================
https://math.stackexchange.com/questions/2987044/how-to-find-the-equation-of-a-rotated-ellipse

===============================
https://math.stackexchange.com/questions/3773593/given-ellipse-of-axes-a-and-b-find-axes-of-tangential-and-concentric-ellips

============================================================
intersection of 2 ellipses

the common point of 2 ellipses are not vertices ( vertex)

https://math.stackexchange.com/questions/1688449/intersection-of-two-ellipses
https://math.stackexchange.com/questions/425366/finding-intersection-of-an-ellipse-with-another-ellipse-when-both-are-rotated/425412#425412

https://math.stackexchange.com/questions/3312747/intersection-area-of-concentric-ellipses
https://math.stackexchange.com/questions/426150/what-is-the-general-equation-of-the-ellipse-that-is-not-in-the-origin-and-rotate/434482#434482
------
xc <- 1 # center x_c or h
yc <- 2 # y_c or k
a <- 5 # major axis length
b <- 2 # minor axis length
phi <- pi/3 # angle of major axis with x axis phi or tau

t <- seq(0, 2*pi, 0.01) 
x <- xc + a*cos(t)*cos(phi) - b*sin(t)*sin(phi)
y <- yc + a*cos(t)*cos(phi) + b*sin(t)*cos(phi)
plot(x,y,pch=19, col='blue')
https://stackoverflow.com/questions/41820683/how-to-plot-ellipse-given-a-general-equation-in-r

===============
Batch file for Maxima CAS
save as a e.mac
run maxima : 
 maxima
and then : 
batch("e.mac");




*/


kill(all);
remvalue(all);
ratprint:false;
numer:true$
display2d:false$


/* 
converts complex number z = x*y*%i 
to the list in a draw format:  
[x,y] 
*/
d(z):=[float(realpart(z)), float(imagpart(z))]$

/* give Draw List from one point*/
dl(z):=points([d(z)])$




/* trigonometric functions in Maxima CAS use radians */
deg2rad(t):= float(t*2*%pi/360)$

GiveImplicit(a,b):=implicit( x^2/(a^2) + (y^2)/(b^2) = 1, x, -4,4, y, -4,4)$

GivePointOfEllipse(a,b, t):= a*cos(t) + b*sin(t)*%i$


/*

xc <- 1 # center x_c or h
yc <- 2 # y_c or k
a <- 5 # major axis length
b <- 2 # minor axis length
phi <- pi/3 # angle of major axis with x axis phi or tau

t <- seq(0, 2*pi, 0.01) 
x <- xc + a*cos(t)*cos(phi) - b*sin(t)*sin(phi)
y <- yc + a*cos(t)*sin(phi) + b*sin(t)*cos(phi)

<math>\mathbf{x} =\mathbf{x}_{\theta}(t) = a\cos\ t\cos\theta - b\sin\ t\sin\theta</math>

<math>\mathbf{y} =\mathbf{y}_{\theta}(t) = a\cos\ t\cos\theta + b\sin\ t\cos\theta</math>


https://stackoverflow.com/questions/65278354/how-to-draw-rotated-ellipse-in-maxima-cas/65294520#65294520
*/

GiveRotatedEllipse(a,b,theta, NumberOfPoints):=block(
	[x, y, zz, t , tmin, tmax, dt, c, s],
	zz:[],
	dt : 1/NumberOfPoints, 
 	tmin: 0, 
 	tmax: 2*%pi,
 	c:float(cos(theta)),
 	s:float(sin(theta)),
 	for t:tmin thru tmax step dt do(
 		x: a*cos(t)*c - b*sin(t)*s,
 		x: float(x), 
 		y: a*cos(t)*s + b*sin(t)*c,
 		y:float(y),
 		zz: cons([x,y],zz)
 	),
 	return (points(zz))
)$

GiveScaledRotatedEllipse(a,b, r,theta, NumberOfPoints):= GiveRotatedEllipse(r*a,r*b,theta, NumberOfPoints)$

GiveEllipseN(a,b,r,n,theta, NumberOfPoints):=GiveRotatedEllipse(a*(r^n),b*(r^n),n*theta, NumberOfPoints)$

Give_N(n):= GiveEllipseN(a,b,r,n,theta, NumberOfPoints)$

GiveEllipses(n):=block(
	[elipses],
	
	ellipses:makelist(i, i, 0, n, 1),
	ellipses:map(Give_N, ellipses),
	return(ellipses)
	


)$

/* 
scale ratio r = a'/a = b'/b

https://math.stackexchange.com/questions/3773593/given-ellipse-of-axes-a-and-b-find-axes-of-tangential-and-concentric-ellips
*/
GiveScaleRatio(a, b, theta):= block(
	[d, r], 
	d: (a/b - b/a)*sin(theta), 
	d:float(d),
	r: sqrt(1+d*d/4) - d/2,
	r:float(r),
	return(r)


)$


compile(all)$

/* compute */

/* angles fo trigonometric functions in radians */
angle: 15$
theta:deg2rad(angle) $  /* theta is the angle between    */
a: 5$
b: 4$
NumberOfPoints : 500$
r:GiveScaleRatio(a, b, theta)$  /* 0.942$ the (axis) scaled ratio r = a'/a = b'/b */


n:70;



ee:GiveEllipses(n)$



path:"~/Dokumenty/ellipse/scaled/s1/"$ /*  pwd, if empty then file is in a home dir , path should end with "/" */

/* draw it using draw package by */

 load(draw); 
/* if graphic  file is empty (= 0 bytes) then run draw2d command again */

 draw2d(
  user_preamble="set key top right; unset mouse",
  terminal  = 'svg,
  file_name = sconcat(path, string(a),"_",string(b), "_",string(theta), "_",string(r),"_", string(n)),
   title = "",  
  dimensions = [1500, 1000],
   axis_top         = false,
  axis_right       = false,
  axis_bottom         = false,
  axis_left       = false,
 ytics  = 'none,
 xtics  = 'none,
  proportional_axes = xy,
  line_width = 1,
  line_type = solid,
  
  fill_color = white,
  point_type=filled_circle,
  points_joined = true,
  point_size = 0.05,
  
    
  key = "",
  color = red,
  ee 
  )$