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