Commit 57ca9f5d authored by jclamber's avatar jclamber
Browse files

compuet gas life of particles

git-svn-id: http://svn.oamp.fr/repos/uns_projects/trunk@89 f264a43e-d52d-4b82-913b-c2bd5215a18a
parent 5ac6575a
#!/usr/bin/python
#
# Find out gas life of particles
# MANDATORY
from py_unsio import * # import py_unsio package (UNSIO)
import numpy as np # arrays are treated as numpy arrays
import math
# -----------------------------------------------------
# Def "run" is a function which loads Ids of particles
def run(comp):
ok,time=uns.getValueF("time")
#get Ids
ok,id = uns.getArrayI(comp,"id")
#return time,id
return time,id[id.argsort()]
#return status, pos,vel,mass,rho,hsml,ienerg,Id
# <codecell>
# -----------------------------------------------------
# M A I N P R O G R A M
# -----------------------------------------------------
simname="mdf001"
#simname="/home/jcl/list.test"
components="gas"
times="all"
verbose=False
# Create a UNSIO object
uns = CunsIn(simname,components,times,verbose)
bits="I" # select properties, "" means all
# <codecell>
comp="gas"
first=True
outf="/home/jcl/gas_life.txt" # output file name
fo = open(outf, "w") # open file for writing
while (uns.nextFrame(bits)): # loop while there is something to read
time,ids=run(comp) # call "run" functionwhich returns time and array of Ids
print time
if first: #The first time....
time_last=time # save current time
ids_ref=np.copy(ids) # make a copy of ids array into ids_ref array
#print ids_ref
first=False
else: #all the other time
#print ids
print ids_ref.size,ids.size
inter=np.in1d(ids_ref,ids) # compute intersection of ird_ref and new id array
#print inter
idfinal=ids_ref[inter==False] # create an array if particles which have vanished into stars
#print "idfinal:",idfinal
fo.write("%f %d\n"%(time_last,idfinal.size)) # write into file, time and #partciles which ended their gas life
fo.write("\n".join(map(str, idfinal))) # write into file, particles ID which ended their gas life
ids_ref=ids_ref[inter==True] # create new ID_ref, ie particles whixh are still gas
if (idfinal.size>1):
fo.write("\n")
print "time [%f] ids_ref [%d] end of life[%d]:"%(time_last,ids_ref.size,idfinal.size)
time_last=time # save last time
fo.close() # close file
print "First must be false, here first = ", first
# <codecell>
print ids.size,ids_ref.size, idfinal.size, time_last
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment