connectorGravityTest.py
You can view and download this file on Github: connectorGravityTest.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test model for ObjectConnectorGravity, realizing gravitational forces between
5# two masses (used for astrospace or small-scale astronomical investigations, e.g., satellites);
6# in this case we simulate a small planetary system, initializing planets with orbital velocity
7#
8# Author: Johannes Gerstmayr
9# Date: 2022-01-30
10#
11# Copyright:This file is part of Exudyn. Exudyn is free software. You can redistribute it and/or modify it under the terms of the Exudyn license. See 'LICENSE.txt' for more details.
12#
13#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
14import sys
15sys.path.append('../TestModels')
16
17import exudyn as exu
18from exudyn.itemInterface import *
19from exudyn.utilities import *
20import numpy as np
21
22useGraphics = True #without test
23#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
24#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
25try: #only if called from test suite
26 from modelUnitTests import exudynTestGlobals #for globally storing test results
27 useGraphics = exudynTestGlobals.useGraphics
28except:
29 class ExudynTestGlobals:
30 pass
31 exudynTestGlobals = ExudynTestGlobals()
32#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
33
34#create an environment for mini example
35SC = exu.SystemContainer()
36mbs = SC.AddSystem()
37
38#some drawing properties, to see the objects ...
39massStar = 1e20
40massSatellite = 1e3
41mass = [massSatellite, massSatellite, massSatellite, massSatellite]
42sizeMass0 = 1e5 #just for drawing
43sizeMass = [2e4,2e4,2e4,2e4] #just for drawing
44rOrbit = [2e5, 4e5, 8e5, 10e5]
45vOrbitEps = [1.,1.,1.,1.] #factors to make non-circular orbits...
46
47background = GraphicsDataCheckerBoard(point=[0,0,-2*sizeMass0], size=2.1*max(rOrbit),
48 color=[0,0,0,1], alternatingColor=[0.05,0,0,1])
49
50oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
51 visualization=VObjectGround(graphicsData=[background])))
52# nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
53node0 = mbs.AddNode(NodePoint(referenceCoordinates = [0,0,0])) #planet
54gMass0 = GraphicsDataSphere(radius=1e5, color=color4blue, nTiles=64)
55oMassPoint0 = mbs.AddObject(MassPoint(nodeNumber = node0, physicsMass=massStar,
56 visualization=VMassPoint(graphicsData=[gMass0])))
57m0 = mbs.AddMarker(MarkerNodePosition(nodeNumber=node0))
58
59#create satellites:
60for i,r in enumerate(rOrbit):
61 G = 6.6743e-11
62 vOrbit = vOrbitEps[i]*np.sqrt(G*massStar/r) #orbital velocity, assumption of heavy star
63 #exu.Print('vOrbit'+str(i)+'=',vOrbit)
64 node1 = mbs.AddNode(NodePoint(referenceCoordinates = [r,0,0],
65 initialVelocities=[0,vOrbit,0])) #satellite
66
67 gMass1 = GraphicsDataSphere(radius=sizeMass[i], color=color4list[i], nTiles=24)
68
69 oMassPoint1 = mbs.AddObject(MassPoint(nodeNumber = node1, physicsMass=mass[i],
70 visualization=VMassPoint(graphicsData=[gMass1])))
71
72 m1 = mbs.AddMarker(MarkerNodePosition(nodeNumber=node1))
73
74 mbs.AddObject(ObjectConnectorGravity(markerNumbers=[m0,m1],
75 mass0 = massStar, mass1=mass[i]))
76
77#assemble and solve system for default parameters
78mbs.Assemble()
79simulationSettings = exu.SimulationSettings()
80
81tEnd = 1e6
82h = 1000
83simulationSettings.solutionSettings.writeSolutionToFile = False
84simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
85simulationSettings.timeIntegration.endTime = tEnd
86simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
87simulationSettings.timeIntegration.newton.useModifiedNewton = True
88
89simulationSettings.displayStatistics = True
90simulationSettings.timeIntegration.verboseMode = 1
91
92# SC.visualizationSettings.nodes.drawNodesAsPoint = False
93
94if useGraphics:
95 exu.StartRenderer() #start graphics visualization
96 mbs.WaitForUserToContinue() #wait for pressing SPACE bar to continue
97
98#start solver:
99# mbs.SolveDynamic(simulationSettings, solverType = exu.DynamicSolverType.TrapezoidalIndex2)
100#gives 7 digits of accuracy for tEnd=1e6, h=1e3:
101mbs.SolveDynamic(simulationSettings, solverType = exu.DynamicSolverType.RK67)
102
103if useGraphics:
104 SC.WaitForRenderEngineStopFlag()#wait for pressing 'Q' to quit
105 exu.StopRenderer() #safely close rendering window!
106
107#check result at default integration time
108#node1 is last node
109pos = mbs.GetNodeOutput(node1, exu.OutputVariableType.Position)
110
111exudynTestGlobals.testResult = pos[0] + pos[1] + pos[2]
112
113exu.Print("result for ObjectConnectorGravity =", exudynTestGlobals.testResult)
114#exudynTestGlobals.testResult = 1014867.2330320379