dommat:=Dom::Matrix():

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
rotate	:= proc( in_mat, x_psi, y_theta, z_phi) // PHI, THETA, PSI
	begin
		PHI		:= z_phi;
		THETA	:= y_theta;
		PSI		:= x_psi;
		mat_zroll :=
		dommat([
			[ cos( PHI),	sin( PHI),	0			],
			[ -sin( PHI),	cos( PHI),	0			],
			[ 0,			0,			1			]
		]):

		mat_yroll :=
		dommat([
			[ cos( THETA),	0,			-sin(THETA)	],
			[ 0,			1,			0			],	
			[ sin( THETA),	0,			cos( THETA)	]
		]):

		mat_xroll :=
		dommat([
			[ 1,			0,			0			],
			[ 0, 			cos( PSI),	sin( PSI)	],
			[ 0, 			-sin( PSI), cos( PSI)	]
		]):

		mat_afterx	:= mat_xroll * in_mat;
		mat_aftery	:= mat_yroll * mat_afterx;
		mat_afterz	:= mat_zroll * mat_aftery;
		return(mat_afterz)
	end_proc:
/*
rotate	:= proc( in_mat, x_psi, y_theta, z_phi)
begin
	PHI		:= z_phi;
	THETA	:= y_theta;
	PSI		:= x_psi;
	rotate_mat	:= dommat(
			[
					[
						cos(PHI) * cos(THETA),
        				sin(PHI) * cos(THETA),
     					-sin(THETA)
					],
					[
						- sin(PHI) * cos(PSI) + cos(PHI) * sin(PSI) * sin(THETA),
        				cos(PHI) * cos(PSI) + sin(PHI) * sin(PSI) * sin(THETA),
      					sin(PSI) * cos(THETA)
					],
					[
						sin(PHI) * sin(PSI) + cos(PHI) * cos(PSI) * sin(THETA),
        				- cos(PHI) * sin(PSI) + sin(PHI) * cos(PSI) * sin(THETA),
        				cos(PSI) * cos(THETA)
					]
			]):
	return( rotate_mat * in_mat);
end_proc:
*/



//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
cnv2p1	:= proc( in_x, in_y, in_z)
	begin
	// Initialize Matrix
	mat_before	:= dommat([
				[ in_x],
				[ in_y],
				[ in_z]
			]);
	// Rotate PI+arccos(1/3) around X axis
	mat_after	:= rotate( mat_before, PI+arccos(1/3), 0, 0);
	mat_before	:= mat_after;

	// Rotate PI/6 around Z axis
	mat_after	:= rotate( mat_before, 0, 0, PI/6);
	mat_before	:= mat_after;

	// Move
	factor_l	:= 1 / ( 2 * cos( PI/6));
	mat_addup	:= dommat([
					[-factor_l * sin( PI/6)],
					[ factor_l * cos( PI/6)],
					[ 0				]
					]);
	mat_after	:= mat_before + mat_addup;
	mat_before	:= mat_after;

	mat_addup	:= dommat([
					[ 0],
					[ 0],
					[ -sqrt(6)/3/3]
					]);

	mat_after	:= mat_before + mat_addup;

	in_x	:= op( mat_after, 1);
	in_y	:= op( mat_after, 2);
	in_z	:= op( mat_after, 3);

	return ( mat_after);
end_proc:

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
mcnv2p1	:= proc( mat_in)
	begin
	mat_in_x	:= op( mat_in, 1);
	mat_in_y	:= op( mat_in, 2);
	mat_in_z	:= op( mat_in, 3);
	return (
		cnv2p1( mat_in_x, mat_in_y, mat_in_z)
	);
end_proc:

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
cnv2p2	:= proc( in_x, in_y, in_z)
	begin

	// Get Plane p1
	mat_after	:= cnv2p1( in_x, in_y, in_z);
	mat_before	:= mat_after;

	// Rotate it 2*PI/3 around Z axis
	mat_after	:= rotate( mat_before, 0, 0, 2*PI/3);

	in_x	:= op( mat_after, 1);
	in_y	:= op( mat_after, 2);
	in_z	:= op( mat_after, 3);

	return ( mat_after);
end_proc:

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
mcnv2p2	:= proc( mat_in)
	begin
	mat_in_x	:= op( mat_in, 1);
	mat_in_y	:= op( mat_in, 2);
	mat_in_z	:= op( mat_in, 3);
	return (
		cnv2p2( mat_in_x, mat_in_y, mat_in_z)
	);
end_proc:

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
cnv2p3	:= proc( in_x, in_y, in_z)
	begin
	// Get Plane 1
	mat_after	:= cnv2p1( in_x, in_y, in_z);
	mat_before	:= mat_after;

	// Rotate it -2*PI/3 around Z axis
	mat_after	:= rotate( mat_before, 0, 0, -2*PI/3);

	in_x	:= op( mat_after, 1);
	in_y	:= op( mat_after, 2);
	in_z	:= op( mat_after, 3);

	return ( mat_after);
end_proc:

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
mcnv2p3	:= proc( mat_in)
	begin
	mat_in_x	:= op( mat_in, 1);
	mat_in_y	:= op( mat_in, 2);
	mat_in_z	:= op( mat_in, 3);
	return (
		cnv2p3( mat_in_x, mat_in_y, mat_in_z)
	);
end_proc:



//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
cnv2p4	:=	proc( in_x, in_y, in_z)
	begin
	// Initialize Matrix
	mat_before	:= dommat([
				[ in_x],
				[ in_y],
				[ in_z]
			]);

	// Rotate -PI/6 around Z axis
	mat_after	:= rotate( mat_before, 0, 0, -PI/6);
	mat_before	:= mat_after;

	// move the vector
	mat_addup	:= dommat( [ -sqrt(3)/6, -1/2, 0]);
	mat_after	:= mat_before + mat_addup;
	mat_before	:= mat_after;

	mat_addup	:= dommat([
					[ 0],
					[ 0],
					[ -sqrt(6)/3/3]
					]);

	mat_after	:= mat_before + mat_addup;

	in_x	:= op( mat_after, 1);
	in_y	:= op( mat_after, 2);
	in_z	:= op( mat_after, 3);

	return ( mat_after);
	end_proc:	



//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
mcnv2p4	:= proc( mat_in)
	begin
	mat_in_x	:= op( mat_in, 1);
	mat_in_y	:= op( mat_in, 2);
	mat_in_z	:= op( mat_in, 3);
	return (
		cnv2p4( mat_in_x, mat_in_y, mat_in_z)
	);
	end_proc:




//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
createTriangle	:=	proc( lstofmatrix, j)
	begin
	ret	:=
		plot::Polygon(
			plot::Point(
				op(lstofmatrix[i],1),
				op(lstofmatrix[i],2),
				op(lstofmatrix[i],3))$i=1..3,
				Closed=TRUE,
				Color=[1-j/15, j/15, 0.9],
				LineWidth=3,
				Filled=TRUE
		);
	return (ret);
	end_proc:



//---------------------------------------------------------------------------
// Definition of Function matlst2pointlst()
// Conversion from list of matrix to list of point( not plot::Pointlist)
//---------------------------------------------------------------------------
matlst2pointlst	:= proc( matlst_in, color_in)
	begin

	matlength	:= nops( matlst_in);
	counter 	:= 1;

	while counter <= matlength do
		ret_pointlst[ counter] := plot::Point(
									[
										op(matlst_in[counter],1),
										op(matlst_in[counter],2),
										op(matlst_in[counter],3)
									],
									Color = color_in,
									PointWidth = 50
								);
		counter := counter + 1;
	end_while;
	
	return ( ret_pointlst);

end_proc:

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
distance	:=	proc( x1, y1, z1, x2, y2, z2)
	begin
	dx	:= x2 - x1;
	dy	:= y2 - y1;
	dz	:= z2 - z1;
	return ( sqrt( dx^2 + dy^2 + dz^2));
end_proc:

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
intersection	:= proc( a, b, c)
	begin
	dist	:= distance( 0, 0, 0, a, b, c);
	ret_a	:= ( a / dist);
	ret_b	:= ( b / dist);
	ret_c	:= ( c / dist);
	mat_ret	:= dommat( [ [ret_a], [ret_b], [ret_c]]);
	return ( mat_ret);
end_proc:

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
mintersection	:= proc( mat_in)
	begin
		return ( intersection(
						op(mat_in,1),
						op(mat_in,2),
						op(mat_in,3)
					)
				);
end_proc:

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
mdotproduct	:= proc( mat_in1, mat_in2)
begin
	ret_val :=	op( mat_in1, 1) * op( mat_in2, 1)
			+	op( mat_in1, 2) * op( mat_in2, 2)
			+	op( mat_in1, 3) * op( mat_in2, 3);
	return ( ret_val)
end_proc:

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
vecLength	:= proc( mat_in)
begin
	dist	:= distance( 
					0,				0,				0,
					op( mat_in, 1),	op( mat_in, 2),	op( mat_in, 3)
				);
	return ( dist);
end_proc:

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
cnv2SphericalCoordinates2	:= proc( in_x, in_y, in_z)
begin
	f_alpha	:= 1;
	f_beta	:= 1;

	unitx	:=	dommat( [ [1], [0], [0]]);
	projxy	:=	dommat( [ [ in_x], [ in_y], [0]]);
	len_projxy	:= vecLength( projxy);

	len_unitx	:= 1; 
	cosalpha	:= mdotproduct( unitx, projxy) / ( len_projxy * len_unitx);
	ret_alpha	:= arccos( cosalpha);

	if float( in_y) < 0
		then	f_alpha	:= -1;
	end_if;

	rad			:= sqrt( in_x^2 + in_y^2 + in_z^2);	
	ret_beta	:= arccos( in_z / rad );	
	ret_mat		:= dommat(
						[
							[ ret_alpha	* f_alpha],
							[ ret_beta	],
							[ 1			]
						]
					);

	return( ret_mat);
end_proc:

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
mcnv2SphericalCoordinates2	:=	proc( in_mat)
begin
	return (
		cnv2SphericalCoordinates2(
			op( in_mat, 1), op( in_mat, 2), op( in_mat, 3)
	)
		);
end_proc:
	
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
cnv2SphericalCoordinates	:= proc( in_x, in_y, in_z)
begin

	rad		:= sqrt( in_x^2 + in_y^2 + in_z^2);
	if float( in_x) = 0
		then	ret_alpha	:= PI/2;
		else	ret_alpha	:= arctan( in_y / in_x);
	end_if;

	ret_beta	:= arccos( in_z / rad );


	if float( in_y) < 0 and float( in_x) < 0 
		then	ret_alpha := ret_alpha - PI / 2;
		else	if float( in_y) > 0 and float( in_x) < 0
					then	ret_alpha := ret_alpha + PI / 2;
				end_if;
  	end_if;

	ret_mat		:= dommat(
						[
							[ ret_alpha	],
							[ ret_beta	],
							[ 1						]
						]
					);

	return( ret_mat);
end_proc:

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
mcnv2SphericalCoordinates	:=	proc( in_mat)
begin
	return (
		cnv2SphericalCoordinates(
			op( in_mat, 1), op( in_mat, 2), op( in_mat, 3)
	)
		);
end_proc:

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
cnv2Mercator	:= proc( in_alpha, in_beta)
begin
	ret_x	:= in_alpha;
	ret_y	:= sin( PI / 2 - in_beta);
	ret_mat	:= dommat(
						[
							[ ret_x], [ ret_y], [ 0]
						]
				);
	return( ret_mat);
end_proc:

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
mcnv2Mercator	:= proc( in_mat)
begin
	return (
		cnv2Mercator( op( in_mat, 1), op( in_mat, 2))
		);
end_proc:

mcnvCartesian2Mercator	:= proc( in_mat)
begin

	ret_mat_after	:= mcnv2SphericalCoordinates2( in_mat);
	ret_mat_before	:= ret_mat_after;

	ret_mat_after	:= mcnv2Mercator( ret_mat_before);
	return( ret_mat_after);
end_proc: