files=['Cl_pulse_1.hdf5','Ar_pulse_1.hdf5','Cl_pulse_2.hdf5','Ar_pulse_2.hdf5','Cl_pulse_3.hdf5','Ar_pulse_3.hdf5','Cl_pulse_4.hdf5','Ar_pulse_4.hdf5']
monolayer = 16 #number of atoms in a monolayer
ML_thickness = 0.136 # in nm
num_ML = 16 #number of monolayers in the substrate
hook=CleanVacuumRegion(fuzz_factor=1.15)
impact=[]
etch=[]
total_impacts=0

for file in files:
    sps=nlread(file,SurfaceProcessSimulation)[-1]
    traj=sps.trajectories()
    for i in range(0,len(traj)):
        conf=traj[i].lastImage()
        hook._apply(conf)
        sym=conf.symbols()
        count=sym.count('Si')
        impact.append(total_impacts+i)
        etch.append(ML_thickness * (num_ML-(count/monolayer))) # gives the thickness removed compared to pristine surace
    total_impacts+=len(traj)
    
# PlotModel with custom units.
model = Plot.PlotModel()
# Set axis label and scale
model.xAxis().setLabel('Impacts')
model.yAxis().setLabel('Etch yield [nm]')
#Line
line = Plot.Line(impact,etch)
# Add line.
model.addItem(line)
# Auto-adjust axis limits.
model.setLimits()
# Show plot and save afterwards to file and as raster.
Plot.show(model)
Plot.save(model, 'etch_plot.hdf5')