M&M schatten

# -*- coding: utf-8 -*-
"""
Created on Thu Dec 13 14:21:46 2018

@author: mmlho
"""

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt

weightPerTen = np.array([23.7,
21.7,
22.8,
22.5,
25.3,
20.9,
22.5,
22.4,
21.8,
23.4])

weightPerTenSig = 0.1 # meetonzekerheid van de weegschaal was 0.1 gram

avweight = (np.sum(weightPerTen))/(10*len(weightPerTen)) #gemiddelde
weightSig = weightPerTenSig/np.sqrt((10**2*len(weightPerTen)))
"""
Als de onzekerheid op een meting van 10 pinda-M&M's 0.1 is, dan geldt 0.1 = sqrt(10*sigma**2) = sqrt(10)*sigma.
Aangezien elke meting exact hetzelfde aantal pinda's en dezelfde onzekerheid heeft, geldt dit ook voor de sigma van een willekeurige pinda-M&M,
gezien over alle metingen. De SDOM is dan gelijk aan sigma/sqrt(100), want er zijn 100 pinda's gewogen in totaal. Dus SDOM = 0.1/sqrt(1000).
"""
print('Gemiddeld gewicht van een pinda-M&M is {:.3f} +/- {:.3f} gram'.format(avweight, weightSig))

height = 9
diameterTop = 24.8
diameterBottom = 18.4
rim = 0.7
lenSig = 0.2 # alle bovenstaande meetwaarden zijn met meetlint gemeten, en hebben vanwege de nauwkeurigheid van het lint en vanwege het aflezen een onzekerheid van 2 mm
radiusTop = diameterTop/2
radiusBottom = diameterBottom/2

"""
Aangezien het verschil in de straal boven en onder (24.8-18.4)/2 = 3.4 cm is, en de hoogte 9 cm is, zou de zijde van de kom, als het om een perfect vlakke, schuine zijde gaat,
sqrt(3.4**2 + 9**2) = 9.62 cm zijn. Er is 9.7 cm gemeten. Dat verschilt niet al te veel. We nemen dus gewoon aan dat de zijde niet gekromd is. De rand is 0.7 cm dik.
Het volume van de kom rekenen we uit als verschil tussen twee geroteerde trapeziums. Die rekenen we uit met integralen. We nemen de buitenwand, de bovenkant
en de onderkant van de kom als grenzen voor de ene trapezium, en de binnenwand, de bovenkant en de bodem aan de binnenkant van de kom als grenzen. Het verschil daartussen
is precies de inhoud van de kom, in de letterlijke zin: de ruimte die gevuld is met materiaal (niet het volume wat opgevuld kan worden).
"""

x = sp.Symbol('x')
outerEdge = radiusTop - x*(radiusTop - radiusBottom)/height
innerEdge = radiusTop - rim - x*(radiusTop - radiusBottom)/height
volumefOuter = sp.integrate(outerEdge**2, x)
volumefInner = sp.integrate(innerEdge**2, x)
volume = np.pi*(volumefOuter.subs(x, height) - volumefOuter.subs(x, 0) - (volumefInner.subs(x, height - rim) - volumefInner.subs(x, 0)))
volumesig = 'TBA' # neem dit iets groter dan de propagatie zoals hij nu zou zijn, ivm het feit dat de wand niet echt gekromd is
print('Het volume van de rand van de bak is {:.1f} cm^3'.format(volume))

weightFilled = 3851.3
weightFilledSig = 0.1 # zoals gemeten door Peter van Capel

density = 2.4
denisitySig = 0.1 # voor de dichtheid van porcelein was te vinden dat het ongeveer 2.3-2.5 g/cm^3 was, dus 2.4 +/- 0.1

weightEmpty = density*volume
weightEmptySig = 'TBA'

peanutNumber = (weightFilled - weightEmpty)/avweight # afgeschat aantal M&M's
print('Er zitten ongeveer {:d} pinda-M&Ms in de bak'.format(int(peanutNumber)))