Calculation in Marshall Mix Design

Hello,
I'm trying to calculate some data from my YADE Marshall Mix Design. I want to get value such as Bulk Volume, Bulk Density, VMA, % Air Void. But the result that I get seem to be negative and incorrect. So, I am confuse on my calculation part whether I set something wrong or my formula is inccorect?

##The Following is my Code##
import random
import math
from yade import geom, pack, utils, plot, ymport
import pandas as pd

# Define material properties
youngModulus = 1e7
poissonRatio = 0.25
density = 2700
frictionAngle =0.7

# Create material
material = O.materials.append(FrictMat(young=youngModulus, poisson=poissonRatio, density=density))
material1 = O.materials.append(FrictMat(young=youngModulus, poisson=poissonRatio, density=density, frictionAngle=frictionAngle))

# Define cylinder with funnel parameters
center = (0, 0, 0)
diameter = 0.102
height = 0.18

# create cylindrical body with radius 0.102 m and height 0.064 m
cylinder = geom.facetCylinder(center=center, radius=diameter/2, height=height, segmentsNumber=80, wallMask=6)

# assign material to each body in the cylinder
for body in cylinder:
    body.bodyMat = material

# add cylinder to simulation
O.bodies.append(cylinder)

# Define cylinder with funnel parameters
center1 = (0,0,height/2)
dBunker = 0.4
dOutput = 0.102
hBunker = 0
hOutput = 0.15
hPipe = 0

# create funnel as a bunker with diameter 0.102 m, height 0.064 m
funnel = geom.facetBunker(center=center1, dBunker=dBunker, dOutput=dOutput, hBunker=hBunker,hOutput=hOutput, hPipe=hPipe, segmentsNumber=80, wallMask=4)

# assign material to each body in the funnel
for body in funnel:
    body.bodyMat = material

# add funnel to simulation
O.bodies.append(funnel)

# define sphere parameters and number of spheres
rMean1 = (0.0125+0.019)/4
rRelFuzz1 = (0.019-0.0125)/4/rMean1
num1 = 13
rMean2 = (0.0095+0.0125)/4
rRelFuzz2 = (0.0125-0.0095)/4/rMean2
num2 = 51
rMean3 = (0.00475+0.0095)/4
rRelFuzz3 = (0.0095-0.00475)/4/rMean3
num3 = 563
rMean4 = (0.00236+0.00475)/4
rRelFuzz4 = (0.00475-0.00236)/4/rMean4
num4 = 5101
rMean5 = (0.00118+0.00236)/4
rRelFuzz5 = (0.00236-0.00118)/4/rMean4
num5 = 18369

#create empty sphere packing
sp = pack.SpherePack()

# generate randomly sphere
sp.makeCloud((-dBunker/4,-dBunker/4,1.3*height),(dBunker/4,dBunker/4,2*height), rMean = rMean1, rRelFuzz = rRelFuzz1, num = num1)
sp.makeCloud((-dBunker/4,-dBunker/4,1.3*height),(dBunker/4,dBunker/4,2*height), rMean = rMean2, rRelFuzz = rRelFuzz2, num = num2)
sp.makeCloud((-dBunker/4,-dBunker/4,1.3*height),(dBunker/4,dBunker/4,2*height), rMean = rMean3, rRelFuzz = rRelFuzz3, num = num3)
sp.makeCloud((-dBunker/4,-dBunker/4,1.3*height),(dBunker/4,dBunker/4,2*height), rMean = rMean4, rRelFuzz = rRelFuzz4, num = num4)
sp.makeCloud((-dBunker/4,-dBunker/4,1.3*height),(dBunker/4,dBunker/4,2*height), rMean = rMean5, rRelFuzz = rRelFuzz5, num = num5)

# add the sphere pack to the simulation
sp.toSimulation(material = material1)

for body in O.bodies:
   if not isinstance(body.shape, Sphere):
       continue
   if body.shape.radius >= rMean1 :
       body.shape.color = (0,0,1) #blue
   if body.shape.radius <= rMean1 and body.shape.radius > rMean2:
       body.shape.color = (0,0,1) #blue
   if body.shape.radius <= rMean2 and body.shape.radius > rMean3:
       body.shape.color = (1,0,0) #red
   if body.shape.radius <= rMean3 and body.shape.radius > rMean4:
       body.shape.color = (0,1,0) #green
   if body.shape.radius <= rMean4 and body.shape.radius > rMean5:
       body.shape.color = (1,1,0) #yellow
   if body.shape.radius <= rMean5 :
       body.shape.color = (1,0,1) #magenta

O.engines = [
        ForceResetter(),
        # sphere, facet, wall
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb(), Bo1_Wall_Aabb()]),
        InteractionLoop(
                # the loading plate is a wall, we need to handle sphere+sphere, sphere+facet, sphere+wall
                [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=(0, 0, -1000), damping=0.3),
        # the label creates an automatic variable referring to this engine
        # we use it below to change its attributes from the functions called
        PyRunner(command='checkUnbalanced()', realPeriod=2, label='checker'),
]
O.dt = PWaveTimeStep()

# enable energy tracking; any simulation parts supporting it
# can create and update arbitrary energy types, which can be
# accessed as O.energy['energyName'] subsequently
O.trackEnergy = True

# the following checkUnbalanced, unloadPlate and stopUnloading functions are all called by the 'checker'
# (the last engine) one after another; this sequence defines progression of different stages of the
# simulation, as each of the functions, when the condition is satisfied, updates 'checker' to call
# the next function when it is run from within the simulation next time
# check whether the gravity deposition has already finished
# if so, add wall on the top of the packing and start the oedometric test
def checkUnbalanced():
 # at the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable
 if O.iter < 25000:
  return
 # add plate at the position on the top of the packing
 # the maximum finds the z-coordinate of the top of the topmost particle
 O.bodies.append(wall(max([b.state.pos[2] + b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)]), axis=2, sense=-1))
 global plate # without this line, the plate variable would only exist inside this function
 plate = O.bodies[-1] # the last particles is the plate
 # Wall objects are "fixed" by default, i.e. not subject to forces
 # prescribing a velocity will therefore make it move at constant velocity (downwards)
 plate.state.vel = (0, 0, -0.8)
 # start plotting the data now, it was not interesting before
 O.engines = O.engines + [PyRunner(command='addPlotData()', iterPeriod=200)]
 # next time, do not call this function anymore, but the next one (unloadPlate) instead
 checker.command = 'unloadPlate()'

def unloadPlate():
 # if the force on plate exceeds maximum load, start unloading
 if abs(O.forces.f(plate.id)[2]) > 1e3:
  plate.state.vel *= -0.8
  # next time, do not call this function anymore, but the next one (stopUnloading) instead
  checker.command = 'stopUnloading()'

def stopUnloading():
    if abs(O.forces.f(plate.id)[2]) < 1e2:
        # calculate the volume of the packing
        volume_packing = 0
        num_spheres = 0
        for b in O.bodies:
            if isinstance(b.shape, yade.wrapper.Sphere):
                volume_packing += 4/3 * math.pi * b.shape.radius**3
                num_spheres += 1
        # print the number of spheres and volume of packking
        print("Number of spheres:", "{:d}".format(num_spheres))
        print("V Packing:", "{:e}".format(volume_packing))
        # print the height of the plate
        top = abs(plate.state.pos[2] + height/2)
        print("Plate height:", top)
        # calculate the volume of the cylinder
        new_volume_cylinder = math.pi * (diameter/2)**2 * top
        print("V Cylinder:", "{:e}".format(new_volume_cylinder))
        # calculate bulk volume
        new_bulk_volume = (new_volume_cylinder - volume_packing)
        print("Bulk Volume:", "{:e}".format(new_bulk_volume))
        # calculate bulk density
        new_bulk_density = (density * volume_packing) / new_bulk_volume
        print("Bulk Density:", "{:.2f}".format(new_bulk_density))
        # calculate VMA
        new_vma = (new_bulk_volume - volume_packing) * 100
        print("VMA:", "{:.2f}%".format(new_vma))
        # calculate % air void
        new_air_void = ((new_bulk_volume -volume_packing)/new_bulk_volume)*100
        print("Air Void:", "{:.2f}%".format(new_air_void))
        # calculate the porosity and porosity percentage
        new_porosity = (new_volume_cylinder - volume_packing) / new_volume_cylinder
        new_porosity_percent = new_porosity * 100
        print("Porosity:", "{:.2f}".format(new_porosity))
        print("Porosity:", "{:.2f}%".format(new_porosity_percent))
        # O.tags can be used to retrieve unique identifiers of the simulation
        # if running in batch, subsequent simulation would overwrite each other's output files otherwise
        # d (or description) is simulation description (composed of parameter values)
        # while the id is composed of time and process number
        plot.saveDataTxt(O.tags['d.id'] + '.txt')
        O.pause()

def addPlotData():
 if not isinstance(O.bodies[-1].shape, Wall):
  plot.addData()
  return
 Fz = O.forces.f(plate.id)[2]
 plot.addData(Fz=Fz, w=plate.state.pos[2] - plate.state.refPos[2], unbalanced=unbalancedForce(), i=O.iter)

# besides unbalanced force evolution, also plot the displacement-force diagram
plot.plots = {'i': ('unbalanced',), 'w': ('Fz',)}
plot.plot()