RA ˇ CUNALNIˇ SKO SIMULIRANJE SIPANJA SVETLOBE V ATMOSFERI

(1)

UNIVERZA V LJUBLJANI PEDAGOˇ SKA FAKULTETA

EVA MARKELJ

RA ˇ CUNALNIˇ SKO SIMULIRANJE SIPANJA SVETLOBE V ATMOSFERI

DIPLOMSKO DELO

LJUBLJANA, 2016

(2)

UNIVERZA V LJUBLJANI PEDAGOˇ SKA FAKULTETA

DVOPREDMETNI U ˇCITELJ: FIZIKA - MATEMATIKA

EVA MARKELJ

MENTOR: izr. prof. dr. BOJAN GOLLI

RA ˇ CUNALNIˇ SKO SIMULIRANJE SIPANJA SVETLOBE V ATMOSFERI

DIPLOMSKO DELO

LJUBLJANA, 2016

(3)

(4)

Zahvala

Iskreno se zahvaljujem svojemu mentorju, izr. prof. dr. Bojanu Golliju, najprej za zanimiva predavanja, zaradi katerih sem se zaˇcela ˇse bolj zanimati za fiziko. Hvala tudi za kvalitetne nasvete, usmeritve in stokovno pomoˇc pri nastajanju diplomskega dela. Cenim vso motivacijo in ˇcas, ki mi ga je profesor namenil.

Rada bi se zahvalila svoji druˇzini, ˇse posebej starˇsema Milanu in Metki, za finanˇcno in moralno podporo.

Hvala tudi soˇsolcem in prijateljem za popestritev ˇstudentskega ˇzivljenja. Posebej bi se rada zahvalila prijatelju Timu za pomoˇc pri programiranju.

Hvala!

(5)

(6)

Povzetek

Diplomsko delo zajema raˇcunalniˇsko simuliranje, v katero je vkljuˇcena metoda Monte Carlo. Metodo najprej predstavim, na kratko povzamem njeno zgodovino in opiˇsem njene sestavne dele. V nadaljevanju opiˇsem nekaj laˇzjih primerov raˇcunalniˇskih simulacij in predstavim njihove rezultate. Ti primeri so: raˇcunanje pribliˇzka ˇstevila π, nakljuˇcna hoja v eni dimenziji, nakljuˇcna hoja v dveh dimenzijah in nevtronski reflektor. Problemi si sledijo po zahtevnosti od najpreprostejˇsega do najzahtevnejˇsega. Nazadnje predstavim ˇse glavni problem dela, to je sipanje svetlobe v atmosferi. Opiˇsem tudi raˇcunalniˇsko simulacijo in njene rezultate. Preko izraˇcunov in grafov pojasnim, zakaj je nebo modro in zakaj je sonce ob sonˇcnem zahodu rdeˇce. Vkljuˇcila sem tudi sliki neba in zahajajoˇcega sonca, ki sta obarvani glede na deleˇze posamezne osnovne barve, pridobljene z raˇcunalniˇsko simulacijo.

Kljuˇcne besede: metoda Monte Carlo, nakljuˇcna ˇstevila, raˇcunalniˇska simulacija, sipanje svetlobe v atmosferi, modro nebo, rdeˇce sonce

(7)

(8)

Abstract

In my diploma thesis I present computer simulation based on the Monte Carlo method.

In the beginning I describe the method and some history behind it. I also present its com- ponents, which are important for using the method. In the following chapter, I present some simple examples of computer simulations. They are organized by difficulty, from basic to the more complicated ones. These examples are: calculating the mathematical constantπ, random walk in one dimension, random walk in two dimensions, and neutron reflector. After that I introduce the main problem of the thesis, light scattering in the atmosphere. I describe the problem and its implementation in the computer simulation in detail. Using the results of my calculations I explain why the sky is blue when the sun is in the zenith and also why the sun is red at the sunset. I also made a picture of the sky and of the sun at the sunset, which are coloured according to the results of simulation.

Key words: Monte Carlo method, random numbers, computer simulation, scattering light in the atmosphere, blue sky, red sun

(9)

(10)

Kazalo

1 Uvod 1

2 Metoda Monte Carlo 3

2.1 Opis metode . . . 3

2.2 Zgodovina . . . 3

2.3 Sestavni deli . . . 4

2.3.1 Verjetnostna porazdelitev . . . 5

2.3.2 Generator nakljuˇcnih ˇstevil . . . 6

2.3.3 Vzorˇcenje . . . 6

2.4 Integracija in napaka . . . 8

3 Preprosti primeri uporabe metode Monte Carlo 11 3.1 Raˇcunanje pribliˇzka ˇstevila π . . . 11

3.1.1 Opis problema . . . 11

3.1.2 Raˇcunalniˇska simulacija . . . 12

3.2 Nakljuˇcna hoja v eni dimenziji . . . 14

3.3 Nakljuˇcna hoja v dveh dimenzijah . . . 18

3.4 Nevtronski reflektor . . . 21

4 Sipanje svetlobe v atmosferi 29 4.1 Opis problema . . . 29

4.2 Rayleighovo sipanje . . . 29

4.2.1 Sipalni presek in prosta pot . . . 30

(11)

KAZALO

4.2.2 Kotna porazdelitvena funkcija . . . 37

4.2.3 Transformacija vektorja premika v prvotni koordinatni sistem . . . 39

4.3 Raˇcunalniˇska simulacija . . . 41

4.3.1 Rezultati raˇcunalniˇske simulacije . . . 42

4.3.2 Odbojnost in prodornost . . . 42

4.3.3 ˇStevilo sipanj . . . 42

4.3.4 Porazdelitve fotonov po kotu v doloˇceni smeri . . . 44

5 Zakljuˇcek 47

A Priloge

A.1 Koda za raˇcunalniˇsko simulacijo sipanja fotonov . . . .

(12)

Slike

2.1 Inverzna transformacijska metoda . . . 7

2.2 Numeriˇcno integriranje . . . 8

3.1 Nakljuˇcne toˇcke znotraj enotskega kvadrata . . . 12

3.2 Grafσ⁻¹(√ N) . . . 13

3.3 Graf porazdelitve delcev po konˇcnih oddaljenostih od izhodiˇsˇca . . . 17

3.4 Grafhri(N) . . . 19

3.5 Grafhr²i(N) . . . 20

3.6 Graf porazdelitve delcev po konˇcnih oddaljenostih od izhodiˇsˇca . . . 20

3.7 Poenostavljen model sipanja nevtrona na jedru . . . 22

3.8 Porazdelitev nevtronov po ˇstevilu sipanj - model 1 . . . 25

3.9 Porazdelitev nevtronov po ˇstevilu sipanj - model 2 . . . 25

3.10 Graf odbojnost(N) . . . 26

3.11 Graf odbojnost(^a_d) . . . 26

4.1 Opis sipanja s sfernim koordinatnim sistemom . . . 31

4.2 Sipalni presek . . . 32

4.3 Graf deleˇza direktne svetlobe v odvisnosti od valovne dolˇzine pri d= 8 km. 35 4.4 Graf deleˇza direktne svetlobe v odvisnosti od valovne dolˇzine prid⁰ = 24 km. 35 4.5 Graf deleˇza direktne svetlobe v odvisnosti od valovne dolˇzine prid⁰⁰ = 320 km. 36 4.6 Graf deleˇza direktne svetlobe v odvisnosti od kota α. . . 36

4.7 Graf deleˇza direktne svetlobe v odvisnosti od kota αv razmerju z deleˇzem direktne rdeˇce svetlobe . . . 37

4.8 Primer metode zavrnitve, kjer toˇcko T₁ zavrnemo,T₂ pa sprejmemo . . . . 39

4.9 Graf porazdelitve fotonov po ˇstevilu sipanj pri debelini atmosfere 8 km . . 43

4.10 Graf porazdelitve fotonov po ˇstevilu sipanj pri debelini atmosfere 24 km . . 43

4.11 Graf porazdelitve fotonov po sipalnem kotu pri debelini atmosfere 8 km . . 44 4.12 Graf porazdelitve fotonov po sipalnem kotu pri debelini atmosfere 24 km . 45 4.13 Graf porazdelitve fotonov po sipalnem kotu v doloˇceni smeri - prvo sipanje 45 4.14 Graf porazdelitve fotonov po sipalnem kotu v doloˇceni smeri - drugo sipanje 46

(13)

SLIKE

4.15 Graf porazdelitve fotonov po sipalnem kotu v doloˇceni smeri - tretje sipanje 46 4.16 Graf porazdelitve fotonov po sipalnem kotu v doloˇceni smeri - ˇcetrto sipanje 46 5.1 Modro nebo . . . 48

(14)

Tabele

3.1 Primer rezultatov simulacije raˇcunanja pribliˇzka ˇstevila π . . . 12

3.2 Primer rezultatov simulacije nakljuˇcne hoje v 1D za milijon delcev . . . 17

3.3 Primer rezultatov simulacije nakljuˇcne hoje v 2D za milijon delcev . . . 19

4.1 Izraˇcuni sipalnega preseka in proste poti . . . 34

4.2 Prilagajanje intervala ˇzrebanja prve opravljene poti za d = 8 km, d⁰ = 240 km ind⁰⁰ = 320 km . . . 34

4.3 Odbojnost in prodornost fotonov pri debelinahd= 8 km in d⁰ = 24 km . . 42

(15)

(16)

Poglavje 1 Uvod

Vsem se zdi samoumevno, da je nebo modre barve. Redko kateri pa zna odgovoriti na vpraˇsanje, zakaj je temu tako. Glavni problem diplomske naloge je sipanje fotonov v Zemljini atmosferi. Ker so premeri molekul zraka manjˇsi v primerjavi z valovno dolˇzino svetlobe, prevladuje pojav, ki ga imenujemo Rayleighovo sipanje. To je elastiˇcno sipanje, pri katerem se svetloba s krajˇso valovno dolˇzino bolj razprˇsi od tiste, ki ima daljˇso. Tako lahko pojasnimo, zakaj je nebo podnevi modro in zakaj sonce rumeno. Modra svetloba se siplje pribliˇzno desetkrat bolj kot rdeˇca, zato torej nebo v smereh proˇc od sonca vidimo modre barve. Ker pa so v sonˇcni svetlobi zastopane vse valovne dolˇzine vidnega spektra, v smeri sonca prihaja svetloba rumene barve. Z istim pojavom lahko razloˇzimo tudi rdeˇce zahajajoˇce sonce. Takrat je pot, ki jo prepotujejo fotoni, daljˇsa, zato se modra svetloba ˇse bolj razprˇsi in v sonˇcnih ˇzarkih ostane rdeˇckasta svetloba.

Izdelala sem raˇcunalniˇski program, pri katerem sem simulirala potovanje posameznega fotona skozi atmosfero. Pri tem sem upoˇstevala, da imajo fotoni z modro valovno dolˇzino na molekulah zraka veˇcji sipalni presek kot fotoni z rdeˇco. Sreˇcanja fotonov z molekulami atmosfere so nakljuˇcna, zato sem v program vkljuˇcila numeriˇcno metodo Monte Carlo, ki za svoje izvajanje uporablja zaporedje nakljuˇcnih ˇstevil.

Metoda Monte Carlo se uporablja v znanosti za preuˇcevanje sistemov, pri katerih ne moremo priti do analitiˇcne reˇsitve. Namenjeno ji je prvo poglavje diplomske naloge. Me- toda je najprej opisana, nato pa je na kratko povzeta tudi njena zgodovina. Navedene in opisane pa so tudi njene glavne sestavine, ki jih je potrebno poznati, ˇce metodo ˇzelimo uporabljati.

Drugo veˇcje poglavje je namenjeno preprostim primerom uporabe metode Monte Carlo. Vsak problem je najprej opisan skupaj s pripadajoˇco raˇcunalniˇsko simulacijo.

1

(17)

2 POGLAVJE 1. UVOD

Predstavljeni so tudi rezultati v obliki tabel in grafov. Prvi preprosti primer je raˇcunanje pribliˇzka matematiˇcne konstanteπ. Ta zelo nazorno predstavi metodo in je zelo zanimiv, saj pokaˇze, da se preko nakljuˇcnih ˇstevil da zelo natanˇcno doloˇciti ˇsteviloπ. Kot drug in tretji primer je obravnavana nakljuˇcna hoja v eni in dveh dimenzijah. Nakljuˇcna hoja je preprost primer, toda osnova za ˇstevilne fizikalne pojave, kot je na primer Brownovo gibanje. Nakljuˇcna ˇstevila so tu uporabljena za izbiro smeri gibanja. ˇCetrti primer pa je nevtronski reflektor. Gre za sipanje nevtronov v snovi, pri ˇcemer so nakljuˇcne njihove proste poti in kot, pod katerem se sipljejo. Ta primer je zelo zanimiv, saj gre za realistiˇcen primer in tudi nekaj, s ˇcimer se ukvarjajo jedrski inˇzenirji.

Zadnje poglavje govori o glavnem problemu, o sipanju svetlobe. Opisana je tudi raˇcunalniˇska simulacija in predstavljena porazdelitev fotonov doloˇcene valovne dolˇzine po sipalnem kotu. Meritve so opravljene za atmosfero debeline 8 kilometrov in ˇse za trikrat veˇcjo. Dalje je zapisana ˇse interpretacija rezultatov in odgovori na raziskovalna vpraˇsanja diplomske naloge. Na koncu diplomske naloge je priloˇzena tudi raˇcunalniˇska koda simulacije sipanja fotonov v atmosferi, napisana v programskem jeziku C.

(18)

Poglavje 2

Metoda Monte Carlo

2.1 Opis metode

Metoda Monte Carlo je vsaka numeriˇcna metoda, ki za svoje izvajanje uporablja zaporedje nakljuˇcnih ˇstevil. Uporablja se pri reˇsevanju problemov, zlasti tistih, pri katerih klasiˇcni pristop odpove. Nekaterih pojavov ne moremo napovedati deterministiˇcno. Pogosto je znana le verjetnost za posamezen dogodek. [1]

Pri metodi izvajamo nakljuˇcne poskuse na mikroskopskem nivoju, pri ˇcemer naredimo veliko ˇstevilo meritev. Te meritve potem statistiˇcno obdelamo in dobimo makroskop- sko reˇsitev problema. Z njeno uporabo simuliramo fizikalne pojave ali matematiˇcne prostore. [2]

2.2 Zgodovina

Ideja o metodi je veliko starejˇsa od iznajdbe raˇcunalnika. Ime je dobila okoli leta 1940 po slavnem kazinoju v monaˇski kneˇzevini zaradi iger na sreˇco in njihovo povezavo z na- kljuˇcnimi dogodki. Poimenoval jo je Nicholas Metropolis, ki je tudi veliko prispeval k njenemu razvoju. Pred tem so jo poznali pod imenom statistiˇcno vzorˇcenje. [3]

Zelo slaven primer uporabe metode brez raˇcunalnika je poskus Buffonova igla iz leta 1733, preko katerega se da doloˇciti matematiˇcno konstanto π. Pribliˇzek je doloˇcen z metanjem igle dolˇzine L na papir, na katerem so narisane vzporedne daljice na razdalji D, pri ˇcemer velja: L ≤ D. Verjetnost, da igla pade in prekriˇza vsaj eno daljico, je P = _πD^2L. Pri poskusu preˇstejemo ˇstevilo vseh metov N in ˇstevilo prekriˇzanjR. Empiriˇcno verjetnost izraˇcunamo tako, da ˇstevilo ugodnih dogodkov delimo s ˇstevilom vseh dogodkov.

3

(19)

4 POGLAVJE 2. METODA MONTE CARLO

Po preoblikovanju zvez dobimo enaˇcbo za izraˇcun ˇstevilaπ:

π = 2LN

RD . (2.1)

Zanimivo pri tem eksperimentu je, da se z nakljuˇcnim poskusom da zelo natanˇcno doloˇciti matematiˇcno konstanto. [3]

Matematiki so metodo dolgo ˇcasa uporabljali kot aproksimacijo za reˇsevanje diferencialnih enaˇcb in integralov, za kar v tistem ˇcasu ni bilo drugih metod, ki bi pripeljale do analitiˇcnih reˇsitev. V teh letih se je metoda kar precej razvila in ljudje so nehali dvomiti o rezultatih, ki jih je ponujala. A vendar ni bila popolnoma izkoriˇsˇcena v znanstvenoraz- iskovalne namene, paˇc pa bolj v didaktiˇcne. [2]

Prvi veˇcji in pomembni izraˇcuni s pomoˇcjo metode Monte Carlo so bili izvedeni pri ˇstudiju sipanja in absorpcije nevtronov. S tem sta se ukvarjala predvsem John von Neu- mann in Stanislaw Marcin Ulam. Delala sta na projektu Manhattan pri izdelavi vodikove bombe v Los Alamosu. V tem obdobju so se zaradi pomembnosti in sredstev projekta razvijali modernejˇsi algoritmi.

Ker eksperimentov niso mogli dejansko izvesti, so jih simulirali. Pojavi, s katerimi so se ukvarjali, so nakljuˇcni, zato je bilo primerno, da so v raˇcunalniˇsko simulacijo vkljuˇcili nakljuˇcna ˇstevila. Ugotovili so, da lahko z uporabo matematiˇcnega modela pridejo do re- alistiˇcnih statistiˇcnih rezultatov. Potrebovali so le dober generator primernih nakljuˇcnih ˇstevil. [3]

Metoda se je po uspeˇsnih uporabah priˇcela ˇsiriti. Uporabili so jo celo pri reˇsevanju Schr¨odingerjeve enaˇcbe in nelinearnih paraboliˇcnih parcialnih diferencialnih enaˇcb ter celo pri razvoju raˇcunalnikov kot hevristiˇcno iskanje algoritmov. [4]

Skupaj z razvojem raˇcunalnikov in raˇcunalniˇskim programiranjem se je razvijala tudi metoda, kar je privedlo do veˇcje preciznosti in prepriˇcljivosti. Danes jo uporabljajo na razliˇcnih podroˇcjih: v fiziki, astronomiji, meteorologiji, matematiki, kemiji, biologiji, inˇzenirstvu, tehnologiji, ekonomiji, raˇcunalniˇstvu, prometu . . .

2.3 Sestavni deli

Glavne komponente metode Monte Carlo so verjetnostna porazdelitev, generator na- kljuˇcnih ˇstevil in vzorˇcenje. Zahtevnejˇse simulacije pa sestavljajo ˇse beleˇzenje izidov,

(20)

2.3. SESTAVNI DELI 5

ocena napake, metoda zmanjˇsanja odklona ter algoritmi, ki omogoˇcajo, da se metoda Monte Carlo lahko izvrˇsuje na zmogljivejˇsih raˇcunalnikih. [2]

2.3.1 Verjetnostna porazdelitev

Nakljuˇcna spremenljivka je spremenljivka, ki lahko zavzame veˇc kot eno vrednost, ki je ni moˇc predvideti vnaprej. Pogosto poznamo le njeno porazdelitev. Ta nam pove, kolikˇsna je verjetnostP, da spremeljivka zavzame doloˇceno vrednost. Za simuliranje pojava z metodo Monte Carlo moramo poznati gostoto verjetnosti f(x). Izraz f(x)dx nam predstavlja verjetnost, da se nakljuˇcna vrednost x⁰ nahaja znotraj intervala [x, x+ dx]: [1]

f(x)dx=P(x⁰ ∈[x, x+ dx]). (2.2) Verjetnost ima vrednost velikosti ploˇsˇcine pod krivuljo, kar zapiˇsemo kot doloˇceni integral. ˇCe nas zanima, kolikˇsna je verjetnost, da se sluˇcajna spremenljivka nahaja na poljubnem intervalu [c, d], jo izraˇcunamo kot:

P(x⁰ ∈[c, d]) = Z d

c

f(x)dx . (2.3)

Za gostoto verjetnosti velja:

• Vedno je pozitivna ali enaka 0:

f(x)≥0. (2.4)

• Je normalizirana:

∞

Z

−∞

f(x)dx= 1. (2.5)

Definirajmo ˇse porazdelitveno funkcijo c(x):

c(x⁰) = Z x⁰

xmin

f(x)dx , (2.6)

ki je monotono naraˇsˇcujoˇca funkcija, definirana na intervalu [0,1].

Ce zgornjo definicijo obrnemo, dobimo zvezo:ˇ f(x) = dc(x)

dx . (2.7)

(21)

2.3.2 Generator nakljuˇ cnih ˇ stevil

Simulacija Monte Carlo za svoje delovanje potrebuje programsko opremo, ki generira ˇcim bolj nakljuˇcna ˇstevila. Zaporedje pravih nakljuˇcnih ˇstevil je nepredvidljivo in zato tudi teˇzko ustvarljivo. Ustvarimo ga lahko preko nakljuˇcnih fizikalnih pojavov, kot sta na primer radioaktivni razpad in termalni hrup v elektronskih napravah. Izkaˇze pa se, da je uporaba takˇsnega naˇcina generiranja nakljuˇcnih ˇstevil nepraktiˇcna, saj zaporedje ni ponovljivo.

Hitrejˇsa in enostavnejˇsa je uporaba psevdo-nakljuˇcnih ˇstevil. To je zaporedje ˇstevil, ki statistiˇcno gledano ni nakljuˇcno, ima pa podobne lastnosti. ˇCleni zaporedja so doloˇceni z matematiˇcno rekurzivno zvezo. Najbolj znan je linearni kongruenˇcni algoritem z rekurzivno zvezo:

x_i+1 = (ax_i+b) modm , (2.8)

kjer so konstantiainb, modul mter zaˇcetni ˇclenx0 pozitivna cela ˇstevila, pri ˇcemer velja ˇse: m > x₀ in m > a. Dobimo ˇstevila, ki so med 0 in m−1 enakomerno porazdeljena.

Pri tem algoritmu ˇzelimo, da je cikel ˇcim daljˇsi, saj se v nasprotnem primeru zaporedje ˇstevil med delovanjem simulacije ponovi. Dolˇzina cikla je odvisna od izbire konstant in od maksimalnega ˇstevila, ki ga generator premore. [1]

2.3.3 Vzorˇ cenje

Pri simulaciji Monte Carlo ˇzelimo dobiti vzorce, ki so porazdeljeni po doloˇceni gostoti verjetnosti f(x). Zelo uporabni metodi vzorˇcenja sta inverzna transformacijska metoda in metoda zavrnitve. [6]

Inverzna transformacijska metoda

1. Najprej izberemo nakljuˇcno ˇstevilo ξ⁰ ∈ [0,1]. Stevilaˇ ξ so na tem intervalu porazdeljena enakomerno.

2. ˇStevilo ξ⁰ izenaˇcimo s porazdelitveno funkcijo:

ξ⁰ =c(x⁰) = Z x⁰

xmin

f(x)dx . (2.9)

3. Spremenljivko x⁰ izrazimo tako, da izraˇcunamo inverz porazdelitvene funkcije:

x⁰ =c⁻¹(ξ⁰). (2.10)

(22)

2.3. SESTAVNI DELI 7

Slika 2.1: Inverzna transformacijska metoda

Metoda zavrnitve

V primeru zapletenega sistema se zgodi, da se inverzna funkcija verjetnostne porazdelitve ne da doloˇciti v analitiˇcni obliki. Takrat uporabimo bolj zapletene metode, s pomoˇcjo katerih vzorˇcimo spremenljivko. Ena izmed njih je metoda zavrnitve, ki je zasnovana na preprosti geometriji. [6]

Najprej nariˇsemo graf gostote verjetnosti f(x). Pod krivuljo te funkcije ˇzelimo generirati enakomerno porazdeljene nakljuˇcne toˇcke. To naredimo z naslednjim postopkom:

1. Doloˇcimo primerjalno funkcijog(x), ki ima pod krivuljo nad osjox konˇcno ploˇsˇcino ter za katero velja: g(x)≥f(x) za ∀x∈[x_min, x_max].

Ta funkcija mora biti analitiˇcna, prav tako tudi njen nedoloˇceni integral in inverz.

Najbolj priroˇcno je izbrati kar konstantno funkcijo. Primerjalna funkcija vedno obstaja, saj je naˇsa gostota verjetnosti normalizirana, kar pomeni, da ima ploˇsˇcino pod krivuljo in nad osjo x konˇcno.

2. Sedaj moramo izbrati nakljuˇcno toˇcko, ki leˇzi nekje v ravnini pod krivuljo primerjalne funkcije g(x).

Oznaˇcimo zA ploˇsˇcino pod krivuljo primerjalne funkcije, na obmoˇcju, kjer je x∈[x_min, x_max]. Ploˇsˇcina pod krivuljo se izraˇcuna kot doloˇceni integral:

A=

Z xmax

xmin

g(x)dx. (2.11)

V naslednjem koraku izˇzrebamo A⁰ ∈ [0, A], pri ˇcemer velja: A⁰ = x⁰_iA, kjer so nakljuˇcna ˇstevilax_i enakomerno porazdeljena na intervalu [0,1].

(23)

Velja:

A⁰ = Z x⁰

xmin

g(x)dx . (2.12)

Po metodi transformacije, ki sem jo opisala zgoraj, lahko sedaj doloˇcimo koordinato x⁰. Izˇzrebati moramo ˇse nakljuˇcno koordinatoy⁰, enakomerno porazdeljeno na intervalu [0, g(x⁰)].

3. V primeru, da toˇcka (x⁰, y⁰) leˇzi zunaj obmoˇcja pod porazdelitveno funkcijo f(x), toˇcko zavrnemo in postopek ponovimo. V tem primeru velja: y⁰ > f(x⁰). V nasprotnem primeru pa toˇcko sprejmemo, kar pomeni, da uporabimo izˇzrebano spremen- ljivkox⁰. Na tak naˇcin dobimo sluˇcajne vrednosti spremenljivkex, ki so porazdeljene po ˇzeleni porazdelitvi.

2.4 Integracija in napaka

Izberimo siN nakljuˇcnih toˇck x₁, x₂, ..., x_N, enakomerno porazdeljenih na intervalu [a, b].

Pri metodi Monte Carlo integral nadomestimo z vsoto: [5]

Z b

a

f(x)dx≈(b−a)1 N

N

X

i=1

f(x_i) = (b−a)f . (2.13) Ker so ˇstevila x_i nakljuˇcna, so tudi vrednosti funkcijef(x_i) nakljuˇcne. Integral lahko aproksimiramo z vsoto nakljuˇcnih vrednosti funkcije.

Slika 2.2: Numeriˇcno integriranje

(24)

2.4. INTEGRACIJA IN NAPAKA 9

Napaka integrala je posledica raˇcunanja povpreˇcja funkcije. ˇCe namesto povpreˇcne vrednosti vzamemo nakljuˇcno vrednost funkcije, v povpreˇcju storimo napako σf. Vemo, da se napaka s ˇstevilom poskusov spreminja kot:

σ²_f(N) = 1

Nσ²_f. (2.14)

Napaka integrala je sorazmerna s to napako:

σI = (b−a)σf(N) = (b−a) 1

√Nσf. (2.15)

Napaka metode Monte Carlo se torej z veˇcanjem ˇstevila ponovitev manjˇsa s korensko odvisnostjo:

σ_I ∝ 1

√N . (2.16)

(25)

(26)

Poglavje 3

Preprosti primeri uporabe metode Monte Carlo

3.1 Raˇ cunanje pribliˇ zka ˇ stevila π

3.1.1 Opis problema

V podpoglavju Zgodovina 2.2 smo ˇze omenili naˇcin, kako preko nakljuˇcnih dogodkov doloˇcimo pribliˇzek matematiˇcne konstante π. Doloˇcimo ga lahko tudi na drugaˇcen naˇcin.

V enotski kvadrat vˇcrtamo ˇcetrtino kroga. Ploˇsˇcina enotskega kvadrata je 1. ˇCetrtina kroga s polmerom dolˇzine 1 pa ima ploˇsˇcino ^π₄. Znotraj kvadrata generiramo toˇcke tako, da izˇzrebamo njihove koordinate. V primeru, da je uporabljen generator dober in izˇzrebamo dovolj veliko ˇstevilo toˇck, bodo izbrane toˇcke po kvadratu razporejene pribliˇzno enakomerno. To lastnost vidimo tudi na slikah 3.1.a) in 3.1.b). ˇCe je to res, bo razmerje ploˇsˇcin ˇ

cetrtine kroga in kvadrata enako razmerju toˇck znotraj ˇcetrtine kroga in vseh generiranih toˇck.

Oznaˇcimo celotno ˇstevilo toˇck s ˇcrko N in ˇstevilo toˇck, ki padejo znotraj ˇcetrtine kroga, s ˇcrko M. V limiti, ko gre ˇstevilo toˇck N proti neskonˇcno, velja:

M

N = S¹

4kroga

S_kvadrata = π

4. (3.1)

Steviloˇ π lahko doloˇcimo kot:

π= 4M

N . (3.2)

11

(27)

12 POGLAVJE 3. PREPROSTI PRIMERI UPORABE METODE MONTE CARLO

Za doloˇcanje matematiˇcnega ˇstevila π obstaja veliko uˇcinkovitejˇsih metod. Opisan naˇcin sicer zelo nazorno predstavi metodo Monte Carlo, vendar je dokaj zamuden.

(a) 100 nakljuˇcnih toˇck (b) 500 nakljuˇcnih toˇck

Slika 3.1: Nakljuˇcne toˇcke znotraj enotskega kvadrata

3.1.2 Raˇ cunalniˇ ska simulacija

V raˇcunalniˇski simulaciji sem generirala N nakljuˇcnih toˇck v ravnini znotraj enotskega kvadrata. Program z ukazom rand() generira nakljuˇcno ˇstevilo z algoritmom (2.8). Ta ˇstevila so enakomerno porazdeljena med ˇstevilom 0 in nekim velikim celim ˇstevilom (RAN D M AX). Vsaka toˇcka T(x, y) ima dve koordinati, ki sta nakljuˇcni ˇstevili iz intervala [0,1]. Ta ˇstevila sem dobila s predpisom: j = rand()/(RAN D M AX + 1), kar pomeni, da sem nakljuˇcno ˇstevilo delila s ˇstevilom vseh moˇznih ˇstevil. Ne smemo pozabiti na ˇstevilo 0. V primeru, da je toˇcka padla znotraj ˇcetrtine kroga, je ustrezala pogoju:

x²+y² ≤1.V tem primeru sem jo ˇstela pod mnoˇzico toˇck M. Na koncu sem izraˇcunala ˇse pribliˇzek ˇstevila π po enaˇcbi (3.2) in doloˇcila odstopanje od prave vrednosti. Preverila sem tudi, v kakˇsni zvezi sta napaka in ˇstevilo toˇck.

Tabela 3.1: Primer rezultatov simulacije raˇcunanja pribliˇzka ˇstevila π N-ˇstevilo korakov π˜ σ √

N σ⁻¹

100 3,14223 0,13102 10,00 7,63250

500 3,14193 0,05875 22,36 17,02139

1000 3,14152 0,04144 31,62 24,13374

1500 3,14152 0,03388 38,73 29,51889

2000 3,14137 0,02937 44,72 34,04301

2500 3,14155 0,02606 50,00 38,37902

Iz rezultatov v tabeli 3.1 opazimo, da se natanˇcnost doloˇcanja konstante s poveˇcevanjem ˇstevila toˇck izboljˇsuje.

(28)

3.1. RA ˇCUNANJE PRIBLI ˇZKA ˇSTEVILAπ 13

Slika 3.2: Graf σ⁻¹(√ N)

Iz grafa razberemo, da je napaka res obratno sorazmerna s korenom ˇstevila toˇck.

Rezultatom se najbolje prilega funkcija: σ⁻¹(√

N) = 0,76072√ N .

(29)

3.2 Nakljuˇ cna hoja v eni dimenziji

3.2.1 Opis problema

Z nakljuˇcno hojo lahko simuliramo veliko pojavov, kot so na primer Brownovo gibanje, finanˇcni modeli in igre na sreˇco.

Delec se ob ˇcasu t = 0 nahaja v izhodiˇsˇcu in koraka po premici. Premika se tako, da nakljuˇcno izbere smer gibanja. Zanima nas, kolikˇsna je verjetnost, da bo po N korakih delec na koordinati x.

Na primer, da se delec premakne le enkrat. Oznaˇcimo s P₊₁ verjetnost, da bo ta premik v desno smer ter s P−1 verjetnost, da bo v levo smer. V naˇsem primeru sta ti dve verjetnosti enaki, poleg tega pa velja ˇse: P₊₁+P₋₁ = 1. Poslediˇcno je njuna vrednost enaka: P₊₁ =P−1 = ¹₂. S k oznaˇcimo vrednost koraka. V naˇsem primeru ima vsak korak dolˇzino 1. Velja torej k ∈ {1,−1}. Vrednosti sta enako verjetni, zato je njuna povpreˇcna vrednost enaka 0. Ker je verjetnost, da se delec nahaja nekje na desni strani, enaka kot verjetnost, da se nahaja na levi, se v povpreˇcju nahaja v izhodiˇsˇcu. Velja: hxi= 0.

Oznaˇcimo z x_N mesto, na katerem smo po N korakih, in s k_N korak naslednjega premika. Velja rekurzivna zveza: x_N =x_N−1 +k_N. Povpreˇcna vrednost oddaljenosti od izhodiˇsˇca po N korakih je ˇse vedno 0, prav tako povpreˇcna vrednosti koraka. Poglejmo si sedaj povpreˇcne vrednosti kvadrata oddaljenosti od izhodiˇsˇca.

hx²_ii=h(xi−1+k_i)²i=hx²_i−1+k²_i + 2xi−1k_ii=hx²_i−1i+hk_i²i+ 2hxi−1k_ii (3.3)

hk_i²i= (±1)² = 1 (3.4)

hx_i−1k_ii=hx_i−1(+1)P₊₁i+hx_i−1(−1)P₋₁i= 0 (3.5) Dobljene zveze za i=N vstavimo v enaˇcbo (3.3):

hx²_Ni=hx²_N−1i+ 2·0 + 1 =hx²_N₋₁i+ 1 =hx²_N₋₂i+ 1 + 1 =. . .=hx²₀i+ 1 + 1 +. . .+ 1 =N . (3.6)

Za standardno deviacijo σ² smo dobili izraz:

σ² =hx²i − hxi² =N −0 =N . (3.7)

(30)

3.2. NAKLJU ˇCNA HOJA V ENI DIMENZIJI 15

Centralni limitni izrek pravi, da se vsota enakomerno porazdeljenih sluˇcajnih spremeljivk v limiti, ko greN proti neskonˇcno, pribliˇzuje normalni porazdelitvi. Tej porazdelitvi pravimo tudi Gaussova porazdelitev s predpisom:

f(x, N) = r 1

2πσ²e⁻^x

2

2σ2. (3.8)

Za standardno deviacijo σ² velja:

σ² =

N

X

i=0

σ²_i =N σ₁² =N . (3.9)

Izraz je enak, ko smo ga izpeljali zgoraj.

Naˇso porazdelitev lahko torej aproksimiramo z normalno porazdelitvijo s predpisom:

f(x, σ²)≈2 r 1

2πσ²e⁻^x

2

2σ2 . (3.10)

Normalizacijska konstanta ima vrednost 2, zato ker je ∆x= 2. Interval dolˇzine ∆xima sre- dino v koordinatixin sega dox−1 na eni strani inx+1 na drugi strani. Njegova ˇsirina je torej 2.

Standarno deviacijo σ² nadomestimo s ˇstevilom korakovN in dobimo:

f(x, N)≈2 r 1

2πNe⁻^x

2

2N . (3.11)

Zveza z difuzijsko enaˇcbo

Predpostavimo, da je ˇcas enega koraka konstanten. V tem primeru dobimo, da se standardna deviacija s ˇcasom linearno poveˇcuje: σ² =N =kt.

Vzeli bomo tudi ˇsirino intervala ∆x= 1,da bo porazdelitev normirana. Dobimo porazdelitev s predpisom:

f(x, t) = r 1

2πkte⁻^x

2

2kt . (3.12)

Vstavimo porazdelitev v difuzijsko enaˇcbo in pokaˇzimo, da ji zadoˇsˇca. Difuzijska enaˇcba ali drugi Fickov zakon ima obliko:

D∂²f

∂x² = ∂f

∂t . (3.13)

(31)

Posebej izraˇcunamo odvode:

∂f

∂t =

r 1 2πkte⁻^x

2 2kt

−1 2t+ x²

2kt²

, (3.14)

∂f

∂x =

r 1 2πkte⁻^x

2 2kt

−x kt

, (3.15)

∂²f

∂x² =

r 1 2πkte⁻^x

2 2kt

x² k²t² − 1

kt

. (3.16)

Izraˇcunane odvode vstavimo v difuzijsko enaˇcbo (3.13).

r 1 2πkte⁻^x

2 2kt

−1 2t + x²

2kt²

= D

r 1 2πkte⁻^x

2 2kt

x² k²t² − 1

kt

−1 2t + x²

2kt²

= D

x² k²t² − 1

kt

D = k

2

Ce vstavimo dobljeno zvezo v izraz za standardno deviacijo:ˇ σ² =kt, dobimo:

σ²(t) = 2Dt . (3.17)

3.2.2 Raˇ cunalniˇ ska simulacija

V simulaciji sem nakljuˇcna ˇstevila uporabila pri izbiri smeri premika. Premik delca sem definirala z enaˇcbo: ∆x = 2m−1. ˇStevilo m je moralo torej imeti vrednosti 0 ali 1. To sem dosegla s predpisom:

m= (int) (2·rand()/(RAN D M AX+ 1)) , (3.18) pri ˇcemer ˇzrebamo nakljuˇcna ˇstevila enakomerno porazdeljena po intervalu [0,2] in vzamemo le celi del tega ˇstevila.

Iz rezultatov v tabeli 3.2 razberemo, da res velja hxi ≈0 inhx²i ≈N.

Povpreˇcna vrednost oddaljenostihxi teoretiˇcno ni 0, ampak se njena vrednost razlikuje od 0 za napako, ki je obratno sorazmerna s ˇstevilom delcev. To smo izpeljali pri zvezi (2.16). V naˇsem primeru je ˇstevilo ponovitev poskusa kar ˇstevilo delcev, zato je napaka sorazmerna z 1/√

10⁶. Za ˇstevilo korakov N = 1500 in ˇstevilo delcev M = 100000 sem narisala graf porazdelitve delcev po konˇcnih oddaljenosti od izhodiˇsˇca. Gaussova krivulja na grafu 3.3 ima predpis f(x,1500) = 2

q 1 3000πe⁻ ^x

2

3000 in se naˇsim rezultatom simulacije dobro prilega.

(32)

3.2. NAKLJU ˇCNA HOJA V ENI DIMENZIJI 17

Tabela 3.2: Primer rezultatov simulacije nakljuˇcne hoje v 1D za milijon delcev N-ˇstevilo korakov hxi hx²i

50 0 50

100 0 99

150 0 149

200 0 200

250 0 250

300 0 299

350 0 349

400 0 400

450 0 449

500 0 499

N-ˇstevilo korakov hxi hx²i

550 0 549

600 0 599

650 0 650

700 0 699

750 0 750

800 0 801

850 0 850

900 0 899

950 0 950

1000 0 1000

Slika 3.3: Graf porazdelitve delcev po konˇcnih oddaljenostih od izhodiˇsˇca

(33)

3.3 Nakljuˇ cna hoja v dveh dimenzijah

3.3.1 Opis problema

Nakljuˇcno hojo v eni dimenziji sem nadgradila ˇse s premikom naprej in nazaj. Delec se ob ˇcasu t= 0 nahaja v izhodiˇsˇcu in se premika po ravnini. Pred vsakim korakom nakljuˇcno izbere eno izmed ˇstirih smeri gibanja.

Teoretiˇcno naj bi bila porazdelitev delcev po oddaljenostih od izhodiˇsˇca Gaussova krivulja pomnoˇzena s spremeniljivkor z enaˇcbo:

f(r, σ²) = 2 σ²re⁻^r

2

2σ2 . (3.19)

Kot pri nakljuˇcni hoji v eni dimenziji ima tudi tu standardna deviacija vrednost N, zato lahko zapiˇsemo porazdelitev ˇse v drugaˇcni obliki:

f(r, N) = 2 Nre⁻^r

2

2N. (3.20)

Doloˇcimo ˇse izraz za povpreˇcno vrednost oddaljenosti in njenega kvadrata.

hri= R∞

0 rf(r, N)dr R∞

0 f(r, N)dr = 1 2

Z ∞ 0

2 Nr²e⁻^r

2 2Ndr =

Z ∞ 0

1 Nr²e⁻^r

2 2Ndr=

rπ

2N ≈1,25√

N (3.21)

hr²i= R∞

0 r²f(r, N)dr R∞

0 f(r, N)dr = 1 2

Z ∞ 0

2 Nr³e⁻^r

2 2Ndr=

Z ∞ 0

1 Nr³e⁻^r

2

2Ndr = 2N (3.22)

3.3.2 Raˇ cunalniˇ ska simulacija

Nakljuˇcna ˇstevila sem uporabila pri izbiri smeri premika. Kot pri hoji v eni dimenziji sem tudi tu s predpisom (3.18) generirarala ˇstevili m ino, ki sta imeli vrednosti 0 ali 1. Premik delca v dveh smereh sem definirala z enaˇcbama: ∆x= 2m−1 in ∆y= 2o−1.

Iz tabele 3.3 razberemo, da povpreˇcni odmik od izhodiˇsˇca nima veˇc vrednosti 0 in se njegov kvadrat poveˇcuje premo sorazmerno s ˇstevilom korakov.

Graf hri(N) na sliki 3.4 nam prikazuje, da je odvisnost korenska funkcija.

Velja torej: hri ∝√ N .

Iz grafa na sliki 3.5 pa razberemo, da je odvisnost hr²i(N) linearna funkcija.

Velja: hr²i ≈2N.

Odvisnosti sta konsistentni s tistima, ki smo ju izpeljali pri izrazih (3.21) in (3.22).

(34)

3.3. NAKLJU ˇCNA HOJA V DVEH DIMENZIJAH 19

Tabela 3.3: Primer rezultatov simulacije nakljuˇcne hoje v 2D za milijon delcev N-ˇstevilo delcev hxi hx²i

50 8,5165 100,05878

100 12,1564 199,7602

150 14,9799 300,1939

200 17,3355 399,7259

250 19,4230 499,6424

300 21,3132 599,7751

350 23,0476 699,9543

400 24,6591 799,5240

450 26,1891 899,8097

500 27,6216 1000,1520

N-ˇstevilo delcev hxi hx²i

550 28,9953 1100,6851

600 30,2931 1199,4837

650 31,5644 1301,2826

700 32,7383 1398,5745

750 33,9225 1500,1906

800 35,0259 1598,7578

850 36,1299 1700,3543

900 37,1811 1799,2988

950 38,2258 1900,6834

1000 39,2389 2001,6472

Slika 3.4: Graf hri(N)

(35)

Slika 3.5: Graf hr²i(N)

Zanimala me je ˇse porazdelitev delcev po konˇcnih oddaljenostih od izhodiˇsˇca. Meritve so bile izvedene za 100000 delcev in za sprehode s 1500 koraki. Krivulja na grafu ima predpis:

f(r,1500) = ₁₅₀₀² re⁻ ^r

2

3000. Rezultati se ji dobro prilegajo.

Slika 3.6: Graf porazdelitve delcev po konˇcnih oddaljenostih od izhodiˇsˇca

(36)

3.4. NEVTRONSKI REFLEKTOR 21

3.4 Nevtronski reflektor

3.4.1 Opis problema

Nevtroni so osnovni delci brez elektriˇcnega naboja. Poslediˇcno se lahko v snovi brez teˇzav prebi- jejo do jeder in povzroˇcijo razliˇcne reakcije. Ena izmed njih je elastiˇcen trk, pri ˇcemer se nevtron proˇzno odbije v razliˇcne smeri. Dogodek je enak, kot ˇce bi opazovali trk dveh trdih krogel, na primer bilijardnih krogel.

Z nevtronskim reflektorjem se sreˇcamo v jedrskih reaktorjih. Njegova vloga je, da doloˇcen deleˇz nevtronov odbije nazaj do jedrskega goriva, kjer sproˇzijo nadaljne reakcije. Tam je veliko ˇstevilo nevtronov (10⁸/cm³) in veliko ˇstevilo jeder (10²²/cm³), zato nas pri takˇsnih in podobnih poja- vih zanima le povpreˇcno obnaˇsanje delcev. Poslediˇcno za simulacijo pojava uporabimo metodo Monte Carlo. [7]

Poglejmo si zelo enostaven primer. Na ploˇsˇco z doloˇceno debelino je pravokotno na njene meje usmerjen curek nevtronov. V ploˇsˇci se nevtroni le sipljejo in niˇc absorbirajo. Velja, da je porazdelitev proste poti pri sipanju nevtronov eksponentna. Takˇsni porazdelitvi pravimo tudi porazdelitev Poissonovega toka in je znaˇcilna tudi za ˇcase radioaktivnih razpadov.

Gostota porazdelitve po prostih poteh je funkcija:

f(r) =Ae⁻^r^a, (3.23)

kjer jer dolˇzina vektorja premika,apovpreˇcna prosta pot inAnormalizacijska konstanta, ki jo doloˇcimo iz normalizacijskega pogoja.

Z ∞ 0

f(r)dr = 1 Z ∞

0

Ae⁻â^rdr = 1 A(−a)(e⁻^∞â −e⁻â⁰) = 1

A = 1

a Gostota porazdelitve ima obliko:

f(r) = 1

ae⁻^r^a. (3.24)

Za vzorˇcenje spremenljivke r uporabimo inverzno transformacijsko metodo, opisano v podpoglavju Vzorˇcenje 2.3.3:

1. Izberemo nakljuˇcno ˇstevilo ξ_i iz intervala [0,1].

2. Nakljuˇcno ˇsteviloξ_i izenaˇcimo s porazdelitveno funkcijo:

ξi =c(ri) = Z ri

0

f(r)dr . (3.25)

(37)

ˇSteviloξi nam predstavlja ploˇsˇcino pod porazdelitveno funkcijoc(r) na intervalu [0, ri] : ξi=c(ri) =

Z ri

0

f(r)dr = Z ri

0

1

ae⁻^r^adr = 1

a(−a)(e⁻^riâ −e⁻⁰â) = 1−e⁻^riâ . (3.26) 3. Izrazimo iskano spremenljivko ri=c⁻¹(ξi):

ri = −aln(1−ξi) =−aln(ξ_i⁰). (3.27) Nakljuˇcno ˇstevilo ξ_i je enakomerno porazdeljeno na intervalu (0,1]. Enako velja tudi za ˇsteviloξ⁰_i= 1−ξ_i.

Nevtron se po trku premakne v smeri sipalnega kota θ glede na prvotno smer curka za dolˇzinor_i. Nas zanima le projekcija na os, ki je vzporedna vpadnemu curku. Premik v smeri te osi je:

x_i=r_icosθ . (3.28)

Za gostoto porazdelitve sipalnega kota bomo vzeli:

f(θ) = cosθ . (3.29)

Pokazati moramo, da je ta porazdelitev po prostorskem kotu enakomerna. Taki porazdelitvi reˇcemo tudi, da je izotropna.

Slika 3.7: Poenostavljen model sipanja nevtrona na jedru

Gostota porazdelitve nevtronov po preˇcnem preseku jedra, ki ima obliko krogle s polmerom R, je konstantna, kar zapiˇsemo kot:

f(S) = dp

dS =C . (3.30)

(38)

Konstanto C izraˇcunamo iz normalizacijskega pogoja.

Z Z

S

f(S)dS = 1 C

Z Z

S

dS = 1 CS = 1

C = 1

S Gostota porazdelitve ima obliko:

f(S) =C= 1 S = 1

πR² . (3.31)

S slike 3.7 lahko razberemo naslednji zvezi:

θ = π−2α , (3.32)

cosα = r

R. (3.33)

Iz zgornjih zvez izrazimo razdaljor:

r =Rcosα=Rcos

π−θ 2

=Rcos π

2 −θ 2

=Rsinθ

2. (3.34)

Zanima nas gostota porazdelitve po prostorskem kotu Ω, torej f(Ω) = _dΩ^dp,ki jo izraˇcunamo s pomoˇcjo transformacije spremenljivk:

dp dΩ = dp

dS

dS dΩ

= dp dΩ

dS dr

dr dΩ

. (3.35)

Diferencial ploˇsˇcine je:

dS

dr = 2πr= 2πRsinθ

2. (3.36)

Diferencial prostorskega kota je:

dΩ = 2πsinθdθ . (3.37)

Izraˇcunajmo ˇse _dΩ^dr : dr

dΩ = dr

2πsinθdθ = 1 2πsinθ

d(Rsin^θ₂)

dθ = 1

2πsinθ R

2 cosθ

2. (3.38)

(39)

Dobljene zveze vstavimo v enaˇcbo (3.35).

f(Ω) = dp

dΩ = 1

πR²

2πRsinθ 2

1 2πsinθ

R 2 cosθ

2

=

sin^θ₂cos^θ₂ 2πsinθ

=

1 2sinθ 2πsinθ

=

= 1

4π

Porazdelitev po prostorskem kotu je enakomerna funkcija. Poglejmo si ˇse, ali to velja za porazdelitev po kosinusu sipalnega kota.

f(cosθ) = dp

d cosθ = dp dΩ

dΩ d cosθ

=

= 1

4π

2πsinθdθ d cosθ

=

= 1

4π

2πd(cosθ) d cosθ

=

= 1

2

Tudi porazdelitev po kosinusu sipalnega kota je enakomerna funkcija. Vrednosti kosinusa sipalnega kota so torej enakomerno porazdeljene po intervalu [−1,1].To bomo dosegli, ˇce bomo nakljuˇcno spremenljivko cosθ vzorˇcili kot:

cosθj = 2ξj−1, (3.39)

kjer je ξj nakljuˇcno ˇstevilo na intervalu [0,1].

3.4.2 Raˇ cunalniˇ ska simulacija

Naredila sem raˇcunalniˇski simulaciji dveh modelov sipanja nevtronov v snovi. Pri prvem modelu so se nevtroni sipali le naprej in nazaj, pri drugem pa izotropno. V simulacijo sem vkljuˇcila na- kljuˇcna ˇstevila za ˇzrebanje dolˇzine premika. Drugi model pa sem nadgradila ˇse z generiranjem nakljuˇcnih vrednosti kosinusa sipalnega kota.

Za povpreˇcno pot a sem najprej doloˇcila vrednost ^d₂, kjer je d debelina ploˇsˇce. Zanimala me je porazdelitev nevtronov po ˇstevilu sipanj. V obeh primerih je bilo v simulacijo vkljuˇcenih milijon nevtronov.

(40)

Slika 3.8: Porazdelitev nevtronov po ˇstevilu sipanj - model 1

Slika 3.9: Porazdelitev nevtronov po ˇstevilu sipanj - model 2

Opazimo, da je porazdelitev pri drugem modelu bolj fina. Nevtroni se v tem primeru tudi veˇckrat sipljejo in so dalj ˇcasa ujeti v ploˇsˇci.

(41)

Pri razmerju ^a_d= 0,5 sem raziskala odvisnost odbojnosti reflektorja od ˇstevila nevtronov.

Slika 3.10: Graf odbojnost(N)

Iz grafa lahko razberemo, da je odbojnost reflektorja pri prvem modelu vedno manjˇsa od odbojnosti pri drugem modelu. Zdi se, da vrednost z veˇcanjem ˇstevila nevtronov konvergira proti vrednosti ³₄ = 0,75.Pri drugem modelu pa se pribliˇzuje vrednosti 0,83.

Zanimala me je tudi odvisnost odbojnosti reflektorja od razmerja ^a_d. V simulacijo je bilo vkljuˇcenih milijon nevtronov.

Slika 3.11: Graf odbojnost(^a_d)

(42)

Pri obeh modelih je odbojnost blizu 1, ko je razmerje ^a_d 1.Proste poti nevtronov so v tem obmoˇcju zelo majhne v primerjavi z debelino ploˇsˇce. Pri prvem modelu odbojnost z naraˇsˇcanjem vrednosti razmerja ^a_d pada hitreje kot pri drugem. ˇCe bi poveˇcevali prosto pot, da bi bila mnogo veˇcja od debeline, bi bila odbojnost reflektorja enaka niˇc, saj bi veˇcina nevtronov ˇze v prvem koraku izstopila iz ploˇsˇce.

(43)

(44)

Poglavje 4

Sipanje svetlobe v atmosferi

4.1 Opis problema

Svetloba je elektromagnetno valovanje, ki potuje skozi snov, pri ˇcemer interagira z atomi in molekulami. Svetlobo sestavljajo majhni energijski paketi oziroma svetlobni kvanti, ki jih imenujemo fotoni. So brez elektriˇcnega naboja, zato lahko prodrejo globoko v snov, pri ˇcemer se razprˇsijo oziroma sipljejo. [8]

4.2 Rayleighovo sipanje

Rayleighovo sipanje je pojav, poimenovan po angleˇskemu fiziku lordu Rayleighu. Teorijo je na- pisal leta 1871. Gre za elastiˇcno sipanje svetlobe ali drugega elektromagnetnega valovanja na delcih, ki so veliko manjˇsi od valovne dolˇzine vpadnega valovanja. Pribliˇzno razmerje, ki doloˇca, za katero vrsto sipanja gre, je ^2πr_λ ,kjer jer polmer delca inλvalovna dolˇzina vpadne svetlobe.

Rayleighovo sipanje nastopi, ko je razmerje mnogo manjˇse od 1, kar pomeni, da je delec manjˇsi od desetine valovne dolˇzine svetlobe. [9]

Pojav opazimo veˇcinoma v plinih, lahko pa tudi v trdnih in tekoˇcih snoveh. Delec svetlobo najprej absorbira, nato pa nihajoˇce elektriˇcno polje valovanja povzroˇci, da se delec vzbudi in zaˇcne nihati z enako frekvenco. Na ta naˇcin postane majhen sevajoˇci dipol z dipolnim momen- tom: [8]

p(t) =α₀E₀e^iωt, (4.1)

kjer je E =E₀e^iωt elektriˇcno polje valovanja inα polarizabilnost delca. Vsak vzbujen delec v razliˇcnih smereh oddaja elektromagnetno valovanje, ki ga imenujemo sipana svetloba. Pojav lahko poenostavimo in interpretiramo kot elastiˇcno sipanje fotona na molekuli.

29

(45)

30 POGLAVJE 4. SIPANJE SVETLOBE V ATMOSFERI

Gostota toka sipanih fotonov ali intenziteta sipane svetlobe na molekuli s polarizabilnostjo α na razdaljiR je: [9]

I =I0α² 2π

λ 4

1 + cos²θ

2R² , (4.2)

kjer je λvalovna dolˇzina vpadne svetlobe in I₀ njena intenziteta.

Iz enaˇcbe razberemo, da je zveza med intenziteto sipane svetlobe in valovno dolˇzino vpadne svetlobe: I ∝λ⁻⁴. To nam pove, da se svetloba krajˇsih valovnih dolˇzin na delcih atmosfere siplje moˇcneje kot svetloba daljˇsih valovnih dolˇzin. Izraˇcun pokaˇze, da je intenziteta modre sipane svetlobe desetkrat veˇcja od intenzitete rdeˇce sipane svetlobe.

Razprˇsena svetloba daje nebu svetlost in barvo. Svetloba krajˇsih valovnih dolˇzin se siplje intenzivnejˇse, zato je nebo videti modro. V smeri sonca pa ostane svetloba daljˇsih valovnih dolˇzin, zato je sonce videti rumeno oziroma oranˇzno. Ob sonˇcnih vzhodih in zahodih pa je pot, ki jo opravijo fotoni v atmosferi, daljˇsa, zato pride do veˇcjega ˇstevila sipanj in se razprˇsi tudi nekaj svetlobe daljˇsih valovnih dolˇzin. V tem primeru pa vidimo sonce rdeˇce.

4.2.1 Sipalni presek in prosta pot

Kam se je foton po odboju sipal, najbolje opiˇse sferni koordinatni sistem na sliki 4.1, z osjo z usmerjeno v smeri vpadne svetlobe. Vektor premika ima koordinate: −→r = −→r(r, θ, ϕ), kjer je r dolˇzina vektorja, θ polarni kot in ϕ azimutni kot. Polarni kot lahko zavzame vrednosti na intervalu [0, π],azimutni pa na [0,2π].

Foton se po sipanju premakne za vektor:

−

→r = (rcosϕsinθ, rsinϕsinθ, rcosθ). (4.3) Interakcijski geometrijski presek tarˇce, ki prestreˇze tok fotonov, imenujemo sipalni presek.

To ni dejanska meritev fizikalnih lastnosti tarˇce.

Definicija sipalnega preseka σ je:

I0σ = Z

IdS = Z

IR²dΩ. (4.4)

R je oddaljenost delca, na katerem se svetloba siplje, od opazovalca in tudi polmer sfere, po kateri integriramo. Diferencial prostorskega kota je: dΩ = sinθdθdϕ.

(46)

4.2. RAYLEIGHOVO SIPANJE 31

Slika 4.1: Opis sipanja s sfernim koordinatnim sistemom

Doloˇcimo sedaj sipalni presek.

σ = Z I

I0

R²dΩ = (4.5)

=

Z I0α²(^2π_λ )^{4 (1+cos}_2R2²^θ)

I₀ R²dΩ = (4.6)

= α² 2π

λ 4

1 2

Z 2π 0

dϕ Z π

0

(1 + cos²θ) sinθdθ . (4.7) Posebej izraˇcunajmo integrala po kotihϕinθ.

Z 2π 0

dϕ= 2π (4.8)

Z _π

0

(1 + cos²θ) sinθdθ = Z _π

0

(sinθ+ cos²θsinθ)dθ= (4.9)

= Z π

0

sinθdθ+ Z π

0

cos²θsinθdθ= 2 +2 3 = 8

3 (4.10)

Vrednost zgornjih dveh integralov vstavimo v enaˇcbo za sipalni presek (4.7) : σ=α²

2π λ

4

1 22π8

3 = 128π⁵α²

3λ⁴ . (4.11)

(47)

Slika 4.2: Sipalni presek

Dober pribliˇzek zgornjega izraza za sipalni presek je: [10]

σ≈ τ0

λ⁴, (4.12)

kjer je τ₀≈4,4·10⁻¹⁶cm²nm⁴ za zrak.

Proste poti fotonov v snovi so kot pri nevtronih porazdeljene eksponentno. Zato jih vzorˇcimo kot:

r_i =−aln(ξ_i), (4.13)

kjer je ξi nakljuˇcno ˇstevilo, enakomerno porazdeljeno po intervalu (0,1].

Povpreˇcna prosta potaje povpreˇcna pot med dvema interakcijama fotona z delci. Zamislimo si, da vsako molekulo zapremo v kvader s stranicoain presekomσ. Njegova prostornina je torej:

V1=aσ. Vseh molekul je N, zato celotni prostor sestavlja N takˇsnih kvadrov. Velja:

V1=aσ= V N = 1

n, (4.14)

kjer smo zn oznaˇcili ˇstevilsko gostoto molekul zraka.

Povpreˇcno prosto pot aizrazimo iz zgornje zveze:

a= 1

nσ. (4.15)

(48)

Izraˇcunajmo ˇstevilsko gostoto zraka pri temperaturi T = 273 K : n=ρ(T)·N_A

A = 1,2922 kg/m³6,022·10²⁶/kmol

28,9644/kmol = 2,687·10²⁵/m³. (4.16) (4.17) V simulaciji bom prevzela, da ima atmosfera na vseh viˇsinah enako temperaturo in tudi gostoto. Zanima me, kolikˇsno debelino ima v okviru te predpostavke.

Zrak je plin, zato uporabimo sploˇsno plinsko enaˇcbo:

pV = m

MRT . (4.18)

Ce v plinsko enaˇˇ cbo vstavimo zvezo za gostotoρ= ^m_V,dobimo:

p= ρ

MRT . (4.19)

Tlak v atmosferi lahko zapiˇsemo kot: p = ^F_S. Ker je F sila teˇze atmosfere debeline d, ki pritiska na povrˇsino S,jo izraˇcunamo kot:

F =mg =ρV g =ρdSg . (4.20)

Izraz za silo vstavimo v definicijo tlaka in v plinsko enaˇcbo (4.19). Dobimo izraz za debelino atmosfere v naˇsem pribliˇzku:

d= RT

gM . (4.21)

Izraˇcun debeline atmosfere pri povpreˇcni temperaturiT = 273 K : d= RT

gM = 8314 J/(K·kmol)·273K

9,81 m/s²·28,9644 kg/kmol = 7988 m = 8 km. (4.22) (4.23) Preko dneva se spreminja viˇsina sonca nad obzorjem in s tem tudi debelina atmosfere, ki je vmes. Njeno debelino lahko iz kosinusnega izreka izraˇcunamo kot:

x=p

R²_zcos²α+ 2R_zd+d²−R_zcosα , (4.24) kjer je kot α kot med potjox in polmerom Zemlje.

Ob sonˇcnem zahodu je kot αenak 90^◦. Ker je cos 90^◦ = 0, to debelino atmosfere izraˇcunamo kot:

d⁰⁰=p

2R_zd+d²=p

2·6378000 m·7988 m + (7988 m)²= 32010 m = 320 km. (4.25)

(49)

Tabela 4.1: Izraˇcuni sipalnega preseka in proste poti barva λ[nm] σ[10⁻²⁷cm²] a[km]

modra 470 9,016 41,278

zelena 550 4,808 77,405

rdeˇca 610 3,178 117,111

Iz tabele razberemo, da se z veˇcanjem valovne dolˇzine sipalni presek manjˇsa, povpreˇcna prosta pot pa poveˇcuje. V primerjavi z debelino atmosfere so povpreˇcne poti mnogo veˇcje, kar pomeni, da gre veˇcina fotonov direktno skozi in se sploh ne siplje. Njihov deleˇz lahko izraˇcunamo tako, da namesto poti vstavimo debelino atmosfere in izraˇcunamo vrednost spremenljivkeξ:

x=−alnξ0, (4.26)

ξ₀= e⁻^x^a. (4.27)

Ce ˇˇ zelimo v simulacijo vkljuˇciti le tiste fotone, ki se vsaj enkrat sipljejo, pri vzorˇcenju prve poti ˇstevilo ξˇzrebamo na intervalu [ξ0,1].

Tabela 4.2: Prilagajanje intervala ˇzrebanja prve opravljene poti za d= 8 km,d⁰ = 240 km ind⁰⁰= 320 km

barva λ[nm] ξ₀ ξ₀⁰ ξ₀⁰⁰ modra 470 0,824 0,559 0,0004 zelena 550 0,902 0,733 0,016 rdeˇca 610 0,934 0,815 0,065

Steviloˇ ξ0 nam torej pove, kolikˇsen deleˇz fotonov gre nemoteno skozi atmosfero. Deleˇz se s krajˇsanjem valovne dolˇzine pomanjˇsuje. Iz tabele vidimo, da so ti deleˇzi kar precejˇsnji ter da se s poveˇcevanjem debeline atmosfere zmanjˇsujejo. Tretja debelina atmosfere d⁰⁰ ustreza debelini atmosfere pri sonˇcnemu zahodu. Opazimo, da se deleˇz fotonov, ki grejo direktno skozi atmosfero, s krajˇsanjem valovne dolˇzine pribliˇzuje ˇstevilu 0. To pomeni, da se skoraj vsa svetloba siplje.

Poglejmo si sploˇsno odvisnost deleˇza ξ0 od valovne dolˇzine pri konstantni debelini d in odvisnost deleˇza ξ0 od poljubne debeline atmosfere x pri doloˇceni valovni dolˇzini svetlobe.

V enaˇcbo (4.27) vstavimo zvezi (4.15) in (4.12) in dobimo:

ξ₀(λ) = e⁻^x^a = exp

−τ₀nx λ⁴

. (4.28)

V enaˇcbo (4.27) vstavimo zvezo (4.24) in dobimo:

ξ0(x) = e⁻^x^a = exp −

pR²_zcos²α+ 2R_zd+d²−R_zcosα a

!

. (4.29)

(50)

Slika 4.3: Graf deleˇza direktne svetlobe v odvisnosti od valovne dolˇzine prid= 8 km.

Slika 4.4: Graf deleˇza direktne svetlobe v odvisnosti od valovne dolˇzine pri d⁰ = 24 km.

(51)

Slika 4.5: Graf deleˇza direktne svetlobe v odvisnosti od valovne dolˇzine pri d⁰⁰ = 320 km.

Slika 4.6: Graf deleˇza direktne svetlobe v odvisnosti od kota α

Rdeˇca krivulja predstavlja rdeˇco svetlobo z valovno dolˇzino 610 nm, zelena zeleno svetlobo z valovno dolˇzino 550 nm in modra modro svetlobo z valovno dolˇzino 470 nm.

(52)

Slika 4.7: Graf deleˇza direktne svetlobe v odvisnosti od kota α v razmerju z deleˇzem direktne rdeˇce svetlobe

Iz grafov na slikah 4.6 in 4.7 razberemo, da je deleˇz direktne rdeˇce svetlobe pri vseh kotihα,torej pri vseh debelinah atmosfere, veˇcji od deleˇzev direktne zelene in modre svetlobe. Seˇstevanje barv v modelu RGB nam da razliˇcne barve. Ob sonˇcnem zahodu, ko gre kotα proti 90^◦, se deleˇza modre in zelene direktne svetlobe zmanjˇsujeta glede na rdeˇco, kar pomeni, da sonce postaja vedno bolj rdeˇce barve.

4.2.2 Kotna porazdelitvena funkcija

Sipanje opisuje kotna porazdelitvena funkcija, ki nam pove verjetnost, da se bo foton sipal v smeri sipalnega kota θ glede na prvotno smer. Gostota porazdelitve po sipalnem kotu pri Rayleighovemu sipanju je: [10]

f(θ) = 2

3π(1 + cos²θ). (4.30)

Vpeljali bomo novo spremenljivko u = cos(θ); u ∈ [−1,1]. Gostota porazdelitve po spre- menljivkiu ima obliko:

f(u) =A(1 +u²), (4.31)

kjer je Anormalizacijska konstanta, ki jo doloˇcimo preko normalizacijskega pogoja.

(53)

Z 1

−1

f(z)du = 1 Z 1

−1

A(1 +u²)du = 1 8

3A = 1

A = 3

8 Dobimo porazdelitev s predpisom:

f(u) = 3

8(1 +u²), (4.32)

ki zavzame vrednosti na intervalu [³₈,³₄].

Te porazdelitve ne moremo enostavno obrniti kot v prejˇsnjih primerih, zato uporabimo metodo zavrnitve, ki je opisana v podpoglavju Vzorˇcenje 2.3.3:

1. Za primerjalno funkcijog(u) bomo vzeli kar konstantno funkcijo, ki ima vrednost najveˇcje vrednosti gostote porazdelitvef(u). Doloˇcimo:

g(u) = 3

4. (4.33)

2. Izˇzrebamo nakljuˇcno vrednost u_i iz intervala [−1,1].

3. Izˇzrebamo nakljuˇcno vrednost v_i iz intervala [0,³₄].Dobili smo nakljuˇcno toˇcko v ravnini s koordinatami (u_i, v_i).

4. • Ce velja:ˇ v_i≤f(v_i), potem vrednostu_isprejmemo in naprej izvajamo sipanje fotona.

• V nasprotnem primeru vrednostu_izavrnemo in postopek ˇzrebanja ponovimo, dokler ne dobimo ustrezne vrednosti.

(54)

Slika 4.8: Primer metode zavrnitve, kjer toˇcko T₁ zavrnemo, T₂ pa sprejmemo

4.2.3 Transformacija vektorja premika v prvotni koordinatni sistem

Vzeli bomo takˇsen koordinatni sistem, ki ima osz usmerjeno v smeri sipanja. Njegovo os x pa izberemo tako, da jeϕ₁= 0,kar pomeni, da se foton po prvem sipanju giblje v ravninix, z.

Po prvem sipanju se foton premakne za vektor:

−

→r1 = (r1sinθ1,0, r1cosθ1), (4.34) kjer je r1 izˇzrebana dolˇzina poti.

V primeru, da je foton ˇse vedno v atmosferi, se siplje ˇse enkrat. Toda foton se sedaj premakne glede na nov koordinatni sistem, ki ima os z usmerjeno v smeri sipanja. V tem sistemu se premakne za vektor:

−

→r₂ = (r₂cosϕ₂sinθ₂, r₂sinϕ₂sinθ₂, r₂cosθ₂). (4.35) Nas zanimata kot in premik glede na prvotni koordinatni sistem, pri katerem je osz usmer- jena v smeri vpadnega curka fotonov. Vektor−→r2 moramo zato zavrteti okoli osiy za sipalni kot prvega sipanja θ1.Uporabimo rotacijsko matriko:

Ry(θ1) =





cosθ1 0 sinθ1

0 1 0

−sinθ₁ 0 cosθ₁



. (4.36)