// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// %                                                                   %
// %   MechSim include file - version 0.3                              %
// %                                                                   %
// %   mechanics simulation include file                               %
// %                                                                   %
// %   written August 2002 - December 2002                             %
// %   by Christoph Hormann <chris_hormann@gmx.de>                     %
// %                                                                   %
// %   A short description of the parameters is given at the macro     %
// %   declarations.  Complete documentation is supplied in a          %
// %   separate file.                                                  %
// %                                                                   %
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// %                                                                   %
// %   This include file is part of the Sim-POV patch.  Like the       %
// %   patch it is subject to the conditions in povlegal.doc that      %
// %   comes with the Sim-POV package. No part of it may be used       %
// %   outside POV-Ray.                                                %
// %                                                                   %
// %   You may contact me if you have ideas for improvements or find   %
// %   errors.                                                         %
// %                                                                   %
// %   This include file like Sim-POV is in experimental state         %
// %   so syntax and used methods could possibly change in future      %
// %   versions.                                                       %
// %                                                                   %
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
#ifdef(MECHSIM_Inc_Temp)
// do nothing
#else
#declare MECHSIM_Inc_Temp = version;

#ifdef(View_POV_Include_Stack)
#   debug "including mechsim.inc\n"
#end

#version unofficial megapov 1.1;

#include "math.inc"
#include "transforms.inc"

#declare MSIM_INCLUDE_VERSION=0.3;
	
#declare MECHSIM_METHOD_EULER=1;
#declare MECHSIM_METHOD_HEUN=2;
#declare MECHSIM_METHOD_RUNGE_KUTTA4=3;
#declare MECHSIM_METHOD_GRADIENT=4;

#declare MECHSIM_COLLISION_NONE=0;
#declare MECHSIM_COLLISION_ALL=1;
#declare MECHSIM_COLLISION_GROUP=2;

#declare MECHSIM_ENV_METHOD_FORCE=1;
#declare MECHSIM_ENV_METHOD_IMPACT=2;

#declare MECHSIM_COLLISION_BOUNDING_NO=0;
#declare MECHSIM_COLLISION_BOUNDING_AUTO=1;
#declare MECHSIM_COLLISION_BOUNDING_BOX=2;
#declare MECHSIM_COLLISION_BOUNDING_HASH=3;

#declare No_Trans=transform{}

#ifndef (MSim_File_Decimals)
	#declare MSim_File_Decimals=12;
#end

#ifndef (MSim_Test_Render)
	#declare MSim_Test_Render=false;
#end

#ifndef (MSim_Tex_N)
	#declare MSim_Tex_N=
	texture {
		pigment {
			color rgb <0.8,0.8,1.0>
		}
		finish{
			diffuse 0.3
			ambient 0.0
			specular 0.6

			#if (!MSim_Test_Render)
				reflection {
					0.8
					metallic
				}
			#end

			conserve_energy
		}
	}
#end

#ifndef (MSim_Tex_C)
	#declare MSim_Tex_C=
	texture {
		pigment {
			color rgb <1, 0.5, 0.1>
		}
		finish{
			diffuse 0.3
			ambient 0.0
			specular 0.6

			#if (!MSim_Test_Render)
				reflection {
					0.8
					metallic
				}
			#end

			conserve_energy
		}
	}
#end

#ifndef (MSim_Tex_Mesh)
	#declare MSim_Tex_Mesh=
	texture {
		pigment {
			color rgb <0.4, 0.3, 1>
		}
		finish{
			diffuse 0.3
			ambient 0.0
			specular 0.6

			#if (!MSim_Test_Render)
				reflection {
					0.8
					metallic
				}
			#end

			conserve_energy
		}
	}
#end

#ifndef (fn_Stress_CM)
	#declare fn_Stress_CM=
	function {
		pigment {
			gradient x
			color_map {
				[ 0.00 color rgb<0.450, 0.750, 0.620> ]
				[ 0.25 color rgb<0.256, 0.430, 1.000> ]
				[ 0.50 color rgb<0.000, 0.000, 0.441> ]
				[ 0.75 color rgb<1.000, 0.456, 0.456> ]
				[ 1.00 color rgb<0.900, 0.800, 0.000> ]
			}
		}
	}
#end


// =======================================================================================
//  Vector_Function()
//
//  generates a pigment function from three float functions
//  representing the x, y and z directions
//
//  Parameters:
//      Fn_X, Fn_Y, Fn_Z - three user defined functions
//
//  The way this macro works makes it necessary to have a limited range for the
//  function values.  Vector_Function() is a wrapper macro for
//  Vector_Function_Range() using a default range of 1e6.  If necessary this can be
//  changed by declaring a different value for MSim_Fn_Range before including
//  'mechsim.inc'.
//
// =======================================================================================

#ifndef (MSim_Fn_Range)
	#declare MSim_Fn_Range=1e6;
#end

#macro Vector_Function_Range(Fn_X, Fn_Y, Fn_Z, Range)

	#local PigX=
	pigment {
		function { Fn_X(x, y, z)*(1/Range) }
		color_map {
			[0 color rgb 0 ]
			[1 color rgb 3*Range*x ]
		}
	}

	#local PigY=
	pigment {
		function { Fn_Y(x, y, z)*(1/Range) }
		color_map {
			[0 color rgb 0 ]
			[1 color rgb 3*Range*y ]
		}
	}

	#local PigZ=
	pigment {
		function { Fn_Z(x, y, z)*(1/Range) }
		color_map {
			[0 color rgb 0 ]
			[1 color rgb 3*Range*z ]
		}
	}

	function {
		pigment {
			average
			pigment_map {
				[1 PigX]
				[1 PigY]
				[1 PigZ]
			}
		}
	}

#end

#macro Vector_Function(Fn_X, Fn_Y, Fn_Z)

Vector_Function_Range(
	function { Fn_X(x, y, z) },
	function { Fn_Y(x, y, z) },
	function { Fn_Z(x, y, z) },
	MSim_Fn_Range
)

#end

// #######################################################################################
//   Object diplay macros
// #######################################################################################

// =======================================================================================
//  MechSim_Show_Objects()
//
//  generates objects representing parts of the simulation topology
//
//  Parameters:
//      Start_Mass_Index       - index of the first mass that should be displayed.
//                               Should be stored when the topology is generated.
//      Start_Connection_Index - index of the first connection
//      Start_Face_Index       - index of the first face
//      End_Mass_Index         - index of the last mass that should be displayed.
//                               if <0 all masses until the last are displayed.
//      End_Connection_Index   - index of the last connection that should be displayed.
//                               if <0 all connections until the last are displayed.
//      End_Face_Index         - index of the last face that should be displayed.
//                               if <0 all faces until the last are displayed.
//      Connect_Rad            - Radius to use for the connections, values <0 lead
//                               to using  0.4*mass radius.
//      Show_Faces             - Only generate a mesh from the faces instead of
//                               cylinders and spheres.
//      Stress_Fact            - if >0 the connections are textured to visualize
//                               the stress.  Values are scaled with this factor.
//      File_Name              - string, if length is >0 the objects are written to a
//                               file with that name in addition
//
// =======================================================================================
#macro MechSim_Show_Objects(Start_Mass_Index, Start_Connection_Index, Start_Face_Index,
                            End_Mass_Index, End_Connection_Index, End_Face_Index,
                            Connect_Rad, Show_Faces, Stress_Fact, File_Name)

object {

	#if (End_Mass_Index<0)
		#local End_Mass=mechsim:mass_count;
	#else
		#local End_Mass=End_Mass_Index;
	#end

	#if (End_Connection_Index<0)
		#local End_Connection=mechsim:connection_count;
	#else
		#local End_Connection=End_Connection_Index;
	#end

	#if (End_Face_Index<0)
		#local End_Face=mechsim:face_count;
	#else
		#local End_Face=End_Face_Index;
	#end

	#if (strlen(File_Name)=0)
		#local Write_File=false;
		#debug "Generating Objects..."
	#else
		#local Write_File=true;

		#fopen FILE File_Name write

		#if (Show_Faces)
			#write (FILE, "// simulation topology objects (",str(End_Face-Start_Face_Index,0,0)," faces)\n")
		#else
			#if (Connect_Rad != 0)
				#write (FILE, "// simulation topology objects (",str(End_Mass-Start_Mass_Index,0,0)," masses, ",
				                                                 str(End_Connection-Start_Connection_Index,0,0)," connections)\n")
			#else
				#write (FILE, "// simulation topology objects (",str(End_Mass-Start_Mass_Index,0,0)," masses)\n")
			#end
		#end
		#write (FILE, "// written by MechSim include file version ",str(MSIM_INCLUDE_VERSION,1,1),"\n\n")

		#if (Stress_Fact>0)
			#write (FILE, "#local Fin_Stress=\n")
			#write (FILE, "  finish {\n")
			#write (FILE, "    specular 0.3\n")
			#write (FILE, "    diffuse 0.7\n")
			#write (FILE, "  }\n\n")
		#end

		#debug concat("Generating Objects (writing to file ",File_Name,")...")

	#end

	#if (!Show_Faces)
		#if ((End_Mass>Start_Mass_Index) | (End_Connection>Start_Connection_Index))

			#if ((End_Mass>Start_Mass_Index) & (End_Connection>Start_Connection_Index))
				union {
					#if (Write_File) #write (FILE, "  union {\n") #end
			#else
				object {
					#if (Write_File) #write (FILE, "  object {\n") #end
			#end

			#if (End_Mass>Start_Mass_Index)

				#if (End_Mass>Start_Mass_Index+1)
					union {
						#if (Write_File) #write (FILE, "  union {\n") #end
				#else
					object {
						#if (Write_File) #write (FILE, "  object {\n") #end
				#end
										
				#local Cnt=Start_Mass_Index;

				#while (Cnt<End_Mass)
					sphere {
						mechsim:mass(Cnt):position,
	          mechsim:mass(Cnt):radius
	        }

					#if (Write_File)
						#write (FILE, "    sphere {<",vstr(3,mechsim:mass(Cnt):position,",",0,MSim_File_Decimals), ">, ",
							                             str(mechsim:mass(Cnt):radius,0,MSim_File_Decimals), " }\n")
					#end

					#local Cnt=Cnt+1;
				#end

						texture { MSim_Tex_N }
					}

				#if (Write_File)
					#write (FILE, "    texture { MSim_Tex_N }\n  }\n")
				#end
			#end

			#if (End_Connection>Start_Connection_Index)

				#if (Connect_Rad != 0)

					#if (End_Connection>Start_Connection_Index+1)
						union {
							#if (Write_File) #write (FILE, "  union {\n") #end
					#else
						object {
							#if (Write_File) #write (FILE, "  object {\n") #end
					#end

					#local Cnt=Start_Connection_Index;

					#while (Cnt<End_Connection)
						#local Dist = vlength(mechsim:mass(mechsim:connection(Cnt):index1):position-
							                    mechsim:mass(mechsim:connection(Cnt):index2):position);
						#if (Dist>0.0001)
							cylinder {
								mechsim:mass(mechsim:connection(Cnt):index1):position,
								mechsim:mass(mechsim:connection(Cnt):index2):position,

								#if (Write_File)
									#write (FILE, "    cylinder {<",vstr(3,mechsim:mass(mechsim:connection(Cnt):index1):position,",",0,MSim_File_Decimals), ">, <",
										                              vstr(3,mechsim:mass(mechsim:connection(Cnt):index2):position,",",0,MSim_File_Decimals), ">, ")
								#end

								#if (Connect_Rad<0)
									#local Radius=mechsim:mass(mechsim:connection(Cnt):index1):radius*0.4;
								#else
									#local Radius=Connect_Rad;
								#end

								Radius

								#if (Write_File)
									#write (FILE, str(Radius,0,MSim_File_Decimals), " ")
								#end

								#if (Stress_Fact>0)

									#local Stress=((Dist-mechsim:connection(Cnt):length)/mechsim:connection(Cnt):length)*Stress_Fact + 0.5;
									#local Stress=max(0.0, min(1.0, Stress));
									texture {
										pigment { color fn_Stress_CM(Stress, 0, 0) }
										finish {
											specular 0.3
											diffuse 0.7
										}
									}

									#if (Write_File)
										#write (FILE, "texture{pigment{ color rgb <",vstr(3,fn_Stress_CM(Stress,0,0),",",0,4), "> }finish{Fin_Stress}} ")
									#end

								#end
							}
							#if (Write_File)
								#write (FILE, "}\n")
							#end

						#end

						#local Cnt=Cnt+1;
					#end

					#if (Stress_Fact<=0)
						  texture { MSim_Tex_C }

						#if (Write_File)
							#write (FILE, "    texture { MSim_Tex_C }\n")
						#end
					#end
						}
					#if (Write_File)
						#write (FILE, "  }\n")
					#end
				#end
			#end
					}
			#if (Write_File)
				#write (FILE, "}\n")
			#end
		#end
	#else
		#if (mechsim:face_count)
			mesh {

				#if (Write_File)
					#write (FILE, "mesh {\n")
				#end

				#local Cnt=Start_Face_Index;

				#while (Cnt<End_Face)

					triangle { mechsim:mass(mechsim:face(Cnt):index1):position,
						         mechsim:mass(mechsim:face(Cnt):index2):position,
						         mechsim:mass(mechsim:face(Cnt):index3):position }

					#if (Write_File)
						#write (FILE, "triangle { <",vstr(3,mechsim:mass(mechsim:face(Cnt):index1):position,",",0,MSim_File_Decimals), ">, <",
							                           vstr(3,mechsim:mass(mechsim:face(Cnt):index2):position,",",0,MSim_File_Decimals), ">, <",
							                           vstr(3,mechsim:mass(mechsim:face(Cnt):index3):position,",",0,MSim_File_Decimals), "> }\n")
					#end

					#local Cnt=Cnt+1;
				#end

				texture { MSim_Tex_Mesh }
			}

			#if (Write_File)
				#write (FILE, "}\n")
			#end

		#end
	#end

	#if (Write_File)
		#fclose FILE
	#end

	}

	#debug "Done\n"

#end

// =======================================================================================
//  MechSim_Show_All_Objects()
//
//  like MechSim_Show_Objects() but showing all elements of the simulation
//
//  Parameters:
//      Connect_Rad - Radius to use for the connections, values <0 lead to using
//                    0.4*mass radius
//      Show_Faces  - Only generate a mesh from the faces instead of
//                    cylinders and spheres
//      Stress_Fact - if >0 the connections are textured to visualize the stress.
//                    values are scaled with this factor
//      File_Name   - string, if length is >0 the objects are written to a
//                    file with that name in addition
//
// =======================================================================================
#macro MechSim_Show_All_Objects(Connect_Rad, Show_Faces, Stress_Fact, File_Name)

	MechSim_Show_Objects(0,0,0, -1, -1, -1, Connect_Rad, Show_Faces, Stress_Fact, File_Name)

#end

// =======================================================================================
//  MechSim_Show_Patch() macro
//
//  generating a mesh corresponding to the topology generated by a
//  MechSim_Generate_Patch() macro
//
//  Parameters:
//      Start_Face_Index - index of the first face that's part of the patch.  Should
//                         be stored when the patch is generated
//      XSize, YSize     - Size of the patch, number of nodes in x- and y-direction
//      Smooth           - boolean value, if true normal vectors are generated
//      UV               - boolean value, if true UV-coordinates are generated
//                         (from 0 to 1 in x-direction)
//      Stress_Fact      - if >0 the mesh is textured to visualize the stress.
//                         values are scaled with this factor
//      File_Name        - string, if length is >0 the mesh is written to a file
//                         with that name in addition
//
// =======================================================================================
#macro MechSim_Show_Patch(Start_Face_Index, XSize, YSize, Smooth, UV, Stress_Fact, File_Name)

	#if ((mechsim:face_count-Start_Face_Index)>0)

		#if (strlen(File_Name)=0)
			#local Write_File=false;
			#debug "Generating Patch Triangle Mesh..."
		#else
			#local Write_File=true;

			#fopen FILE File_Name write

			#write (FILE, "// rectangular patch mesh (",str(XSize,0,0),"x",str(YSize,0,0)," vertices)\n")
			#write (FILE, "// written by MechSim include file version ",str(MSIM_INCLUDE_VERSION,1,1),"\n\n")
			#write (FILE, "mesh2 {\n\n")
			#debug concat("Generating Patch Triangle Mesh (writing to file ",File_Name,")...")

		#end

		#local Start_Mass_Index = mechsim:face(Start_Face_Index):index1;

		#if (abs(Stress_Fact)>0)

			#local Vertex_Mult=array[XSize*YSize]
			#local Vertex_Stress=array[XSize*YSize]
			#local Vertex_Text=array[XSize*YSize]

			#if (Write_File)
				#write (FILE, "  #local Vertex_Text=array[",str(XSize*YSize,0,0),"]\n")
				#write (FILE, "  #local Fin_Stress=\n")
				#write (FILE, "    finish {\n")
				#write (FILE, "      specular 0.3\n")
				#write (FILE, "      diffuse 0.7\n")
				#write (FILE, "    }\n\n")
			#end

			#local Cnt=0;
			#while (Cnt < XSize*YSize)
				#local Vertex_Mult[Cnt]=0;
				#local Vertex_Stress[Cnt]=0;
				#local Cnt=Cnt+1;
			#end

			#local Cnt=0;
			#while (Cnt < mechsim:connection_count)

				#if ((mechsim:connection(Cnt):index1 >= Start_Mass_Index) &
					   (mechsim:connection(Cnt):index2 >= Start_Mass_Index) &
					   (mechsim:connection(Cnt):index1 < Start_Mass_Index+(XSize*YSize)) &
					   (mechsim:connection(Cnt):index2 < Start_Mass_Index+(XSize*YSize)))

					#local Dist = vlength(mechsim:mass(mechsim:connection(Cnt):index1):position-
						                    mechsim:mass(mechsim:connection(Cnt):index2):position);

					#local Stress=((Dist-mechsim:connection(Cnt):length)/mechsim:connection(Cnt):length)*Stress_Fact + 0.5;
					#local Vertex_Stress[mechsim:connection(Cnt):index1-Start_Mass_Index]=
					       Vertex_Stress[mechsim:connection(Cnt):index1-Start_Mass_Index]+max(0.0, min(1.0, Stress));
					#local Vertex_Mult[mechsim:connection(Cnt):index1-Start_Mass_Index]=
					       Vertex_Mult[mechsim:connection(Cnt):index1-Start_Mass_Index]+1;
					#local Vertex_Stress[mechsim:connection(Cnt):index2-Start_Mass_Index]=
					       Vertex_Stress[mechsim:connection(Cnt):index2-Start_Mass_Index]+max(0.0, min(1.0, Stress));
					#local Vertex_Mult[mechsim:connection(Cnt):index2-Start_Mass_Index]=
					       Vertex_Mult[mechsim:connection(Cnt):index2-Start_Mass_Index]+1;

				#end
				#local Cnt=Cnt+1;
			#end

		#end

		#if (Write_File)
			#write (FILE, "  vertex_vectors {\n")
			#write (FILE, "    ",str(XSize*YSize,0,0))
		#end

		mesh2 {
			vertex_vectors {
				(XSize*YSize)

				#local Cnt=Start_Mass_Index;
				#while (Cnt < XSize*YSize + Start_Mass_Index)
					,
					mechsim:mass(Cnt):position

					#if (Write_File)
						#write (FILE, ",\n    <",vstr(3,mechsim:mass(Cnt):position,",",0,MSim_File_Decimals), ">")
					#end

					#local Cnt=Cnt+1;
				#end
			}

			#if (Write_File)
				#write (FILE, "\n  }\n")
			#end

			#if (Smooth)

				#if (Write_File)
					#write (FILE, "  normal_vectors {\n")
					#write (FILE, "    ",str(XSize*YSize,0,0))
				#end

				normal_vectors {
					(XSize*YSize)

					#local Cnt=Start_Mass_Index;
					#while (Cnt < XSize*YSize + Start_Mass_Index)
						,

						#local Pos=<mod(Cnt-Start_Mass_Index, XSize), div(Cnt-Start_Mass_Index, XSize)>;

						#if (Pos.x=0)
							#if (Pos.y=0)
								#local Norm=vcross(mechsim:mass(Cnt+1):position-mechsim:mass(Cnt):position,
									                 mechsim:mass(Cnt+XSize):position-mechsim:mass(Cnt):position);
							#else
								#if (Pos.y=YSize-1)
									#local Norm=vcross(mechsim:mass(Cnt+1):position-mechsim:mass(Cnt):position,
										                 mechsim:mass(Cnt):position-mechsim:mass(Cnt-XSize):position);
								#else /* Pos.y middle */
									#local Norm=vcross(mechsim:mass(Cnt+1):position-mechsim:mass(Cnt):position,
										                 mechsim:mass(Cnt+XSize):position-mechsim:mass(Cnt-XSize):position);
								#end
							#end
						#else
							#if (Pos.x=XSize-1)
								#if (Pos.y=0)
									#local Norm=vcross(mechsim:mass(Cnt):position-mechsim:mass(Cnt-1):position,
										                 mechsim:mass(Cnt+XSize):position-mechsim:mass(Cnt):position);
								#else
									#if (Pos.y=YSize-1)
										#local Norm=vcross(mechsim:mass(Cnt):position-mechsim:mass(Cnt-1):position,
											                 mechsim:mass(Cnt):position-mechsim:mass(Cnt-XSize):position);

									#else /* Pos.y middle */
										#local Norm=vcross(mechsim:mass(Cnt):position-mechsim:mass(Cnt-1):position,
											                 mechsim:mass(Cnt+XSize):position-mechsim:mass(Cnt-XSize):position);
									#end
								#end
							#else /* Pos.x middle */
								#if (Pos.y=0)
									#local Norm=vcross(mechsim:mass(Cnt+1):position-mechsim:mass(Cnt-1):position,
										                 mechsim:mass(Cnt+XSize):position-mechsim:mass(Cnt):position);
								#else
									#if (Pos.y=YSize-1)
										#local Norm=vcross(mechsim:mass(Cnt+1):position-mechsim:mass(Cnt-1):position,
											                 mechsim:mass(Cnt):position-mechsim:mass(Cnt-XSize):position);
									#else /* Pos.y and Pos.x middle */
										#local Norm=vcross(mechsim:mass(Cnt+1):position-mechsim:mass(Cnt-1):position,
											                 mechsim:mass(Cnt+XSize):position-mechsim:mass(Cnt-XSize):position);
									#end
								#end
							#end
						#end

						Norm

						#if (Write_File)
							#write (FILE, ",\n    <",vstr(3,Norm,",",0,MSim_File_Decimals), ">")
						#end

						#local Cnt=Cnt+1;
					#end
				}

				#if (Write_File)
					#write (FILE, "\n  }\n")
				#end

			#end

			#if (UV)

				#if (Write_File)
					#write (FILE, "  uv_vectors {\n")
					#write (FILE, "    ",str(XSize*YSize,0,0))
				#end

				uv_vectors {
					(XSize*YSize)

					#local Cnt=Start_Mass_Index;
					#while (Cnt < XSize*YSize + Start_Mass_Index)
						,
						#local UV_Pos=<mod(Cnt-Start_Mass_Index, XSize), div(Cnt-Start_Mass_Index, XSize)>/(XSize-1);

						UV_Pos

						#if (Write_File)
							#write (FILE, ",\n    <",vstr(2,UV_Pos,",",0,MSim_File_Decimals), ">")
						#end

						#local Cnt=Cnt+1;
					#end

				}

				#if (Write_File)
					#write (FILE, "\n  }\n")
				#end

			#end

			#if (abs(Stress_Fact)>0)

				#if (Write_File)
					#write (FILE, "  texture_list {\n")
					#write (FILE, "    ",str(XSize*YSize,0,0))
				#end

				texture_list {
					(XSize*YSize)

					#local Cnt=0;
					#while (Cnt < XSize*YSize)

						#local Vertex_Text[Cnt]=
						texture {
							pigment { color fn_Stress_CM(Vertex_Stress[Cnt]/Vertex_Mult[Cnt], 0, 0) }
							finish {
								specular 0.3
								diffuse 0.7
							}
						}

						,

						texture { Vertex_Text[Cnt] }

						#if (Write_File)
							#write (FILE, "\n    #local Vertex_Text[",str(Cnt,0,0),"]=texture{pigment{ color rgb <",
								            vstr(3, fn_Stress_CM(Vertex_Stress[Cnt]/Vertex_Mult[Cnt],0,0),",",0,4),
								            "> }finish{Fin_Stress}}\n")

							#write (FILE, "    , texture { Vertex_Text[",str(Cnt,0,0),"] }")
						#end

						#local Cnt=Cnt+1;
					#end
				}

				#if (Write_File)
					#write (FILE, "\n  }\n")
				#end

			#end

			#if (Write_File)
				#write (FILE, "  face_indices {\n")
				#write (FILE, "    ",str((XSize-1)*(YSize-1)*2,0,0))
			#end

			face_indices {
				(XSize-1)*(YSize-1)*2

				#local Cnt=Start_Face_Index;

				#while (Cnt< (XSize-1)*(YSize-1)*2 + Start_Face_Index)
					,
					<mechsim:face(Cnt):index1-Start_Mass_Index,
					 mechsim:face(Cnt):index2-Start_Mass_Index,
					 mechsim:face(Cnt):index3-Start_Mass_Index>
					
					#if (Write_File)
						#write (FILE, ",\n    <",vstr(3,<mechsim:face(Cnt):index1-Start_Mass_Index,
							                               mechsim:face(Cnt):index2-Start_Mass_Index,
							                               mechsim:face(Cnt):index3-Start_Mass_Index>,",",0,0), ">")
					#end

					#if (Stress_Fact>0)
						,
						mechsim:face(Cnt):index1-Start_Mass_Index,
						mechsim:face(Cnt):index2-Start_Mass_Index,
						mechsim:face(Cnt):index3-Start_Mass_Index

						#if (Write_File)
							#write (FILE, ", ",vstr(3,<mechsim:face(Cnt):index1-Start_Mass_Index,
								                         mechsim:face(Cnt):index2-Start_Mass_Index,
								                         mechsim:face(Cnt):index3-Start_Mass_Index>,",",0,0))
						#end

					#end

					#local Cnt=Cnt+1;
				#end
			}

			//texture { MSim_Tex_Mesh }
		}

		#if (Write_File)
			#write (FILE, "\n  }\n}\n")
			#fclose FILE
		#end

	#end

	#debug "Done\n"

#end

// =======================================================================================
//  MechSim_Show_Grid() macro
//
//  generating a mesh for the outer surface of a grid generated with the
//  MechSim_Generate_Grid() macro (fith 'Faces' parameter true)
//
//  Parameters:
//      Start_Face_Index - index of the first face that's part of the patch.  Should
//                         be stored when the patch is generated
//      XSize, YSize     - Size of the patch, number of nodes in x- and y-direction
//      Smooth           - boolean value, if true normal vectors are generated
//      UV               - boolean value, if true UV-coordinates are generated
//                         (from 0 to 1 in x-direction)
//      Stress_Fact      - if >0 the mesh is textured to visualize the stress.
//                         values are scaled with this factor
//      File_Name        - string, if length is >0 the mesh is written to a file
//                         with that name in addition
//
// =======================================================================================
#macro MechSim_Show_Grid(Start_Face_Index, XSize, YSize, ZSize, Smooth, UV, Stress_Fact, File_Name)

	#local FCount=(XSize-1)*(YSize-1)*4 + (XSize-1)*(ZSize-1)*4 + (ZSize-1)*(YSize-1)*4;
					
	#if ((mechsim:face_count-Start_Face_Index)>=FCount)

		#if (strlen(File_Name)=0)
			#local Write_File=false;
			#debug "Generating Grid Surface Triangle Mesh...\n"
		#else
			#local Write_File=true;

			#fopen FILE File_Name write

			#write (FILE, "// grid mesh (",str(XSize,0,0),"x",str(YSize,0,0),"x",str(ZSize,0,0)," vertices)\n")
			#write (FILE, "// written by MechSim include file version ",str(MSIM_INCLUDE_VERSION,1,1),"\n\n")
			#debug concat("Generating Grid Triangle Mesh (writing to file ",File_Name,")...\n")

		#end

		#local Start_Mass_Index = mechsim:face(Start_Face_Index):index1;
		#local Count=(XSize*YSize*ZSize);
		#local EPS=1e-6;

		#if (abs(Stress_Fact)>0)

			#local Vertex_Mult=array[XSize*YSize*ZSize]
			#local Vertex_Stress=array[XSize*YSize*ZSize]
			#local Vertex_Text=array[XSize*YSize*ZSize]

			#if (Write_File)
				#write (FILE, "  #local Vertex_Text=array[",str(XSize*YSize*ZSize,0,0),"]\n")
				#write (FILE, "  #local Fin_Stress=\n")
				#write (FILE, "    finish {\n")
				#write (FILE, "      specular 0.3\n")
				#write (FILE, "      diffuse 0.7\n")
				#write (FILE, "    }\n\n")
			#end

			#local Cnt=0;
			#while (Cnt < Count)
				#local Vertex_Mult[Cnt]=0;
				#local Vertex_Stress[Cnt]=0;
				#local Cnt=Cnt+1;
			#end

			#local Cnt=0;
			#while (Cnt < mechsim:connection_count)

				#if ((mechsim:connection(Cnt):index1 >= Start_Mass_Index) &
					   (mechsim:connection(Cnt):index2 >= Start_Mass_Index) &
					   (mechsim:connection(Cnt):index1 < Start_Mass_Index+Count) &
					   (mechsim:connection(Cnt):index2 < Start_Mass_Index+Count))

					#local Dist = vlength(mechsim:mass(mechsim:connection(Cnt):index1):position-
						                    mechsim:mass(mechsim:connection(Cnt):index2):position);

					#local Stress=((Dist-mechsim:connection(Cnt):length)/mechsim:connection(Cnt):length)*Stress_Fact + 0.5;
					#local Vertex_Stress[mechsim:connection(Cnt):index1-Start_Mass_Index]=
					       Vertex_Stress[mechsim:connection(Cnt):index1-Start_Mass_Index]+max(0.0, min(1.0, Stress));
					#local Vertex_Mult[mechsim:connection(Cnt):index1-Start_Mass_Index]=
					       Vertex_Mult[mechsim:connection(Cnt):index1-Start_Mass_Index]+1;
					#local Vertex_Stress[mechsim:connection(Cnt):index2-Start_Mass_Index]=
					       Vertex_Stress[mechsim:connection(Cnt):index2-Start_Mass_Index]+max(0.0, min(1.0, Stress));
					#local Vertex_Mult[mechsim:connection(Cnt):index2-Start_Mass_Index]=
					       Vertex_Mult[mechsim:connection(Cnt):index2-Start_Mass_Index]+1;

				#end
				#local Cnt=Cnt+1;
			#end

			#local Cnt=0;
			#while (Cnt < Count)

				#local Vertex_Text[Cnt]=
				texture {
					pigment { color fn_Stress_CM(Vertex_Stress[Cnt]/Vertex_Mult[Cnt], 0, 0) }
					finish {
						specular 0.3
						diffuse 0.7
					}
				}

				#if (Write_File)
					#write (FILE, "  #local Vertex_Text[",str(Cnt,0,0),"]=texture{pigment{ color rgb <",
						      vstr(3, fn_Stress_CM(Vertex_Stress[Cnt]/Vertex_Mult[Cnt],0,0),",",0,4),
						      "> }finish{Fin_Stress}}\n")
				#end

				#local Cnt=Cnt+1;
			#end
			
		#end

		mesh {
			
			#if (Write_File)
				#write (FILE, "\nmesh {\n")			
			#end
			
			#if (Smooth)
				
				#local Normals=array[Count][3]
				
				#local Cnt=0;
				#while (Cnt<Count)
					#local Normals[Cnt][0]=<0,0,0>;
					#local Normals[Cnt][1]=<0,0,0>;
					#local Normals[Cnt][2]=<0,0,0>;
					#local Cnt=Cnt+1;
				#end
					
				#local Cnt=Start_Face_Index;

				#while (Cnt<Start_Face_Index+FCount)
					
					#local Norm=vnormalize(vcross(mechsim:mass(mechsim:face(Cnt):index1):position-mechsim:mass(mechsim:face(Cnt):index2):position,
						                            mechsim:mass(mechsim:face(Cnt):index1):position-mechsim:mass(mechsim:face(Cnt):index3):position));
					
					#local Coord1=<mod(mechsim:face(Cnt):index1-Start_Mass_Index, XSize),
					               div(mod(mechsim:face(Cnt):index1-Start_Mass_Index, XSize*YSize), XSize),
					               div(mechsim:face(Cnt):index1-Start_Mass_Index, XSize*YSize)>;
					
					#local Coord2=<mod(mechsim:face(Cnt):index2-Start_Mass_Index, XSize),
					               div(mod(mechsim:face(Cnt):index2-Start_Mass_Index, XSize*YSize), XSize),
					               div(mechsim:face(Cnt):index2-Start_Mass_Index, XSize*YSize)>;
					
					#local Coord3=<mod(mechsim:face(Cnt):index3-Start_Mass_Index, XSize),
					               div(mod(mechsim:face(Cnt):index3-Start_Mass_Index, XSize*YSize), XSize),
					               div(mechsim:face(Cnt):index3-Start_Mass_Index, XSize*YSize)>;					

					#if ((abs(Coord1.x-Coord2.x)<EPS)&(abs(Coord2.x-Coord3.x)<EPS))
						// ------ x-negative/positive Face ------
						#local Normals[mechsim:face(Cnt):index1-Start_Mass_Index][0]=
						       Normals[mechsim:face(Cnt):index1-Start_Mass_Index][0] + Norm;
						#local Normals[mechsim:face(Cnt):index2-Start_Mass_Index][0]=
						       Normals[mechsim:face(Cnt):index2-Start_Mass_Index][0] + Norm;
						#local Normals[mechsim:face(Cnt):index3-Start_Mass_Index][0]=
						       Normals[mechsim:face(Cnt):index3-Start_Mass_Index][0] + Norm;
					#end

					#if ((abs(Coord1.y-Coord2.y)<EPS)&(abs(Coord2.y-Coord3.y)<EPS))
						// ------ y-negative/positive Face ------					
						#local Normals[mechsim:face(Cnt):index1-Start_Mass_Index][1]=
						       Normals[mechsim:face(Cnt):index1-Start_Mass_Index][1] + Norm;
						#local Normals[mechsim:face(Cnt):index2-Start_Mass_Index][1]=
						       Normals[mechsim:face(Cnt):index2-Start_Mass_Index][1] + Norm;
						#local Normals[mechsim:face(Cnt):index3-Start_Mass_Index][1]=
						       Normals[mechsim:face(Cnt):index3-Start_Mass_Index][1] + Norm;
					#end

					#if ((abs(Coord1.z-Coord2.z)<EPS)&(abs(Coord2.z-Coord3.z)<EPS))
						// ------ z-negative/positive Face ------	
						#local Normals[mechsim:face(Cnt):index1-Start_Mass_Index][2]=
						       Normals[mechsim:face(Cnt):index1-Start_Mass_Index][2] + Norm;
						#local Normals[mechsim:face(Cnt):index2-Start_Mass_Index][2]=
						       Normals[mechsim:face(Cnt):index2-Start_Mass_Index][2] + Norm;
						#local Normals[mechsim:face(Cnt):index3-Start_Mass_Index][2]=
						       Normals[mechsim:face(Cnt):index3-Start_Mass_Index][2] + Norm;
					#end
										
					#local Cnt=Cnt+1;
				#end
				
				#local Cnt=Start_Face_Index;

				#while (Cnt<Start_Face_Index+FCount)				
					
					#if ( // --- is x-Face? ---
						   (vlength(Normals[mechsim:face(Cnt):index1-Start_Mass_Index][0])>0.5)&
						   (vlength(Normals[mechsim:face(Cnt):index2-Start_Mass_Index][0])>0.5)&
						   (vlength(Normals[mechsim:face(Cnt):index3-Start_Mass_Index][0])>0.5)					
					    )
						smooth_triangle { 
							mechsim:mass(mechsim:face(Cnt):index1):position, 
							vnormalize(Normals[mechsim:face(Cnt):index1-Start_Mass_Index][0]),
							mechsim:mass(mechsim:face(Cnt):index2):position, 
							vnormalize(Normals[mechsim:face(Cnt):index2-Start_Mass_Index][0]),
							mechsim:mass(mechsim:face(Cnt):index3):position, 
							vnormalize(Normals[mechsim:face(Cnt):index3-Start_Mass_Index][0]) 
							#if (abs(Stress_Fact)>0)
								texture_list {
									Vertex_Text[mechsim:face(Cnt):index1]
									Vertex_Text[mechsim:face(Cnt):index2]
									Vertex_Text[mechsim:face(Cnt):index3]
								}									
							#end
						}
	
						#if (Write_File)
							#write (FILE, "smooth_triangle { <",vstr(3,mechsim:mass(mechsim:face(Cnt):index1):position,",",0,MSim_File_Decimals), ">, <",
								                                  vstr(3,vnormalize(Normals[mechsim:face(Cnt):index1-Start_Mass_Index][0]),",",0,MSim_File_Decimals), ">, <",
								                                  vstr(3,mechsim:mass(mechsim:face(Cnt):index2):position,",",0,MSim_File_Decimals), ">, <",
								                                  vstr(3,vnormalize(Normals[mechsim:face(Cnt):index2-Start_Mass_Index][0]),",",0,MSim_File_Decimals), ">, <",
								                                  vstr(3,mechsim:mass(mechsim:face(Cnt):index3):position,",",0,MSim_File_Decimals), ">, <",
								                                  vstr(3,vnormalize(Normals[mechsim:face(Cnt):index3-Start_Mass_Index][0]),",",0,MSim_File_Decimals), ">")
						
							#if (abs(Stress_Fact)>0)
								#write (FILE, "\n  texture_list { Vertex_Text[", 
									      str(mechsim:face(Cnt):index1,0,0), "] Vertex_Text[", 
									      str(mechsim:face(Cnt):index2,0,0), "] Vertex_Text[", 
									      str(mechsim:face(Cnt):index3,0,0), "] } }\n"
							#else
								#write (FILE, " }\n")
							#end						
						#end						
						
					#end

					#if ( // --- is y-Face? ---
						   (vlength(Normals[mechsim:face(Cnt):index1-Start_Mass_Index][1])>0.5)&
						   (vlength(Normals[mechsim:face(Cnt):index2-Start_Mass_Index][1])>0.5)&
						   (vlength(Normals[mechsim:face(Cnt):index3-Start_Mass_Index][1])>0.5)					
					    )
						smooth_triangle { 
							mechsim:mass(mechsim:face(Cnt):index1):position, 
							vnormalize(Normals[mechsim:face(Cnt):index1-Start_Mass_Index][1]),
							mechsim:mass(mechsim:face(Cnt):index2):position, 
							vnormalize(Normals[mechsim:face(Cnt):index2-Start_Mass_Index][1]),
							mechsim:mass(mechsim:face(Cnt):index3):position, 
							vnormalize(Normals[mechsim:face(Cnt):index3-Start_Mass_Index][1]) 
							#if (abs(Stress_Fact)>0)
								texture_list {
									Vertex_Text[mechsim:face(Cnt):index1]
									Vertex_Text[mechsim:face(Cnt):index2]
									Vertex_Text[mechsim:face(Cnt):index3]
								}									
							#end						
						}
						
						#if (Write_File)
							#write (FILE, "smooth_triangle { <",vstr(3,mechsim:mass(mechsim:face(Cnt):index1):position,",",0,MSim_File_Decimals), ">, <",
								                                  vstr(3,vnormalize(Normals[mechsim:face(Cnt):index1-Start_Mass_Index][1]),",",0,MSim_File_Decimals), ">, <",
								                                  vstr(3,mechsim:mass(mechsim:face(Cnt):index2):position,",",0,MSim_File_Decimals), ">, <",
								                                  vstr(3,vnormalize(Normals[mechsim:face(Cnt):index2-Start_Mass_Index][1]),",",0,MSim_File_Decimals), ">, <",
								                                  vstr(3,mechsim:mass(mechsim:face(Cnt):index3):position,",",0,MSim_File_Decimals), ">, <",
								                                  vstr(3,vnormalize(Normals[mechsim:face(Cnt):index3-Start_Mass_Index][1]),",",0,MSim_File_Decimals), ">")
							#if (abs(Stress_Fact)>0)
								#write (FILE, "\n  texture_list { Vertex_Text[", 
									      str(mechsim:face(Cnt):index1,0,0), "] Vertex_Text[", 
									      str(mechsim:face(Cnt):index2,0,0), "] Vertex_Text[", 
									      str(mechsim:face(Cnt):index3,0,0), "] } }\n"
							#else
								#write (FILE, " }\n")
							#end							
						#end
						
					#end
						
					#if ( // --- is z-Face? ---
						   (vlength(Normals[mechsim:face(Cnt):index1-Start_Mass_Index][2])>0.5)&
						   (vlength(Normals[mechsim:face(Cnt):index2-Start_Mass_Index][2])>0.5)&
						   (vlength(Normals[mechsim:face(Cnt):index3-Start_Mass_Index][2])>0.5)					
					    )
						smooth_triangle { 
							mechsim:mass(mechsim:face(Cnt):index1):position, 
							vnormalize(Normals[mechsim:face(Cnt):index1-Start_Mass_Index][2]),
							mechsim:mass(mechsim:face(Cnt):index2):position, 
							vnormalize(Normals[mechsim:face(Cnt):index2-Start_Mass_Index][2]),
							mechsim:mass(mechsim:face(Cnt):index3):position, 
							vnormalize(Normals[mechsim:face(Cnt):index3-Start_Mass_Index][2]) 
							#if (abs(Stress_Fact)>0)
								texture_list {
									Vertex_Text[mechsim:face(Cnt):index1]
									Vertex_Text[mechsim:face(Cnt):index2]
									Vertex_Text[mechsim:face(Cnt):index3]
								}									
							#end						
						}
						
						#if (Write_File)
							#write (FILE, "smooth_triangle { <",vstr(3,mechsim:mass(mechsim:face(Cnt):index1):position,",",0,MSim_File_Decimals), ">, <",
								                                  vstr(3,vnormalize(Normals[mechsim:face(Cnt):index1-Start_Mass_Index][2]),",",0,MSim_File_Decimals), ">, <",
								                                  vstr(3,mechsim:mass(mechsim:face(Cnt):index2):position,",",0,MSim_File_Decimals), ">, <",
								                                  vstr(3,vnormalize(Normals[mechsim:face(Cnt):index2-Start_Mass_Index][2]),",",0,MSim_File_Decimals), ">, <",
								                                  vstr(3,mechsim:mass(mechsim:face(Cnt):index3):position,",",0,MSim_File_Decimals), ">, <",
								                                  vstr(3,vnormalize(Normals[mechsim:face(Cnt):index3-Start_Mass_Index][2]),",",0,MSim_File_Decimals), ">")
							#if (abs(Stress_Fact)>0)
								#write (FILE, "\n  texture_list { Vertex_Text[", 
									      str(mechsim:face(Cnt):index1,0,0), "] Vertex_Text[", 
									      str(mechsim:face(Cnt):index2,0,0), "] Vertex_Text[", 
									      str(mechsim:face(Cnt):index3,0,0), "] } }\n"
							#else
								#write (FILE, " }\n")
							#end							
						#end						
						
					#end					
					#local Cnt=Cnt+1;
					
				#end
				
			#else // not Smooth
			
				#local Cnt=Start_Face_Index;

				#while (Cnt<Start_Face_Index+FCount)

					triangle { 
						mechsim:mass(mechsim:face(Cnt):index1):position,
						mechsim:mass(mechsim:face(Cnt):index2):position,
						mechsim:mass(mechsim:face(Cnt):index3):position 
						#if (abs(Stress_Fact)>0)
							texture_list {
								Vertex_Text[mechsim:face(Cnt):index1]
								Vertex_Text[mechsim:face(Cnt):index2]
								Vertex_Text[mechsim:face(Cnt):index3]
							}									
						#end						
					}

					#if (Write_File)
						#write (FILE, "triangle { <",vstr(3,mechsim:mass(mechsim:face(Cnt):index1):position,",",0,MSim_File_Decimals), ">, <",
							                           vstr(3,mechsim:mass(mechsim:face(Cnt):index2):position,",",0,MSim_File_Decimals), ">, <",
							                           vstr(3,mechsim:mass(mechsim:face(Cnt):index3):position,",",0,MSim_File_Decimals), ">")
						#if (abs(Stress_Fact)>0)
							#write (FILE, "\n  texture_list { Vertex_Text[", 
								      str(mechsim:face(Cnt):index1,0,0), "] Vertex_Text[", 
									    str(mechsim:face(Cnt):index2,0,0), "] Vertex_Text[", 
									    str(mechsim:face(Cnt):index3,0,0), "] } }\n"
						#else
							#write (FILE, " }\n")
						#end						
					#end

					#local Cnt=Cnt+1;
				#end
				
			#end

			texture { MSim_Tex_Mesh }
		}

		#if (Write_File)
			#write (FILE, "}\n")
		#end		
		
	#end

	#debug "Done\n"

#end					
					
// =======================================================================================
//  MechSim_Show_Sphere() macro
//
//  generating a mesh for the outer surface of a grid generated with the
//  MechSim_Generate_Sphere() macro (fith 'Faces' parameter true)
//
//  Parameters:
//      Start_Face_Index - index of the first face that's part of the patch.  Should
//                         be stored when the patch is generated
//      XSize, YSize     - Size of the patch, number of nodes in x- and y-direction
//      Smooth           - boolean value, if true normal vectors are generated
//      UV               - boolean value, if true UV-coordinates are generated
//                         (from 0 to 1 in x-direction)
//      Stress_Fact      - if >0 the mesh is textured to visualize the stress.
//                         values are scaled with this factor
//      File_Name        - string, if length is >0 the mesh is written to a file
//                         with that name in addition
//
// =======================================================================================
#macro MechSim_Show_Sphere(Start_Face_Index, Face_Count, Smooth, UV, Stress_Fact, File_Name)

	#if ((mechsim:face_count-Start_Face_Index)>=Face_Count)

		#if (strlen(File_Name)=0)
			#local Write_File=false;
			#debug "Generating Sphere Surface Triangle Mesh...\n"
		#else
			#local Write_File=true;

			#fopen FILE File_Name write

			#write (FILE, "// sphere mesh (",str(XSize,0,0),"x",str(YSize,0,0),"x",str(ZSize,0,0)," vertices)\n")
			#write (FILE, "// written by MechSim include file version ",str(MSIM_INCLUDE_VERSION,1,1),"\n\n")
			#write (FILE, "mesh {\n\n")
			#debug concat("Generating Sphere Triangle Mesh (writing to file ",File_Name,")...\n")

		#end

		mesh {

			#if (Write_File)
				#write (FILE, "mesh {\n")
			#end

			#if (Smooth)
							
				#local Center_Pos=mechsim:mass(mechsim:face(Start_Face_Index):index1-1):position;
				
				#local Cnt=Start_Face_Index;

				#while (Cnt<Start_Face_Index+Face_Count)

					smooth_triangle { mechsim:mass(mechsim:face(Cnt):index1):position, 
						                vnormalize(mechsim:mass(mechsim:face(Cnt):index1):position-Center_Pos),
						                mechsim:mass(mechsim:face(Cnt):index2):position, 
						                vnormalize(mechsim:mass(mechsim:face(Cnt):index2):position-Center_Pos),
						                mechsim:mass(mechsim:face(Cnt):index3):position, 
						                vnormalize(mechsim:mass(mechsim:face(Cnt):index3):position-Center_Pos) }								

					#if (Write_File)
						#write (FILE, "triangle { <",vstr(3,mechsim:mass(mechsim:face(Cnt):index1):position,",",0,MSim_File_Decimals), ">, <",
							                           vstr(3,mechsim:mass(mechsim:face(Cnt):index2):position,",",0,MSim_File_Decimals), ">, <",
							                           vstr(3,mechsim:mass(mechsim:face(Cnt):index3):position,",",0,MSim_File_Decimals), "> }\n")
					#end

					#local Cnt=Cnt+1;
				#end				
				
			#else // not Smooth
			
				#local Cnt=Start_Face_Index;

				#while (Cnt<Start_Face_Index+Face_Count)

					triangle { mechsim:mass(mechsim:face(Cnt):index1):position,
						         mechsim:mass(mechsim:face(Cnt):index2):position,
						         mechsim:mass(mechsim:face(Cnt):index3):position }

					#if (Write_File)
						#write (FILE, "triangle { <",vstr(3,mechsim:mass(mechsim:face(Cnt):index1):position,",",0,MSim_File_Decimals), ">, <",
							                           vstr(3,mechsim:mass(mechsim:face(Cnt):index2):position,",",0,MSim_File_Decimals), ">, <",
							                           vstr(3,mechsim:mass(mechsim:face(Cnt):index3):position,",",0,MSim_File_Decimals), "> }\n")
					#end

					#local Cnt=Cnt+1;
				#end
				
			#end

			texture { MSim_Tex_Mesh }
		}

		#if (Write_File)
			#write (FILE, "}\n")
		#end		
		
	#end

	#debug "Done\n"

#end										
			
// #######################################################################################
//   Grid topology generation macros
// #######################################################################################					
					
// =======================================================================================
//  MechSim_Generate_Grid_Fn() macro
//
//  topology generation macro building a 3D grid of masses connected and with optional
//  faces on the outside.  density, stiffness and damping can be controlled by functions
//
//  Parameters:
//      Velocity     - common velocity of all masses
//      Radius       - radius of the masses
//      fn_Density   - function defining the density of the masses.
//                     x, y and z in [0..1] range
//      fn_Stiffness - function defining the stiffness of the connections
//      fn_Damping   - function defining the damping of the connections
//      Faces        - Boolean value controlling whether faces are generated for the
//                     surface.
//      Cube_Scale   - Grid distance vector, defining the size of an elemetary
//                     cell in the grid.
//      Grid_Size    - Vector with number of masses in x, y and z direction
//      Transf       - Transform applied to the mass positions
//      Connect_Arr  - Array containing weights for connection stiffness and damping
//
// =======================================================================================
#macro MechSim_Generate_Grid_Fn(Velocity, Radius, fn_Density, fn_Stiffness, fn_Damping, Faces, Cube_Scale, Grid_Size, Transf, Connect_Arr)

	#debug "Generating Grid..."

	#local Vel=Velocity+<0,0,0>;
	#local CubeScale=Cube_Scale+<0,0,0>;
	#local GridSize=Grid_Size+<0,0,0>;

	#local EPS=0.001*CubeScale.x;

	#local StIdx=mechsim:mass_count;
	#local Count=(GridSize.x*GridSize.y*GridSize.z);
	#local Diagonal=dimension_size(Connect_Arr, 1);
	#if (Connect_Arr[0]=0)
		#local Diagonal=0;
	#end

	#local Cnt=0;

	#while (Cnt<Count)

		#local Coord=<mod(Cnt, GridSize.x),
		              div(mod(Cnt, GridSize.x*GridSize.y), GridSize.x),
		              div(Cnt, GridSize.x*GridSize.y)>;


		#local Coord2=<Coord.x*CubeScale.x,
		               Coord.y*CubeScale.y,
		               Coord.z*CubeScale.z>;

    #local CoordF=<Coord.x/(GridSize.x-1),
                   Coord.y/(GridSize.y-1),
                   Coord.z/(GridSize.z-1)>;

		mass { vtransform(Coord2, Transf), Vel, Radius density fn_Density(CoordF.x, CoordF.y, CoordF.z) }

		#local Cnt=Cnt+1;
	#end

	#if (Faces)
		
		#local Cnt=0;

		#while (Cnt<Count)

			#local Coord=<mod(Cnt, GridSize.x),
			              div(mod(Cnt, GridSize.x*GridSize.y), GridSize.x),
			              div(Cnt, GridSize.x*GridSize.y)>;

			#local UpB=GridSize-1-EPS;

			#if (((Coord.x<EPS)|(Coord.x>UpB.x))&
				    (Coord.y<UpB.y)&(Coord.z<UpB.z))  // ------ x-negative/positive Face ------

				face { StIdx+Cnt, StIdx+Cnt+GridSize.x, StIdx+Cnt+GridSize.x*GridSize.y+GridSize.x }
				face { StIdx+Cnt, StIdx+Cnt+GridSize.x*GridSize.y+GridSize.x, StIdx+Cnt+GridSize.x*GridSize.y }
			#end

			#if (((Coord.y<EPS)|(Coord.y>UpB.y))&
				    (Coord.x<UpB.x)&(Coord.z<UpB.z))  // ------ y-negative/positive Face ------

				face { StIdx+Cnt, StIdx+Cnt+1, StIdx+Cnt+GridSize.x*GridSize.y+1 }
				face { StIdx+Cnt, StIdx+Cnt+GridSize.x*GridSize.y+1, StIdx+Cnt+GridSize.x*GridSize.y}
			#end

			#if (((Coord.z<EPS)|(Coord.z>UpB.z))&
				    (Coord.x<UpB.x)&(Coord.y<UpB.y))  // ------ z-negative/positive Face ------

				face { StIdx+Cnt, StIdx+Cnt+1, StIdx+Cnt+GridSize.x+1 }
				face { StIdx+Cnt, StIdx+Cnt+GridSize.x+1, StIdx+Cnt+GridSize.x }
			#end

			#local Cnt=Cnt+1;

		#end
	#end


	#if (Diagonal>0)
		#local Cnt=0;

		#while (Cnt<Count)

			#local Coord=<mod(Cnt, GridSize.x),
			              div(mod(Cnt, GridSize.x*GridSize.y), GridSize.x),
			              div(Cnt, GridSize.x*GridSize.y)>;

			#local CoordF=<Coord.x/(GridSize.x-1),
			               Coord.y/(GridSize.y-1),
			               Coord.z/(GridSize.z-1)>;

			#local Stiffness=fn_Stiffness(CoordF.x, CoordF.y, CoordF.z);
			#local Damping=fn_Damping(CoordF.x, CoordF.y, CoordF.z);

			#if (Coord.x<GridSize.x-1)
				connection { StIdx+Cnt, StIdx+Cnt+1 stiffness Stiffness*Connect_Arr[0] damping Damping*Connect_Arr[0] }
				#if (Diagonal>1)
					#if (Coord.y<GridSize.y-1)
						connection { StIdx+Cnt, StIdx+Cnt+GridSize.x+1 stiffness Stiffness*Connect_Arr[1] damping Damping*Connect_Arr[1] }
						#if (Diagonal>2)
							connection { StIdx+Cnt+1, StIdx+Cnt+GridSize.x stiffness Stiffness*Connect_Arr[2] damping Damping*Connect_Arr[2] }
						#end
					#end
				#end
			#end
			#if (Coord.y<GridSize.y-1)
				connection { StIdx+Cnt, StIdx+Cnt+GridSize.x stiffness Stiffness*Connect_Arr[0] damping Damping*Connect_Arr[0] }
				#if (Diagonal>1)
					#if (Coord.z<GridSize.z-1)
						connection { StIdx+Cnt, StIdx+Cnt+GridSize.x*(GridSize.y+1) stiffness Stiffness*Connect_Arr[1] damping Damping*Connect_Arr[1] }
						#if (Diagonal>2)
							connection { StIdx+Cnt+GridSize.x, StIdx+Cnt+GridSize.x*GridSize.y stiffness Stiffness*Connect_Arr[2] damping Damping*Connect_Arr[2] }
						#end
					#end
				#end
			#end
			#if (Coord.z<GridSize.z-1)
				connection { StIdx+Cnt, StIdx+Cnt+GridSize.x*GridSize.y stiffness Stiffness*Connect_Arr[0] damping Damping*Connect_Arr[0] }
				#if (Diagonal>1)
					#if (Coord.x<GridSize.x-1)
						connection { StIdx+Cnt, StIdx+Cnt+GridSize.x*GridSize.y+1 stiffness Stiffness*Connect_Arr[1] damping Damping*Connect_Arr[1] }
						#if (Diagonal>2)
							connection { StIdx+Cnt+1, StIdx+Cnt+GridSize.x*GridSize.y stiffness Stiffness*Connect_Arr[2] damping Damping*Connect_Arr[2] }
						#end
					#end
				#end
			#end

			#if (Diagonal>3)
				#if ((Coord.x<GridSize.x-1) &
					   (Coord.y<GridSize.y-1) &
					   (Coord.z<GridSize.z-1))

					connection { StIdx+Cnt, StIdx+Cnt+GridSize.x*(GridSize.y+1)+1 stiffness Stiffness*Connect_Arr[3] damping Damping*Connect_Arr[3] }
					connection { StIdx+Cnt+1, StIdx+Cnt+GridSize.x*(GridSize.y+1) stiffness Stiffness*Connect_Arr[3] damping Damping*Connect_Arr[3] }
					connection { StIdx+Cnt+GridSize.x, StIdx+Cnt+GridSize.x*GridSize.y+1 stiffness Stiffness*Connect_Arr[3] damping Damping*Connect_Arr[3] }
					connection { StIdx+Cnt+GridSize.x+1, StIdx+Cnt+GridSize.x*GridSize.y stiffness Stiffness*Connect_Arr[3] damping Damping*Connect_Arr[3] }
				#end
			#end

			#local Cnt=Cnt+1;
		#end
	#end

	#debug "Done\n"

#end

// =======================================================================================
//  MechSim_Generate_Grid() macro
//
//  topology generation macro building a 3D grid of masses connected and with optional
//  faces on the outside.  density, stiffness and damping as fixed values
//
//  Parameters:
//      Velocity    - common velocity of all masses
//      Radius      - radius of the masses
//      Density     - density of the masses
//      Stiffness   - stiffness of the connections
//      Damping     - damping of the connections
//      Faces       - Boolean value controlling whether faces are generated for the
//                    surface.
//      Cube_Scale  - Grid distance vector, defining the size of an elemetary
//                    cell in the grid.
//      Grid_Size   - Vector with number of masses in x, y and z direction
//      Transf      - Transform applied to the mass positions
//      Connect_Arr - Array containing weights for connection stiffness and damping
//
// =======================================================================================
#macro MechSim_Generate_Grid(Velocity, Radius, Density, Stiffness, Damping, Faces, Cube_Scale, Grid_Size, Transf, Connect_Arr)

	#local fn_Density=function(x, y, z) {Density}
	#local fn_Stiffness=function(x, y, z) {Stiffness}
	#local fn_Damping=function(x, y, z) {Damping}

	MechSim_Generate_Grid_Fn(Velocity, Radius, fn_Density, fn_Stiffness, fn_Damping, Faces, Cube_Scale, Grid_Size, Transf, Connect_Arr)

#end

// =======================================================================================
//  MechSim_Generate_Grid_Std() macro
//
//  topology generation macro building a 3D grid of masses connected and with optional
//  faces on the outside.  density, stiffness and damping as fixed values.
//  Instead of the connection weight array a Diagonal parameter is used
//
//  Parameters:
//      Velocity    - common velocity of all masses
//      Radius      - radius of the masses
//      Density     - density of the masses
//      Stiffness   - stiffness of the connections
//      Damping     - damping of the connections
//      Faces       - Boolean value controlling whether faces are generated for the
//                    surface.
//      Cube_Scale  - Grid distance vector, defining the size of an elemetary
//                    cell in the grid.
//      Grid_Size   - Vector with number of masses in x, y and z direction
//      Transf      - Transform applied to the mass positions
//      Diagonal   - controling which diagonal connections are generated in the grid.
//
// =======================================================================================
#macro MechSim_Generate_Grid_Std(Velocity, Radius, Density, Stiffness, Damping, Faces, Cube_Scale, Grid_Size, Transf, Diagonal)

	#local fn_Density=function(x, y, z) {Density}
	#local fn_Stiffness=function(x, y, z) {Stiffness}
	#local fn_Damping=function(x, y, z) {Damping}

	#local Connect_Arr=array[max(1, Diagonal)]

	#local Cnt=0;
	#while (Cnt<max(1, Diagonal))
		#local Connect_Arr[Cnt]=1;
		#local Cnt=Cnt+1;
	#end

	#if (Diagonal=0)
		#local Connect_Arr[0]=0;
	#end

	MechSim_Generate_Grid_Fn(Velocity, Radius, fn_Density, fn_Stiffness, fn_Damping, Faces, Cube_Scale, Grid_Size, Transf, Connect_Arr)

#end

// =======================================================================================
//  MechSim_Generate_Box() macro
//
//  topology generation macro building a 3D grid of masses connected and with optional
//  faces on the outside.  Mass of the whole body given.
//
//  Parameters:
//      Velocity   - common velocity of all masses
//      Radius     - radius of the masses
//      Mass       - mass of the whole body
//      Stiffness  - stiffness of the connections
//      Damping    - damping of the connections
//      Faces      - Boolean value controlling whether faces are generated for the
//                   surface.
//      Cube_Scale - Grid distance vector, defining the size of an elemetary
//                   cell in the grid.
//      Grid_Size  - Vector with number of masses in x, y and z direction
//      Transf     - Transform applied to the mass positions
//      Diagonal   - controling which diagonal connections are generated in the grid.
//
// =======================================================================================
#macro MechSim_Generate_Box(Velocity, Radius, Mass, Stiffness, Damping, Faces, Cube_Scale, Grid_Size, Transf, Diagonal)

	#local Density=(Mass/(Grid_Size.x*Grid_Size.y*Grid_Size.z))/((4/3)*pi*pow(Radius,3));

	#local fn_Density=function(x, y, z) {Density}
	#local fn_Stiffness=function(x, y, z) {Stiffness}
	#local fn_Damping=function(x, y, z) {Damping}

	#local Connect_Arr=array[max(1, Diagonal)]

	#local Cnt=0;
	#while (Cnt<max(1, Diagonal))
		#local Connect_Arr[Cnt]=1;
		#local Cnt=Cnt+1;
	#end

	#if (Diagonal=0)
		#local Connect_Arr[0]=0;
	#end

	MechSim_Generate_Grid_Fn(Velocity, Radius, fn_Density, fn_Stiffness, fn_Damping, Faces, Cube_Scale, Grid_Size, Transf, Connect_Arr)

#end

// #######################################################################################
//   Patch topology generation macros
// #######################################################################################					

// =======================================================================================
//  MechSim_Generate_Patch() macro
//
//  topology generation macro building a rectangular patch of masses connected and with
//  optional faces.
//
//  Parameters:
//      Velocity    - common velocity of all masses
//      Radius      - radius of the masses
//      Density     - density of the masses
//      Stiffness   - stiffness of the connections
//      Damping     - damping of the connections
//      Faces       - Boolean value controlling whether faces are generated for the
//                    surface.
//      Rect_Scale  - 2D Grid distance vector, defining the size of an elemetary
//                    rectangle in the grid.
//      Grid_Size   - 2D Vector with number of masses in x and y direction
//      Transf      - Transform applied to the mass positions
//      fn_Fixed    - function controlling which nodes are fixed
//      Connect_Arr - Array containing weights for connection stiffness and damping
//
// =======================================================================================
#macro MechSim_Generate_Patch(Velocity, Radius, Density, Stiffness, Damping, Faces, Rect_Scale, Grid_Size, Transf, fn_Fixed, Connect_Arr)

	#debug "Generating Patch..."

	#local Vel=Velocity+<0,0,0>;
	#local CubeScale=Rect_Scale+<0,0,0>;
	#local GridSize=Grid_Size+<0,0,0>;

	#local EPS=0.001*CubeScale.x;

	#local StIdx=mechsim:mass_count;
	#local Count=(GridSize.x*GridSize.y);
	#local Con_Cnt=dimension_size(Connect_Arr, 1);

	#local Cnt=0;

	#while (Cnt<Count)

		#local Coord=<mod(Cnt, GridSize.x),
		              div(Cnt, GridSize.x)>;


		#local Coord2=<Coord.x*CubeScale.x,
		               Coord.y*CubeScale.y, 0>;


		#local CoordF=<Coord.x/(GridSize.x-1),
		               Coord.y/(GridSize.y-1)>;

		mass {
			vtransform(Coord2, Transf), Vel, Radius density Density

			#if ( fn_Fixed(CoordF.x, CoordF.y, 0)>0 )
				fixed on
			#end
		}

		#local Cnt=Cnt+1;
	#end

	#if (Faces)

		#local Cnt=0;

		#while (Cnt<Count)

			#local Coord=<mod(Cnt, GridSize.x),
			              div(Cnt, GridSize.x)>;

			#local UpB=GridSize-1-EPS;

			#if ((Coord.x<UpB.x)&(Coord.y<UpB.y))

				face { StIdx+Cnt, StIdx+Cnt+1, StIdx+Cnt+GridSize.x+1 }
				face { StIdx+Cnt, StIdx+Cnt+GridSize.x, StIdx+Cnt+GridSize.x+1 }

			#end

			#local Cnt=Cnt+1;

		#end
	#end


	#if (Con_Cnt>0)

		#local Cnt=0;

		#while (Cnt<Count)

			#local Coord=<mod(Cnt, GridSize.x),
			              div(Cnt, GridSize.x)>;

			#local UpB=GridSize-1-EPS;

			#if (Coord.x<UpB.x)
				connection { StIdx+Cnt, StIdx+Cnt+1 stiffness Stiffness*Connect_Arr[0] damping Damping*Connect_Arr[0] length CubeScale.x }
			#end
			#if (Coord.y<UpB.y)
				connection { StIdx+Cnt, StIdx+Cnt+GridSize.x stiffness Stiffness*Connect_Arr[0] damping Damping*Connect_Arr[0] length CubeScale.y }
			#end

			#if (Con_Cnt>1)
				#if ((Coord.x<UpB.x)&(Coord.y<UpB.y))
					connection { StIdx+Cnt+GridSize.x, StIdx+Cnt+1 stiffness Stiffness*Connect_Arr[1] damping Damping*Connect_Arr[1]
						           length sqrt(CubeScale.x*CubeScale.x + CubeScale.y*CubeScale.y) }
					connection { StIdx+Cnt, StIdx+Cnt+GridSize.x+1 stiffness Stiffness*Connect_Arr[1] damping Damping*Connect_Arr[1]
						           length sqrt(CubeScale.x*CubeScale.x + CubeScale.y*CubeScale.y) }
				#end

				#if (Con_Cnt>2)

					#local UpB2=GridSize-2-EPS;

					#if (Coord.x<UpB2.x)
						connection { StIdx+Cnt, StIdx+Cnt+2 stiffness Stiffness*Connect_Arr[2] damping Damping*Connect_Arr[2]
							           length 2*CubeScale.x }
					#end
					#if (Coord.y<UpB2.y)
						connection { StIdx+Cnt, StIdx+Cnt+GridSize.x*2 stiffness Stiffness*Connect_Arr[2] damping Damping*Connect_Arr[2]
							           length 2*CubeScale.y }
					#end

					#if (Con_Cnt>3)
						#if ((Coord.x<UpB2.x)&(Coord.y<UpB.y))
							connection { StIdx+Cnt+GridSize.x, StIdx+Cnt+2 stiffness Stiffness*Connect_Arr[3] damping Damping*Connect_Arr[3]
								           length sqrt(CubeScale.x*CubeScale.x*4 + CubeScale.y*CubeScale.y) }
							connection { StIdx+Cnt, StIdx+Cnt+GridSize.x+2 stiffness Stiffness*Connect_Arr[3] damping Damping*Connect_Arr[3]
								           length sqrt(CubeScale.x*CubeScale.x*4 + CubeScale.y*CubeScale.y) }
						#end

						#if ((Coord.x<UpB.x)&(Coord.y<UpB2.y))
							connection { StIdx+Cnt+GridSize.x*2, StIdx+Cnt+1 stiffness Stiffness*Connect_Arr[3] damping Damping*Connect_Arr[3]
								           length sqrt(CubeScale.x*CubeScale.x + CubeScale.y*CubeScale.y*4) }
							connection { StIdx+Cnt, StIdx+Cnt+GridSize.x*2+1 stiffness Stiffness*Connect_Arr[3] damping Damping*Connect_Arr[3]
								           length sqrt(CubeScale.x*CubeScale.x + CubeScale.y*CubeScale.y*4) }
						#end

						#if (Con_Cnt>4)
							#if ((Coord.x<UpB2.x)&(Coord.y<UpB2.y))
								connection { StIdx+Cnt+GridSize.x*2, StIdx+Cnt+2 stiffness Stiffness*Connect_Arr[4] damping Damping*Connect_Arr[4]
									           length 2*sqrt(CubeScale.x*CubeScale.x + CubeScale.y*CubeScale.y) }
								connection { StIdx+Cnt, StIdx+Cnt+GridSize.x*2+2 stiffness Stiffness*Connect_Arr[4] damping Damping*Connect_Arr[4]
									           length 2*sqrt(CubeScale.x*CubeScale.x + CubeScale.y*CubeScale.y) }
							#end

							#if (Con_Cnt>5)

								#local UpB3=GridSize-3-EPS;

								#if ((Coord.x<UpB3.x)&(Coord.y<UpB3.y))
									connection { StIdx+Cnt, StIdx+Cnt+3 stiffness Stiffness*Connect_Arr[5] damping Damping*Connect_Arr[5]
										           length 3*CubeScale.x }
									connection { StIdx+Cnt, StIdx+Cnt+GridSize.x*3 stiffness Stiffness*Connect_Arr[5] damping Damping*Connect_Arr[5]
										           length 3*CubeScale.y }
								#end
							#end /* Connect>5 */
						#end /* Connect>4 */
					#end /* Connect>3 */
				#end /* Connect>2 */
			#end /* Connect>1 */

			#local Cnt=Cnt+1;

		#end

	#end /* Connect>0 */

	#debug "Done\n"

#end

// =======================================================================================
//  MechSim_Generate_Patch_Std() macro
//
//  topology generation macro building a rectangular patch of masses connected and with
//  optional faces. No nodes fixed and same stiffness and damping for all connections.
//
//  Parameters:
//      Velocity    - common velocity of all masses
//      Radius      - radius of the masses
//      Density     - density of the masses
//      Stiffness   - stiffness of the connections
//      Damping     - damping of the connections
//      Faces       - Boolean value controlling whether faces are generated for the
//                    surface.
//      Rect_Scale  - 2D Grid distance vector, defining the size of an elemetary
//                    rectangle in the grid.
//      Grid_Size   - 2D Vector with number of masses in x and y direction
//      Transf      - Transform applied to the mass positions
//      Connect     - Which connections to generate
//
// =======================================================================================
#macro MechSim_Generate_Patch_Std(Velocity, Radius, Density, Stiffness, Damping, Faces, Rect_Scale, Grid_Size, Transf, Connect)

	#local fn_Fixed=function(x, y, z) { -1 }
	#local Connect_Arr=array[Connect]

	#local Cnt=0;
	#while (Cnt<Connect)
		#local Connect_Arr[Cnt]=1;
		#local Cnt=Cnt+1;
	#end

	MechSim_Generate_Patch(Velocity, Radius, Density, Stiffness, Damping, Faces, Rect_Scale, Grid_Size, Transf, fn_Fixed, Connect_Arr)

#end

// #######################################################################################
//   Line topology generation macros
// #######################################################################################					
					
// =======================================================================================
//  MechSim_Generate_Line() macro
//
//  topology generation macro building a straight line of masses connected.
//
//  Parameters:
//      Velocity    - common velocity of all masses
//      Radius      - radius of the masses
//      Density     - density of the masses
//      Stiffness   - stiffness of the connections
//      Damping     - damping of the connections
//      Spacing     - Distance between the masses
//      Count       - Number of masses
//      Direction   - Direction vector of the line
//      Transf      - Transform applied to the mass positions
//      fn_Fixed    - function controlling which nodes are fixed
//      Connect_Arr - Array containing weights for connection stiffness and damping
//
// =======================================================================================
#macro MechSim_Generate_Line(Velocity, Radius, Density, Stiffness, Damping, Spacing, Count, Direction, Transf, fn_Fixed, Connect_Arr)

	#debug "Generating Line..."

	#local Vel=Velocity+<0,0,0>;
	#local Dir=vnormalize(Direction);

	#local StIdx=mechsim:mass_count;
	#local Con_Cnt=dimension_size(Connect_Arr, 1);

	#if (Con_Cnt>0)

		#local Cnt=0;

		#while (Cnt<Count)

			mass {
				vtransform(Dir*Cnt*Spacing, Transf), Velocity, Radius density Density

				#if ( fn_Fixed(Cnt/(Count-1), 0, 0)>0 )
					fixed on
				#end
			}

			#local Cnt=Cnt+1;
		#end


		#if (Con_Cnt>0)

			#local Cnt=0;

			#while (Cnt<Count-1)

				#local BCnt=0;
				#while ((BCnt<Con_Cnt) & ((Cnt+BCnt)<Count-1))

					connection {
						StIdx+Cnt,
						StIdx+Cnt+BCnt+1
						stiffness Stiffness*Connect_Arr[BCnt]
						damping Damping*Connect_Arr[BCnt]
						length Spacing*(BCnt+1)
					}

					#local BCnt=BCnt+1;
				#end

				#local Cnt=Cnt+1;
			#end

		#end

	#end

	#debug "Done\n"

#end

// =======================================================================================
//  MechSim_Generate_Line_Std() macro
//
//  topology generation macro building a straight line of masses connected.
//  Only simple connections and not fixed
//
//  Parameters:
//      Velocity    - common velocity of all masses
//      Radius      - radius of the masses
//      Density     - density of the masses
//      Stiffness   - stiffness of the connections
//      Damping     - damping of the connections
//      Spacing     - Distance between the masses
//      Count       - Number of masses
//      Direction   - Direction vector of the line
//      Transf      - Transform applied to the mass positions
//
// =======================================================================================
#macro MechSim_Generate_Line_Std(Velocity, Radius, Density, Stiffness, Damping, Spacing, Count, Direction, Transf)

	#local fn_Fixed=function(x, y, z) { -1 }
	#local Connect_Arr=array[1]
	#local Connect_Arr[0]=1;

	MechSim_Generate_Line(Velocity, Radius, Density, Stiffness, Damping, Spacing, Count, Direction, Transf, fn_Fixed, Connect_Arr)

#end
			
// #######################################################################################
//   Ball topology generation macros
// #######################################################################################					
		
// --- Icosahedron data ---					
					
#declare Icosahedron_Vtx = array[12]

#local ZPos=1-4/(5+sqrt(5));
#local Rad=sqrt(1-ZPos*ZPos);
#local Cnt=1;
#while (Cnt<6)   
	#declare Icosahedron_Vtx[  Cnt]=<Rad*cos(Cnt*2*pi/5),Rad*sin(Cnt*2*pi/5),-ZPos>;
	#declare Icosahedron_Vtx[5+Cnt]=<Rad*cos((Cnt-7.5)*2*pi/5),Rad*sin((Cnt-7.5)*2*pi/5),ZPos>;
	#declare Cnt=Cnt+1;
#end  
#declare Icosahedron_Vtx[0]=<0,0,-1>;
#declare Icosahedron_Vtx[11]=<0,0,1>;

#declare Icosahedron_Face = array[20]

#declare Icosahedron_Face[ 0]=< 0, 1, 2>;
#declare Icosahedron_Face[ 1]=< 0, 2, 3>;
#declare Icosahedron_Face[ 2]=< 0, 3, 4>;
#declare Icosahedron_Face[ 3]=< 0, 4, 5>;
#declare Icosahedron_Face[ 4]=< 0, 5, 1>;
#declare Icosahedron_Face[ 5]=< 1, 2, 9>;
#declare Icosahedron_Face[ 6]=< 2, 3,10>;
#declare Icosahedron_Face[ 7]=< 3, 4, 6>;
#declare Icosahedron_Face[ 8]=< 4, 5, 7>;
#declare Icosahedron_Face[ 9]=< 5, 1, 8>;
#declare Icosahedron_Face[10]=< 1, 8, 9>;
#declare Icosahedron_Face[11]=< 9,10, 2>;
#declare Icosahedron_Face[12]=<10, 6, 3>;
#declare Icosahedron_Face[13]=< 6, 7, 4>;
#declare Icosahedron_Face[14]=< 7, 8, 5>;
#declare Icosahedron_Face[15]=<11, 6, 7>;
#declare Icosahedron_Face[16]=<11, 7, 8>;
#declare Icosahedron_Face[17]=<11, 8, 9>;
#declare Icosahedron_Face[18]=<11, 9,10>;
#declare Icosahedron_Face[19]=<11,10, 6>;

// --- Builds a triangle mesh from Vertex and Face arrays ---						
					
#macro MechSim_Ico_Build_Mesh(Arr_V, Arr_F)

	#local Cnt_F=dimension_size(Arr_F, 1);
					
	#local Cnt=0;
	#while (Cnt<Cnt_F)	
	
		triangle { Arr_V[Arr_F[Cnt].x], Arr_V[Arr_F[Cnt].y], Arr_V[Arr_F[Cnt].z] }
		
		#local Cnt=Cnt+1;
	#end

#end

// --- Builds a mechsim topology from Vertex and Face arrays ---						
					
#macro MechSim_Ico_Build_Grid(Cnt_E, Arr_V, Arr_F, Velocity, Radius, Density, Stiffness, Damping, Density2, Stiffness2, Damping2, Faces, BRadius, Transf)

	#local StIdx=mechsim:mass_count;
	#local Cnt_V=dimension_size(Arr_V, 1);
	#local Cnt_F=dimension_size(Arr_F, 1);		
					
	#local Arr_E=array[Cnt_E]
	#local Is=array[3];
	
	#local FCnt=0;
	#local EnCnt=0;

	#while (FCnt<Cnt_F)	
	
		#local Cnt=0;
		#local Is[0]=false;
		#local Is[1]=false;
		#local Is[2]=false;
		
		#while (Cnt<EnCnt)
				
				#if (((Arr_E[Cnt].x=Arr_F[FCnt].x)&(Arr_E[Cnt].y=Arr_F[FCnt].y))|
					   ((Arr_E[Cnt].y=Arr_F[FCnt].x)&(Arr_E[Cnt].x=Arr_F[FCnt].y)))
					#local Is[0]=true;
				#end
				#if (((Arr_E[Cnt].x=Arr_F[FCnt].y)&(Arr_E[Cnt].y=Arr_F[FCnt].z))|
					   ((Arr_E[Cnt].y=Arr_F[FCnt].y)&(Arr_E[Cnt].x=Arr_F[FCnt].z)))
					#local Is[1]=true;
				#end
				#if (((Arr_E[Cnt].x=Arr_F[FCnt].z)&(Arr_E[Cnt].y=Arr_F[FCnt].x))|
					   ((Arr_E[Cnt].y=Arr_F[FCnt].z)&(Arr_E[Cnt].x=Arr_F[FCnt].x)))
					#local Is[2]=true;
				#end
				
			#local Cnt=Cnt+1;
		#end
		
		#if (Is[0]=false)
			#local Arr_E[EnCnt]=<Arr_F[FCnt].x, Arr_F[FCnt].y>;
			#local EnCnt=EnCnt+1;
		#end
		#if (Is[1]=false)
			#local Arr_E[EnCnt]=<Arr_F[FCnt].y, Arr_F[FCnt].z>;
			#local EnCnt=EnCnt+1;
		#end
		#if (Is[2]=false)
			#local Arr_E[EnCnt]=<Arr_F[FCnt].z, Arr_F[FCnt].x>;
			#local EnCnt=EnCnt+1;
		#end
		
		#local FCnt=FCnt+1;
	#end
	
	// center mass
	#local Center_Pos=vtransform(<0,0,0>, Transf);
	mass {
		Center_Pos, Velocity, Radius density Density2
	}						
					
	#local Cnt=0;
	#while (Cnt<Cnt_V)		
	
		#local MPosition=vtransform(Arr_V[Cnt]*BRadius, Transf);
		
		mass {
			MPosition, Velocity, Radius density Density
		}				
		
		connection { 
			StIdx, StIdx+Cnt+1 stiffness Stiffness2 damping Damping2
			length vlength(MPosition-Center_Pos) 
		}			
		
		#local Cnt=Cnt+1;
	#end					
					
	#local Cnt=0;
	#while (Cnt<EnCnt)
		
		#local MPosition1=vtransform(Arr_V[Arr_E[Cnt].x]*BRadius, Transf);
		#local MPosition2=vtransform(Arr_V[Arr_E[Cnt].y]*BRadius, Transf);
		
		connection { 
			StIdx+Arr_E[Cnt].x+1, StIdx+Arr_E[Cnt].y+1 stiffness Stiffness damping Damping
			length vlength(MPosition1-MPosition2)
		}			
	
		#local Cnt=Cnt+1;
	#end
	
	#if (Faces)
		#local Cnt=0;
		#while (Cnt<Cnt_F)			
			
			face { 				
				StIdx+Arr_F[Cnt].x+1,
				StIdx+Arr_F[Cnt].y+1,
				StIdx+Arr_F[Cnt].z+1
			}
			
			#local Cnt=Cnt+1;
		#end		
	#end
		
#end

//  Pnt6 indices:
//
//      0
//     / \
//    5---3
//   / \ / \
//  2---4---1

#macro MechSim_Subdivide_Ico_Rec(Depth, Cnt_V, Cnt_F, Cnt_E, Arr_V, Arr_F, Velocity, Radius, Density, Stiffness, Damping, Density2, Stiffness2, Damping2, Faces, BRadius, Transf)

	#if (Depth>0)

		#local Cnt_Vn=Cnt_V+Cnt_E;
		#local Cnt_Fn=4*Cnt_F;
		#local Cnt_En=2*Cnt_E + 3*Cnt_F;
		
		#local Arr_Vn=array[Cnt_Vn]
		#local Arr_Fn=array[Cnt_Fn]
		#local Pnt6=array[6]
		
		#local FCnt=0;
		#local FnCnt=0;	
		#local VnCnt=0;		
		
		#while (FCnt<Cnt_F)			
						
			#local Cnt=0;
			#while (Cnt<6)	
				#local Pnt6[Cnt]=-1;
				#local Cnt=Cnt+1;
			#end
			
			#local Mid1 = vnormalize(Arr_V[Arr_F[FCnt].x]+Arr_V[Arr_F[FCnt].y]);
			#local Mid2 = vnormalize(Arr_V[Arr_F[FCnt].y]+Arr_V[Arr_F[FCnt].z]);
			#local Mid3 = vnormalize(Arr_V[Arr_F[FCnt].z]+Arr_V[Arr_F[FCnt].x]);	
			
			#local Cnt=0;
			#while (Cnt<VnCnt)
				
				#if (vlength(Arr_Vn[Cnt]-Arr_V[Arr_F[FCnt].x])<1e-6)					
					#local Pnt6[0]=Cnt;
				#end
				#if (vlength(Arr_Vn[Cnt]-Arr_V[Arr_F[FCnt].y])<1e-6)	
					#local Pnt6[1]=Cnt;
				#end
				#if (vlength(Arr_Vn[Cnt]-Arr_V[Arr_F[FCnt].z])<1e-6)	
					#local Pnt6[2]=Cnt;
				#end
				#if (vlength(Arr_Vn[Cnt]-Mid1)<1e-6)	
					#local Pnt6[3]=Cnt;
				#end
				#if (vlength(Arr_Vn[Cnt]-Mid2)<1e-6)	
					#local Pnt6[4]=Cnt;
				#end
				#if (vlength(Arr_Vn[Cnt]-Mid3)<1e-6)	
					#local Pnt6[5]=Cnt;
				#end
				
				#local Cnt=Cnt+1;
			#end
			
			#if (Pnt6[0]<0)
				#local Arr_Vn[VnCnt]=Arr_V[Arr_F[FCnt].x];
				#local Pnt6[0]=VnCnt;
				#local VnCnt=VnCnt+1;	
			#end
			#if (Pnt6[1]<0)
				#local Arr_Vn[VnCnt]=Arr_V[Arr_F[FCnt].y];
				#local Pnt6[1]=VnCnt;
				#local VnCnt=VnCnt+1;	
			#end
			#if (Pnt6[2]<0)
				#local Arr_Vn[VnCnt]=Arr_V[Arr_F[FCnt].z];
				#local Pnt6[2]=VnCnt;
				#local VnCnt=VnCnt+1;	
			#end
			#if (Pnt6[3]<0)
				#local Arr_Vn[VnCnt]=Mid1;
				#local Pnt6[3]=VnCnt;
				#local VnCnt=VnCnt+1;	
			#end
			#if (Pnt6[4]<0)
				#local Arr_Vn[VnCnt]=Mid2;
				#local Pnt6[4]=VnCnt;
				#local VnCnt=VnCnt+1;	
			#end
			#if (Pnt6[5]<0)
				#local Arr_Vn[VnCnt]=Mid3;
				#local Pnt6[5]=VnCnt;
				#local VnCnt=VnCnt+1;	
			#end
			
			#local Arr_Fn[FnCnt]=<Pnt6[0], Pnt6[3], Pnt6[5]>;
			#local FnCnt=FnCnt+1;
			#local Arr_Fn[FnCnt]=<Pnt6[1], Pnt6[4], Pnt6[3]>;
			#local FnCnt=FnCnt+1;
			#local Arr_Fn[FnCnt]=<Pnt6[2], Pnt6[5], Pnt6[4]>;
			#local FnCnt=FnCnt+1;
			#local Arr_Fn[FnCnt]=<Pnt6[3], Pnt6[4], Pnt6[5]>;
			#local FnCnt=FnCnt+1;
				
			#local FCnt=FCnt+1;
		#end
		
		#if (Depth>1)
	
			MechSim_Subdivide_Ico_Rec(Depth-1, Cnt_Vn, Cnt_Fn, Cnt_En, Arr_Vn, Arr_Fn, Velocity, Radius, Density, Stiffness, Damping, Density2, Stiffness2, Damping2, Faces, BRadius, Transf)
		
		#else
			
			MechSim_Ico_Build_Grid(Cnt_En, Arr_Vn, Arr_Fn, Velocity, Radius, Density, Stiffness, Damping, Density2, Stiffness2, Damping2, Faces, BRadius, Transf)
			
		#end
		
	#end

#end
		
#macro MechSim_Subdivide_Ico(Depth, Velocity, Radius, Density, Stiffness, Damping, Density2, Stiffness2, Damping2, Faces, Division, BRadius, Transf)

  #local Cnt_V=12;
  #local Cnt_F=20;
  #local Cnt_E=30;

  #local Arr_V=Icosahedron_Vtx;
  #local Arr_F=Icosahedron_Face;

  MechSim_Subdivide_Ico_Rec(Division, Cnt_V, Cnt_F, Cnt_E, Arr_V, Arr_F, Velocity, Radius, Density, Stiffness, Damping, Density2, Stiffness2, Damping2, Faces, BRadius, Transf)

#end					
					
// =======================================================================================
//  MechSim_Generate_Sphere() macro
//
//  topology generation macro building a sphere.
//  Only simple connections and not fixed
//
//  Parameters:
//      Velocity    - common velocity of all masses
//      Radius      - radius of the masses
//      Density     - density of the surface masses
//      Stiffness   - stiffness of the surface connections
//      Damping     - damping of the surface connections
//      Density2    - density of the center mass
//      Stiffness2  - stiffness of the center connections
//      Damping2    - damping of the center connections
//      Division    - Icosahedron subdivision depth
//      BRadius     - Radius of the ball
//      Transf      - Transform applied to the mass positions
//
// =======================================================================================
#macro MechSim_Generate_Sphere(Velocity, Radius, Density, Stiffness, Damping, Density2, Stiffness2, Damping2, Faces, Division, BRadius, Transf)

	#debug "Generating Sphere...\n"

	#local Vel=Velocity+<0,0,0>;

	 MechSim_Subdivide_Ico(Division, Vel, Radius, Density, Stiffness, Damping, Density2, Stiffness2, Damping2, Faces, Division, BRadius, Transf)

	#debug "Done\n"
					
#end					


#version MECHSIM_Inc_Temp;
#end
