Compare commits

..

No commits in common. "master" and "a7a8e4fbcbcb5b0f4f7645a20eadb1f7758d053d" have entirely different histories.

69 changed files with 687 additions and 7952 deletions

20
.gitignore vendored
View File

@ -1,20 +0,0 @@
/modfizlab/fig1.png
*.png
*.pdf
*.log
*.html
/kutinf/elso
/kutinf/nagyhazi/field
/kutinf/nagyhazi/field.o
/kutinf/geometric-algebra
/kutinf/newton-cotes
/kutinf/root-finder
/modfizlab/221115.tar.xz
/modfizlab/221115/
/kutinf/nagyhazi/results.txt
*.npy
*.aux
*.out
*.pickle
/digilab/ccd/ccd.tex
*.xz

View File

@ -3,32 +3,4 @@
+ fiznum2/Intel.c -> experiment to apply the rk4 to guess cpu frequencies based on previous models.
+ fiznum2/lapakk.c -> program to benchmark lapack svd against home-grown gauss elimination
+ modfizlab/labor1.scm -> script to generate latex tables from measured values for a laboratory
+ kutinf -> programs for the korszám course
+ matrix:
#+begin_src bash
make matrix
./matrix
#+end_src
+ vector:
#+begin_src bash
make vector
./vector
#+end_src
+ newton-cotes:
#+begin_src bash
make newton-cotes
./newton-cotes
#+end_src
+ root-finder:
#+begin_src bash
make root-finder
./root-finder
#+end_src
+ nagyhazi/field:
You'll need GNU parallel to do this
#+begin_src bash
make field
make results.txt
make temp-magnetization.png
make temp-energy.png
#+end_src
+ kutinf -> a program to simulate electrostatic potential and fields

View File

@ -1,340 +0,0 @@
#+OPTIONS: tex:t
#+AUTHOR: Barna Zsombor
#+DATE: Mérés időpontja: 2023. 10. 09.
#+TITLE: CCD
#+LATEX_CLASS: report
#+LATEX_CLASS_OPTIONS: [10pt,pdftex]
#+LATEX_HEADER: \usepackage[magyar]{babel}
#+LATEX_HEADER: \hypersetup{colorlinks=true,linkcolor=blue, linktoc=all, filecolor=magenta, urlcolor=cyan, pdfstartview=FitB,}
#+LATEX_HEADER: \usepackage[T1]{fontenc}
#+LATEX_HEADER: \usepackage[utf8]{inputenc}
#+LATEX_HEADER: \usepackage{graphicx}
#+LATEX_HEADER: \usepackage{amsmath}
#+LATEX_HEADER: \usepackage{amssymb}
#+LATEX_HEADER: \usepackage{booktabs}
#+LATEX_HEADER: \usepackage{indentfirst}
#+LATEX_HEADER: \usepackage[a4paper]{geometry}
#+LATEX_HEADER: \usepackage{url}
#+LATEX_HEADER: \usepackage{titling}
#+LATEX_HEADER: \usepackage{multirow}
#+LATEX_HEADER: \usepackage{braket}
#+LATEX_HEADER: \urlstyle{same}
#+LATEX_HEADER: \usepackage{lmodern}
#+LATEX_HEADER: \usepackage{blindtext}
#+LATEX_HEADER: \usepackage{xcolor}
#+LATEX_HEADER: \definecolor{titlebg}{RGB}{128,128,128}
#+OPTIONS: toc:nil
#+BEGIN_EXPORT
\begin{titlepage}
\pagenumbering{gobble}
\maketitle
\end{titlepage}
\renewcommand{\thesection}{\arabic{section}}
\renewcommand*{\maketitle}{\noindent\colorbox{titlebg}{%
\parbox{\dimexpr\linewidth-2\fboxsep}{\color{white}%
\parbox{.4\linewidth}{\fontsize{24}{28}\selectfont\sffamily\bfseries\thetitle}\hfill%
\parbox{.4\linewidth}{\fontsize{12}{14}\selectfont\raggedleft\today\\\theauthor%
}}}\vskip 1em}
\maketitle
#+END_EXPORT
#+name: imports
#+begin_src python :session generic
import pickle
import numpy as np
import scipy.optimize
import scipy.stats
import matplotlib.pyplot as plt
from itertools import takewhile
#+end_src
Érdemes az első paragrafusban leírni, hogy mi a mérés célja, röviden vázolni, milyen műszereket használsz. Ne felejtsd el rögzíteni, a mérés körülményeinek fontos változóit.
\blindtext
Elképzelhető, hogy folyószövegben mutat jobban a tudnivalók felsorolása, de dönthetünk úgy is, hogy pontokba szedve soroljuk fel mondanivalónkat:
- két csatornás oszcilloszkóp,
- asztali digitális mérőműszer,
- dots
* 1. feladat: a CCD termosztátjának vizsgálata
** a) feladatrész
A kamera termisztorának egyszeri leolvasása:
#+BEGIN_SRC python
# Hőmérséklet kiolvasása
return st5.get_temperature()
#+END_SRC
#+RESULTS:
: -9.98321164661452
Mérés egy percen keresztül, másodpercenként, majd ábra rajzolása
#+BEGIN_SRC python
# Mintavételezés
temps = list()
for _ in range(60):
temps.append(st5.get_temperature())
time.sleep(1)
# Ábrarajzolás
temps_arr = np.array(temps)
t = linspace(0, len(temps_arr)-1, len(temps_arr))
plot(t, temps_arr, marker="x")
xlabel('t [s]')
ylabel("T [°C]")
#+END_SRC
#+ATTR-LATEX: :float hb! :width .9\textwidth
#+CAPTION: Sajnos a mérés végénél elvesztettem az erdeti adatokat, így csak egy stabil $-10°C$-os hőmérsékletről tudtam új ábrát készíteni.
#+NAME: fig:stable
[[./temps-stable.png]]
A következő adatokat azonban az még az eredeti adatsorból generáltam. Statisztika kiszámítása a mért adatokból:
#+BEGIN_SRC python :results verbatim
return f"""Átlagos hőmérséklet: {scipy.stats.tmean(temps)}°C
Szórás: {scipy.stats.tstd(temps)}°C"""
#+END_SRC
#+RESULTS:
: Átlagos hőmérséklet: 17.48172534645875 °C
: Szórás: 0.008123205502034303 °C
** b) feladatrész
#+begin_src python
# mérés
st5.set_temperature(-5)
temps_cooling = list()
for _ in range(90):
temps_cooling.append(st5.get_temperature())
time.sleep(1)
temps_cooling_arr = np.array(temps_cooling)
# ábra
t_cooling = linspace(0, len(temps_cooling_arr)-1, len(temps_cooling_arr))
plot(t_cooling, temps_cooling_arr, marker="x")
xlabel('t [s]')
ylabel("T [°C]")
#+end_src
#+ATTR-LATEX: :float hb! :width .9\textwidth
#+CAPTION: A hőmérséklet előszőr a cél alá megy, majd elulról konvergál -5°C-hoz
#+NAME: fig:cooling
[[./temps-cooling.png]]
A detektor hőmérséklete ránézésre kb. 80 másodperc elteltével lesz stabil.
Amit az ábrán láthatunk, az a túlhűtés jelensége. Vélhetően az az oka, hogy a detektálás és a hűtés helye között fizikai távolság van, és a hőterjedéshez idő kell. Ezt vezérlőelektronika úgy próbálja meg áthidalni, hogy nagyobbat hűt rá, hogy hamarabb közelítsen a célhőmérséklethez. Ellenkező esetben a konvergáláshoz sokkal többet kellene várni.
#+begin_src python
# mérés
st5.reset_temperature()
temps_free = list()
for _ in range(180):
temps_free.append(st5.get_temperature())
time.sleep(1)
temps_free_arr = np.array(temps_free)
# ábra
t_free = linspace(0, len(temps_free_arr)-1, len(temps_free_arr))
plot(t_free, temps_free_arr, marker="x")
xlabel('t [s]')
ylabel("T [°C]")
#+end_src
Az idei adataim egy inaktivitás miatti leállítás áldozatai lettek, ezért elővettem a tavaly mért adatokat, hiszen görbét illeszteni arra is lehet. ( mégha haszontalan is a mérés szempontjából, tanulni lehet belőle )
#+begin_src python :session generic
last = np.load("last-years.bin.npy")[10:]
t_values = np.arange(len(last), dtype=float)
def fit_exp(y):
expfunc = lambda t, A, K, limit: A * np.exp(K*t) + limit
guess_limit = np.array([np.max(y)]*len(y))
values = (guess_limit - y)
for index, element in enumerate(values):
if element == 0.0:
break
filtered = np.log(values[:index])
times = np.arange(len(filtered), dtype=float)
fit_data = np.polyfit(times, filtered, 1, cov=True)
K, A_log = fit_data[0]
cov = fit_data[1:]
A = np.exp(A_log)
return A, K, guess_limit[0], cov
A, K, limit, cov = fit_exp(last)
fitfunc = lambda t: -A * np.exp(K*t) + limit
plt.clf()
plt.plot(t, last, marker="x")
plt.plot(t, fitfunc(t))
plt.xlabel('t [s]')
plt.ylabel("T [°C]")
plt.savefig("convergence.png")
[round(A,2), round(K,2), round(limit,2), cov]
#+end_src
Az illesztett függvény tehát a következő
#+name: eq:fit
#+BEGIN_equation
T(t) = 19.73°C \cdot e^{-0.04t} + 24.54°C
#+END_equation
A kovarianciamátrix:
| 1.02603451e-07 | -8.25957783e-06 |
| -8.25957783e-06 | 0.000889281213 |
#+ATTR-LATEX: :float hb! :width .9\textwidth
#+CAPTION: A nyers hőmérsékleti adatok, és az illesztett függvény
#+NAME: fig:conv
[[./convergence.png]]
Az idei adatokból készült kép szerencsére nem veszett el, így meg tudjuk becsülni a szobahőmérsékletet, amin a következő feladatrész képe készült.
#+ATTR-LATEX: :float hb! :width .9\textwidth
#+CAPTION: A hőmérséklet bekonvergál a szobahőmérséklet közelébe, 22.5°C-hoz
#+NAME: fig:free
[[./temps-free.png]]
** Második feladat
4 másodperces záridővel, az előző feladatban "megállapított" határértéken mértem:
#+begin_src python
# mérés
st5.take_image(4)
img_black = st5.read_image()
# kép ábrázolása
np.save("img_black.bin", img_black)
imshow(img_black, cmap='gray')
#+end_src
#+ATTR-LATEX: :float hb! :width .9\textwidth
#+CAPTION: A hőmérséklet bekonvergál a szobahőmérséklet közelébe, 22°C-hoz
#+NAME: fig:black
[[./black.png]]
** b feladatrész
#+begin_src python :session generic
img_black = np.load("img_black.bin.npy")
# hisztogram elkészítése
hist = scipy.ndimage.histogram(img_black, min=0, max=2**14-1, bins=2**14)
# hisztogram ábra
plt.plot(hist)
plt.xlabel('pixelintenzitás [1]')
plt.ylabel("pixelszám [db]")
plt.xlim(150, 550)
plt.savefig("hist.png")
#+end_src
#+ATTR-LATEX: :float hb! :width .9\textwidth
#+CAPTION: A hisztogram egy részlete, a 150-es és az 550-es intenzitású binek közt. A széleken kívül alig esett pixel.
#+NAME: fig:hist
[[./hist.png]]
#+begin_src python :session generic :results verbatim
# átlagos érték és zaj
f"""Átlagos fényesség: {round(scipy.stats.tmean(img_black), 2)}
Szórás: {round(scipy.stats.tstd(img_black, axis=None),2)}"""
#+end_src
#+RESULTS:
: Átlagos fényesség: 298.69
: Szórás: 54.33
** c feladatrész
Az adatgyűjtéshez az alábbi segédfüggvényeket használtam.
#+begin_src python :session generic
images = dict()
def set_temp(t):
st5.set_temperature(t)
time.sleep(85)
print(st5.get_temperature())
def take_pix(at_temp):
set_temp(at_temp)
images[at_temp] = dict()
for exp_t in range(100,1200,200):
st5.take_image(exp_t)
images[at_temp][exp_t] = st5.read_image()
#+end_src
Az idei adatokat az újraindítás jóvoltából elvesztettem, azonban a tavalyikat az archívumból ki tudtam bányászni.
#+begin_src python :session generic
with open("temp-images.pickle", "rb") as fh:
d = pickle.load(fh)
#+end_src
#+RESULTS:
Lehet matematikai kifejezést is betenni a szövegbe:
\[
x_i(t) = a_{i1} s_1(t) + a_{i2} s_2(t).
\]
Ugyanakkor sokszor hasznos a kifejezésre a folyószövegből hivatkozni. Tehát a fenti kifejezés megegyezik alábbi .~egyenlettel:
#+NAME: eq:bla
#+BEGIN_equation
x_i(t) = a_{i1} s_1(t) + a_{i2} s_2(t).
#+END_equation
És természetesen a szövegben illik elmondani, hogy melyik matematikai kifejezéssel, mit jelölünk. Ha ismernénk $a_{ij}$-t, klasszikus módszerekkel meghatározhatnánk $s_j(t)$-t.
Ne felejtsétek el, ha referenciákat tartalmaz a dokumentum, akkor a fordítást több menetben kell megejtenetek. A \verb,pdflatex, kétszeri-háromszori futtatására van szükség, mert elsőre a referenciák helyközét számolja ki a kompiler, és táblázatot készít a behelyettesítendő mennyiségekről, és csak az újabb fordításkor helyettesíti be a korrekt értékeket. Ha erre nem figyelünk, akkor =??= jelet fogunk a dokumentumban látni.
Álljon itt még [[tab:bla]].~példatáblázat, amire szépen a folyó szövegből is hivatkozunk.
#+ATTR-LATEX: :align l|r
#+CAPTION: Ez egy példatáblázat. Az első oszlopban a mérés száma szerepel, a második oszlop légbőlkapott feszültségértékeket sorol fel.
#+NAME: tab:bla
| $n$ | $U [V]$ |
|-----+-------------|
| 1 | 12 $\pm$ .1 |
| 2 | 11 $\pm$ .1 |
Ha fordítás során a jegyzőkönyvbe illesztenő képet nem találja a fordító, tipikusan [[fig:bla]].~ábrán látható hibába botlunk.
#+ATTR-LATEX: :float hb! :width .9\textwidth
#+CAPTION: Mindig kell írni a kép alá értelmező szöveget. Itt azt lehet látni, milyen hibát mutat a fordító.
#+NAME: fig:bla
[[./valami.png]]
* Tapasztalatok összegzése
Hasznos a jegyzőkönyvet egy pár soros lezárással befejezni. Leírható, mit tanított a mérés, milyen ötleteket ébresztett, ha azok megosztásra érdemesek.
A \LaTeX-ben ebben a sablonban is használt =figure= és =table= környezetek úszó objektumokat valósítanak meg. Ha sok van belőlük, és a kiszedett folyószöveg relatív kevés területet foglal el az oldalon, akkor úgy járhatunk, hogy a táblázatok, és ábrák jelentős része a jegyzőkönyv végére csúszva sorakoznak fel. Ezen javít, ha billen a szövegterjedelem javára a mérleg. Azonban ha minden szükséges magyarázatot tartalmaz a jegyzőkönyv, és továbbra is fennáll ez az esztétikai probléma, ezen túllépünk, nem von le a munka értékéből.
* Felhasznált források
- [[https://bla.com][bla.com]]

View File

@ -1,97 +0,0 @@
\documentclass[10pt,pdftex]{report}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage[magyar]{babel}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{multirow}
\usepackage{booktabs}
\usepackage{indentfirst}
\usepackage[a4paper]{geometry}
\usepackage{url}
\usepackage{titling}
\usepackage{xcolor}
\definecolor{titlebg}{RGB}{128,128,128}
\usepackage{lmodern}
\usepackage{blindtext}
\renewcommand{\thesection}{\arabic{section}}
\title{Mérés címe}
\author{Gipsz Jakab}
\renewcommand*{\maketitle}{\noindent\colorbox{titlebg}{%
\parbox{\dimexpr\linewidth-2\fboxsep}{\color{white}%
\parbox{.4\linewidth}{\fontsize{24}{28}\selectfont\sffamily\bfseries\thetitle}\hfill%
\parbox{.4\linewidth}{\fontsize{12}{14}\selectfont\raggedleft\today\\\theauthor%
}}}\vskip 1em}
\begin{document}
\maketitle
Érdemes az első paragrafusban leírni, hogy mi a mérés célja, röviden vázolni, milyen műszereket használsz. Ne felejtsd el rögzíteni, a mérés körülményeinek fontos változóit.
\blindtext
Elképzelhető, hogy folyószövegben mutat jobban a tudnivalók felsorolása, de dönthetünk úgy is, hogy pontokba szedve soroljuk fel mondanivalónkat:
\begin{itemize}
\item{}két csatornás oszcilloszkóp,
\item{}asztali digitális mérőműszer,
\item{}\dots
\end{itemize}
\section{Első feladat}
Lehet matematikai kifejezést is betenni a szövegbe:
\[
x_i(t) = a_{i1} s_1(t) + a_{i2} s_2(t).
\]
Ugyanakkor sokszor hasznos a kifejezésre a folyószövegből hivatkozni. Tehát a fenti kifejezés megegyezik \aref{eq:bla}.~egyenlettel.
\begin{equation}
x_i(t) = a_{i1} s_1(t) + a_{i2} s_2(t).\label{eq:bla}
\end{equation}
És természetesen a szövegben illik elmondani, hogy melyik matematikai kifejezéssel, mit jelölünk. Ha ismernénk $a_{ij}$-t, klasszikus módszerekkel meghatározhatnánk $s_j(t)$-t.
Ne felejtsétek el, ha referenciákat tartalmaz a dokumentum, akkor a fordítást több menetben kell megejtenetek. A \verb,pdflatex, kétszeri-háromszori futtatására van szükség, mert elsőre a referenciák helyközét számolja ki a kompiler, és táblázatot készít a behelyettesítendő mennyiségekről, és csak az újabb fordításkor helyettesíti be a korrekt értékeket. Ha erre nem figyelünk, akkor \verb_??_ jelet fogunk a dokumentumban látni.
Álljon itt még \aref{tab:bla}.~példatáblázat, amire szépen a folyó szövegből is hivatkozunk.
\begin{table}[hbt!]
\begin{center}
\begin{tabular}{l|r}
$n$ & $U [V]$ \\
\hline
1 & 12 $\pm$ .1 \\
2 & 11 $\pm$ .1 \\
\end{tabular}
\caption{\label{tab:bla}Ez egy példatáblázat. Az első oszlopban a mérés száma szerepel, a második oszlop légbőlkapott feszültségértékeket sorol fel.}
\end{center}
\end{table}
Ha fordítás során a jegyzőkönyvbe illesztenő képet nem találja a fordító, tipikusan \aref{fig:bla}.~ábrán látható hibába botlunk.
\begin{figure}[hb!]
\includegraphics[width=.9\textwidth]{valami.png}
\caption{\label{fig:bla}Mindig kell írni a kép alá értelmező szöveget. Itt azt lehet látni, milyen hibát mutat a fordító.}
\end{figure}
\section*{Tapasztalatok összegzése}
Hasznos a jegyzőkönyvet egy pár soros lezárással befejezni. Leírható, mit tanított a mérés, milyen ötleteket ébresztett, ha azok megosztásra érdemesek.
A \LaTeX-ben ebben a sablonban is használt \verb_figure_ és \verb.table. környezetek úszó objektumokat valósítanak meg. Ha sok van belőlük, és a kiszedett folyószöveg relatív kevés területet foglal el az oldalon, akkor úgy járhatunk, hogy a táblázatok, és ábrák jelentős része a jegyzőkönyv végére csúszva sorakoznak fel. Ezen javít, ha billen a szövegterjedelem javára a mérleg. Azonban ha minden szükséges magyarázatot tartalmaz a jegyzőkönyv, és továbbra is fennáll ez az esztétikai probléma, ezen túllépünk, nem von le a munka értékéből.
% rutinosabban bibtexelhetnek, de ide most ez overkill
\section*{Felhasznált források}
\begin{itemize}
\item{}\url{https://bla.com}
\end{itemize}
\end{document}

View File

@ -1,27 +1,34 @@
# right now only gcc supports float variables as template parameters
CXX = clang++-15
CXX = g++ # egyelőre csakis gcc a beágyazott helper függvények miatt
DBGFLAGS = -Wall -Wextra -g
DBGFLAGS = -Wall -Wextra
#-DDEBUG -g
#PARAFLAGS= -DTHREADED -fopenmp -D_REENTRANT
#OPTIMFLAGS= -O1 -march=znver1 -mprefer-vector-width=256 -flto -fdevirtualize-at-ltrans -fno-plt -fno-common -fipa-pta -fno-semantic-interposition -fgraphite-identity -floop-nest-optimize
PARAFLAGS= -DVLEN=4 -DTHREADED -fopenmp -I/usr/include/SDL2 -D_REENTRANT
OPTIMFLAGS= -O1
# -march=znver1 -mprefer-vector-width=256 -flto -fdevirtualize-at-ltrans -fno-plt -fno-common -fipa-pta -fno-semantic-interposition -fgraphite-identity -floop-nest-optimize
CXXFLAGS = -std=c++20 $(DBGFLAGS)
# $(OPTIMFLAGS) $(PARAFLAGS)
LDFLAGS =
# -flto=6
CXXFLAGS = --std=c++17 $(DBGFLAGS) $(OPTIMFLAGS) $(PARAFLAGS)
LDFLAGS = -flto=6 -pthread -lSDL2
SRCS = $(wildcard *.cc)
OBJS = $(SRCS:.cc=.o)
HDRS = $(SRCS:.cc=.h)
all: field
# Map each object file to a program name (e.g. foo.o -> foo).
PROJECTS = $(basename $(OBJS))
field: $(OBJS)
$(CXX) $(CXXFLAGS) -o $@ $(OBJS) $(LDFLAGS)
%.o: %.cc
$(CXX) $(CXXFLAGS) -o $@ -c $< $(LDFLAGS)
all: $(PROJECTS)
TAGS: $(SRCS)
ctags -Re $(SRCS)
.PHONY: clean
clean:
rm -fv $(OBJS) $(PROJECTS)
rm -fv $(OBJS) field TAGS
.PHONY: lint
lint:
clang-tidy $(SRCS) $(HDRS) -- $(CXXFLAGS)

114
kutinf/field.cc Normal file
View File

@ -0,0 +1,114 @@
#include <iostream>
#include <array>
#include <SDL2/SDL.h>
#define RESY 600
#define RESX 800
using namespace std;
typedef array<array<double, 3>, 3> Stencil;
// we add a line of zeroes around the area to simplify the computation kernel
// this we won't need a series of ifs in the loop or breaking up the loop
typedef array<array<double, RESX+2>, RESY+2> Field;
inline void applyStencil(Field &current, Field &next, Stencil &stencil) {
double sum;
#pragma omp parallel for
for(unsigned int y=1; y < current.size()-1 ; y++) {
for(unsigned int x=1; x < current[y].size()-1 ; x++) {
// TODO this is slow and a bad abstraction
// TODO this is numerically bad
sum = 0.0;
sum+=current[y-1][x-1]*stencil[0][0];
sum+=current[y-1][x+0]*stencil[0][1];
sum+=current[y-1][x+1]*stencil[0][2];
sum+=current[y+0][x-1]*stencil[1][0];
sum+=current[y+0][x+0]*stencil[1][1];
sum+=current[y+0][x+1]*stencil[1][2];
sum+=current[y+1][x-1]*stencil[2][0];
sum+=current[y+1][x+0]*stencil[2][1];
sum+=current[y+1][x+1]*stencil[2][2];
next[y][x] = sum/4.0;
}
}
}
Field fields[2];
int _current_field=0;
#define current fields[_current_field]
#define next fields[(_current_field+1)%2]
inline void swap_fields() {
_current_field++;
_current_field%=2;
}
inline void add_points() {
current[10][10]=250.0;
current[59][79]=250.0;
current[34][18]=250.0;
current[35][18]=250.0;
current[34][19]=250.0;
current[35][19]=250.0;
current[400][500]=250.0;
}
int main(void) {
int maxiter = 2000;
Stencil laplacian;
laplacian[0]={0.0, 1.0, 0.0};
laplacian[1]={1.0, 0.0, 1.0};
laplacian[2]={0.0, 1.0, 0.0};
// laplacian[0]={0.0, 1.0, 0.0};
// laplacian[1]={1.0, 0.0, 1.0};
// laplacian[2]={0.0, 1.0, 0.0};
add_points();
for(int iter=0; iter < maxiter; iter++) {
applyStencil(current, next, laplacian);
swap_fields();
add_points();
}
SDL_Window *window;
SDL_Renderer *renderer;
SDL_Event event;
SDL_Init(SDL_INIT_VIDEO);
SDL_CreateWindowAndRenderer(RESX, RESY, 0, &window, &renderer);
SDL_SetRenderDrawColor(renderer, 255, 255, 255, 0);
SDL_RenderClear(renderer);
for(unsigned int y=1; y < current.size()-1; y++){
for(unsigned int x=1; x < current[y].size()-1; x++){
SDL_SetRenderDrawColor(renderer, 255-(Uint8)(current[y][x]) , 255-(Uint8)(current[y][x]), 255, 0);
SDL_RenderDrawPoint(renderer, x-1, y-1); //Renders on middle of screen.
}
}
SDL_RenderPresent(renderer);
while (1) {
if (SDL_PollEvent(&event) && event.type == SDL_QUIT)
break;
}
SDL_DestroyRenderer(renderer);
SDL_DestroyWindow(window);
SDL_Quit();
cout << current[400][500] << endl;
cout << current[401][500] << endl;
cout << current[402][500] << endl;
cout << current[403][500] << endl;
return EXIT_SUCCESS;
}

Binary file not shown.

View File

@ -1,210 +0,0 @@
#include <iostream>
#include <functional>
#include <cmath>
#include <array>
#include <complex>
#include "matrix.h"
using namespace std;
template<class T>
Matrix<T>& Matrix<T>::operator=(Matrix<T> const& cpy)
{
if(&cpy == this)
{
return *this;
}
data = cpy.data;
siz = cpy.siz;
return *this;
}
template<class T>
Matrix<T>&& Matrix<T>::operator=(Matrix<T>&& mv)
{
if(&mv == this)
{
return *this;
}
siz = mv.siz;
data = std::move(mv.data);
return *this;
}
template<class T>
Matrix<T> operator+(Matrix<T> const& left, Matrix<T> const& right)
{
return Matrix<T>([&](unsigned int x){ return left[x] + right[x]; },
left.size());
}
template<class T>
Matrix<T>&& operator+(Matrix<T>&& left, Matrix<T> const& right)
{
std::transform(left.data.begin(),
left.data.end(),
right.data.begin(),
left.data.begin(),
[](T const& x, T const& y){ return x + y; });
return std::move(left);
}
template<class T>
Matrix<T>& Matrix<T>::operator+=(Matrix<T> const& other)
{
std::transform(this->data.begin(),
this->data.end(),
other.data.begin(),
this->data.begin(),
[](T const& x, T const& y){ return x + y; });
return *this;
}
template<class T>
Matrix<T> operator-(Matrix<T> const& left, Matrix<T> const& right)
{
return Matrix<T>([&](unsigned int x){ return left[x] - right[x]; },
left.size());
}
template<class T>
Matrix<T>&& operator-(Matrix<T>&& left, Matrix<T> const& right)
{
std::transform(left.data.begin(),
left.data.end(),
right.data.begin(),
left.data.begin(),
[](T const& x, T const& y){ return x - y; });
return std::move(left);
}
template<class T>
Matrix<T>& Matrix<T>::operator-=(Matrix<T> const& other)
{
std::transform(this->data.begin(),
this->data.end(),
other.data.begin(),
this->data.begin(),
[](T const& x, T const& y){ return x - y; });
return *this;
}
template<class T>
Matrix<T> operator/(Matrix<T> const& A, T const& s)
{
return Matrix<T>([&](unsigned int i)
{
return A[i] / s;
},
A.size());
}
template<class T>
Matrix<T>&& operator/(Matrix<T>&& A, T const& s)
{
std::for_each(A.data.begin(),
A.data.end(),
[&](T & x){ return x /= s; });
return std::move(A);
}
template<class T>
Matrix<T>& Matrix<T>::operator/=(T const& s)
{
std::for_each(this->data.begin(),
this->data.end(),
[&](T & x){ return x /= s; });
return *this;
}
template<class T>
Matrix<T> operator*(Matrix<T> const& A, T const& s)
{
return Matrix<T>([&](unsigned int i)
{
return A[i] * s;
},
A.size());
}
template<class T>
Matrix<T>&& operator*(Matrix<T>&& A, T const& s)
{
std::for_each(A.data.begin(),
A.data.end(),
[&](T & x){ return x *= s; });
return std::move(A);
}
template<class T>
Matrix<T>& Matrix<T>::operator*=(T const& s)
{
std::for_each(this->data.begin(),
this->data.end(),
[&](T & x){ return x *= s; });
return *this;
}
template<class T>
Matrix<T> operator*(Matrix<T> const& left, Matrix<T> const& right)
{
auto N = left.size();
return
Matrix<T>([&](int i, int j)
{
T sum = static_cast<T>(0.0);
for(unsigned int k = 0; k < N; k++)
{
sum += left(i, k) * right(k, j);
}
return sum;
},
N);
}
template<class T>
std::ostream& operator<<(std::ostream& out, Matrix<T> const& arg)
{
auto N = arg.size();
for(unsigned int i = 0; i < N; i++)
{
for(unsigned int j = 0; j < N; j++)
{
out << arg(j, i) << "\t";
}
out << std::endl;
}
return out;
}
int main(int, char**)
{
auto M = Matrix<int>([](int i)
{
return i+1;
},
2);
auto N = Matrix<int>([](int i)
{
return i-1;
},
2);
cout << (M * 2) - (M / 2) * 2 << std::endl;
cout << (M * M) << std::endl;
cout << (M *= 2) << std::endl;
cout << (M -= N) << std::endl;
auto oof = Matrix<int>([](int i, int j)
{
return (i-1)*(j-1);
},
2) + N;
cout << N << std::endl;
cout << oof << std::endl;
return 0;
}

View File

@ -1,128 +0,0 @@
#include <vector>
#include <functional>
#ifndef WHOOM
#define WHOOM
template<class T>
class Matrix;
template<class T>
Matrix<T> operator* (Matrix<T> const& left, Matrix<T> const& right);
template<class T>
Matrix<T> operator+ (Matrix<T> const& left, Matrix<T> const& right);
template<class T>
Matrix<T>&& operator+ (Matrix<T>&& left, Matrix<T> const& right);
template<class T>
Matrix<T> operator- (Matrix<T> const& left, Matrix<T> const& right);
template<class T>
Matrix<T>&& operator- (Matrix<T>&& left, Matrix<T> const& right);
template<class T>
Matrix<T> operator* (Matrix<T> const& left, T const& right);
template<class T>
Matrix<T>&& operator* (Matrix<T>&& left, T const& right);
template<class T>
Matrix<T> operator/ (Matrix<T> const& left, T const& right);
template<class T>
Matrix<T>&& operator/ (Matrix<T>&& left, T const& right);
template<class T>
class Matrix
{
std::vector<T> data;
unsigned int siz;
public:
Matrix(unsigned int N)
{
this->data = std::vector<T>(N*N,
static_cast<T>(0));
this->siz = N;
}
Matrix(Matrix<T> const& cpy)
{
siz = cpy.siz;
data = cpy.data;
}
Matrix(Matrix<T>&& mv)
{
siz = mv.siz;
std::swap(data, mv.data);
}
Matrix(std::function<T(unsigned int)> f, unsigned int N)
{
data = std::vector<T>(N*N);
siz = N;
for(unsigned int k=0; k < (N*N); k++)
{
data[k] = f(k);
}
}
Matrix(std::function<T(unsigned int, unsigned int)> f, unsigned int N)
{
data = std::vector<T>(N*N);
siz = N;
for(unsigned int k = 0; k < N; k++)
{
for(unsigned int l = 0; l < N; l++)
{
data[k + N * l] = f(k, l);
}
}
}
Matrix<T>& operator=(Matrix<T> const& cpy);
Matrix<T>&& operator=(Matrix<T>&& mv);
//~Matrix();
T& operator[](int index)
{
return data[index];
}
T const& operator[](int index) const
{
return data[index];
}
T& operator()(unsigned int x, unsigned int y)
{
return data[x + this->siz * y];
}
T operator()(unsigned int x, unsigned int y) const
{
return data[x + this->siz * y];
}
unsigned int size() const
{
return siz;
}
friend Matrix operator* <>(Matrix const& left, Matrix const& right);
friend Matrix operator+ <>(Matrix const& left, Matrix const& right);
friend Matrix&& operator+ <>(Matrix&& left, Matrix const& right);
Matrix<T>& operator+=(const Matrix<T>&);
friend Matrix operator- <>(Matrix const& left, Matrix const& right);
friend Matrix&& operator- <>(Matrix&& left, Matrix const& right);
Matrix<T>& operator-=(const Matrix<T>&);
// scalar products
friend Matrix operator* <>(Matrix const& left, T const& right);
friend Matrix&& operator* <>(Matrix&& left, T const& right);
Matrix<T>& operator*=(T const&);
friend Matrix operator/ <>(Matrix const& left, T const& right);
friend Matrix&& operator/ <>(Matrix&& left, T const& right);
Matrix<T>& operator/=(T const&);
};
#endif

View File

@ -1,43 +0,0 @@
CXX = g++
DBGFLAGS = -Wall -Wextra
#-g
LIBFLAGS = -I/usr/include/SDL2 -D_REENTRANT
OPTIMFLAGS= -march=native -O2 -flto -fdevirtualize-at-ltrans -fno-plt -fno-common -fipa-pta -fno-semantic-interposition -fgraphite-identity -floop-nest-optimize
# -mprefer-vector-width=256
CXXFLAGS = --std=c++17 $(DBGFLAGS) $(LIBFLAGS) $(OPTIMFLAGS)
LDFLAGS = -pthread -lSDL2
#-flto=6
SRCS = $(wildcard *.cc)
OBJS = $(SRCS:.cc=.o)
HDRS = $(SRCS:.cc=.h)
all: field
field: $(OBJS)
$(CXX) $(CXXFLAGS) -o $@ $(OBJS) $(LDFLAGS)
%.o: %.cc
$(CXX) $(CXXFLAGS) -o $@ -c $< $(LDFLAGS)
results.txt: field
fish run.fish
temp-magnetization.png: results.txt
gnuplot magn.gnuplot
temp-energy.png: results.txt
gnuplot energy.gnuplot
.PHONY: clean
clean:
rm -fv $(OBJS) field
.PHONY: lint
lint:
clang-tidy $(SRCS) $(HDRS) -- $(CXXFLAGS)

View File

@ -1,10 +0,0 @@
#!/usr/bin/env gnuplot
set encoding utf8
set terminal png size 800,600 enhanced font 'IBM Plex Sans Text,14'
set xlabel 'Temperature'
set ylabel 'Total energy'
set output "temp-energy.png"
plot "results.txt" using 1:2 title "simulated energies" with points ps 1 pt 2 lw 2 lt rgb "#007faf"

View File

@ -1,170 +0,0 @@
#include <iostream>
#include <array>
#include <random>
#include <cmath>
#include <SDL2/SDL.h>
#include <utility>
#include <cstdlib>
#define RESY 100
#define RESX 100
#define MAXITER (RESX*RESY*iter_multiplier)
using namespace std;
static double T = 245.0;
#define phys_beta (1.0/T)
static double J = 1.0;
static double h = 0.0;
enum Spin {
up = 1,
down = -1,
};
inline void flip(Spin &spin)
{
spin = (spin == up) ? down : up;
}
// we add a line of zeroes around the area to simplify the computation kernel
// this way, we won't need a series of ifs in the loop or breaking up the loop
typedef array<array<Spin, RESX>, RESY> Field;
double get_energy_difference_for_flip(const Field& field, unsigned i, unsigned j)
{
int sum = 0;
// one up
sum += (i != 0) ? field[i - 1][j] : 0;
// one down
sum += (i != RESX - 1) ? field[i + 1][j] : 0;
// one left
sum += (j != 0) ? field[i][j - 1] : 0;
// one right
sum += (j != RESY - 1) ? field[i][j + 1] : 0;
return 2 * J * field[i][j] * sum + 2 * h * field[i][j];
}
pair<double, double> get_results(const Field& field)
{
int mag = 0;
int pairwise = 0;
// everthing except the last column and row
for(int y = 0; y < (RESY - 1); y++)
{
for(int x = 0; x < (RESX - 1); x++)
{
pairwise += field[y][x]*field[y + 1][x];
pairwise += field[y][x]*field[y][x + 1];
mag += field[y][x];
}
}
// last column
for(int y = 0; y < (RESY - 1); y++)
{
pairwise += field[y][RESX - 1]*field[y + 1][RESX - 1];
mag += field[y][RESX - 1];
}
// last row
for(int x = 0; x < (RESX - 1); x++)
{
pairwise += field[RESY - 1][x]*field[RESY - 1][x + 1];
mag += field[RESY - 1][x];
}
// last field
mag += field[RESY - 1][RESX - 1];
double magnetization = static_cast<double>(mag) / RESX / RESY;
return pair(- J * pairwise + h * mag, magnetization);
}
int main(int argc, char *argv[] ) {
T = atof(argv[1]);
int iter_multiplier = atoi(argv[2]);
//h = atof(argv[1]);
// stolened right from the reference
random_device rd;
mt19937 gen(rd());
uniform_int_distribution<> idis(0, RESX * RESY - 1);
uniform_real_distribution<> fdis(0.0, 1.0);
uniform_int_distribution<> spindis(0, 1);
Field field;
for(int y = 0; y < RESY; y++)
{
for(int x = 0; x < RESX; x++)
{
field[y][x] = ((spindis(gen) == 1) ? up : down);
}
}
double dE;
int xy, x, y;
double threshold;
for(int i = 0; i < MAXITER; i++)
{
// maybe i could use 2 random numbers but I can't be arsed to created 2 uniform generators
xy = idis(gen);
x = xy % RESX;
y = xy / RESX;
dE = get_energy_difference_for_flip(field, y, x);
if(dE < 0)
{
flip(field[y][x]);
continue;
}
threshold = fdis(gen);
if(exp(- phys_beta * dE) < threshold)
{
flip(field[y][x]);
continue;
}
}
// SDL_Window *window;
// SDL_Renderer *renderer;
// SDL_Event event;
// SDL_Init(SDL_INIT_VIDEO);
// SDL_CreateWindowAndRenderer(RESY, RESX, 0, &window, &renderer);
// SDL_SetRenderDrawColor(renderer, 255, 255, 255, 0);
// SDL_RenderClear(renderer);
// for(unsigned int y = 0; y < RESY; y++){
// for(unsigned int x = 0; x < RESX; x++){
// SDL_SetRenderDrawColor(renderer,
// (field[y][x] == up) ? 255 : 0,
// (field[y][x] == up) ? 0 : 255,
// 0,
// 0);
// SDL_RenderDrawPoint(renderer, y, x); //Renders on middle of screen.
// }
// }
// SDL_RenderPresent(renderer);
// while (1) {
// if (SDL_PollEvent(&event) && event.type == SDL_QUIT)
// break;
// }
// SDL_DestroyRenderer(renderer);
// SDL_DestroyWindow(window);
// SDL_Quit();
pair<double, double> result = get_results(field);
cout << T << "\t" << get<0>(result) << "\t" << get<1>(result) << endl;
return EXIT_SUCCESS;
}

View File

@ -1,11 +0,0 @@
#!/usr/bin/env gnuplot
set encoding utf8
set terminal png size 800,600 enhanced font 'IBM Plex Sans Text,14'
set xlabel 'Temperature'
set ylabel 'Magnetization'
set output "temp-magnetization.png"
plot "results.txt" using 1:3 title "simulated magnetizations" with points ps 1 pt 2 lw 2 lt rgb "#af007f"

View File

@ -1,5 +0,0 @@
#!/usr/bin/env fish
parallel ./field {} 10 >> results.txt ::: (python -c "for x in range(1, 501):
print(x/100)")

View File

@ -1,78 +0,0 @@
#include <iostream>
#include <functional>
#include <tuple>
#include <cmath>
using namespace std;
template<class T, unsigned int N=2u>
T adaptive_newton_cotes(function<T(T)> integrand,
function<T(function<T(T)>, T, T)> approximation,
T lower,
T upper,
T tolerance)
{
static_assert(N > 1, "N must at least 2");
T sum = static_cast<T>(0.0);
T h = (upper - lower) / N;
T interval_lower, interval_upper, interval_middle;
T interval_integral, lower_half_integral, upper_half_integral;
// divide the whole interval into N subintervals
for(unsigned int i = 0; i < N; i++)
{
// divide the subinterval into 2 subdivisions
interval_lower = lower + i * h;
interval_middle = lower + (i+0.5) * h;
interval_upper = lower + (i+1) * h;
interval_integral = approximation(integrand, interval_lower, interval_upper);
lower_half_integral = approximation(integrand, interval_lower, interval_middle);
upper_half_integral = approximation(integrand, interval_middle, interval_upper);
cout << "Margin of error on interval [" << interval_lower << "\t.. " << interval_upper << "\t]: " <<
fabs(interval_integral-(lower_half_integral+upper_half_integral)) << endl;
// compare the subinterval to the 2 subdivisions
if(fabs(interval_integral-(lower_half_integral+upper_half_integral)) <= tolerance)
{
sum += interval_integral;
}
else
{
cout << "Dividing the interval..." << endl;
sum += adaptive_newton_cotes<T>(integrand,
approximation,
interval_lower,
interval_upper,
tolerance);
}
}
// dunno why but it calculates double...
return sum;
}
template<class T>
T simpson_approximation(function<T(T)> f, T lower, T upper)
{
return 1.0/6.0 * (upper - lower) * (f(lower) + 4*f((upper+lower)/2) + f(upper));
}
template<class T>
T trapezoidal_approximation(function<T(T)> f, T lower, T upper)
{
return (upper - lower) * (f(lower) + f(upper)) / 2.0;
}
int main(int, char**)
{
auto poly = [](double z)->double {return -z*z*z - z*z + 4*z - 6;};
//auto test = [](double z)->double {return z*z*(z + 1);};
auto integral = adaptive_newton_cotes<double, 20u>(poly, trapezoidal_approximation<double>, 0.1, 1.0, 0.1);
cout << "The integral on [0.1, 1.0] is: " << integral << endl;
return 0;
}

View File

@ -1,62 +0,0 @@
#include <iostream>
#include <functional>
#include <cmath>
#include <limits>
const double epsilon = 0.01;
using namespace std;
template <class T, int max_iterations = 1000>
T newton_iteration(function<T(T)> f, function<T(T)> derivative, function<bool(function<T(T)>, T, T)> stop_predicate, T x0)
{
T x=x0;
// technically we don't need this
T prev_x = numeric_limits<T>::infinity();
int i = 1;
// I wanted the loop to run at least once
do
{
if(i > max_iterations) throw domain_error("Iteration didn't reach a stable point.");
if(derivative(x) == 0) throw domain_error("Function is constant at the point.");
prev_x = x;
// if f is linear, then dx*f(x)/dx = f(x) and f(x) instantly becomes zero
x-=f(x)/derivative(x);;
cout << "Iteration " << i << ":\t" << "x = " << x << "\tdfx = " << derivative(x) << endl;
i++;
} while(stop_predicate(f,x, prev_x));
return x;
}
template <class T>
bool metric_predicate(function<T(T)> f, T prev, T curr)
{
return (fabs(f(curr)-f(prev)) > epsilon);
}
double sqrt_newton(double k)
{
auto sqrt_finder = [&k](double z)->double {return z*z - k;};
auto derivative = [](double z)->double {return 2*z;};
return newton_iteration<double, 1000>(sqrt_finder, derivative, metric_predicate<double>, +40.0);
}
int main(void)
{
// this has root at z=-3 and no other real roots
auto poly = [](double z)->double {return -z*z*z - z*z + 4*z - 6;};
auto dpoly = [](double z)->double {return -3*z*z - 2*z + 4;};
double x0 = newton_iteration<double, 1000>(poly, dpoly, metric_predicate<double>, +40.0);
cout << "Function nears zero around x = " << x0 << endl;
double k = 321;
double root = sqrt_newton(k);
cout << "The square root of " << k << " is " << root << endl;
return 0;
}

View File

@ -1,101 +0,0 @@
#include <iostream>
#include <functional>
#include <cmath>
#include <array>
#include <complex>
#include "vector.h"
using namespace std;
template<class T>
Vector2D<T>& Vector2D<T>::operator+=(const Vector2D<T>& other)
{
coords[0] += other.coords[0];
coords[1] += other.coords[1];
return *this;
};
template<class T>
Vector2D<T>& Vector2D<T>::operator-=(const Vector2D<T>& other)
{
coords[0] -= other.coords[0];
coords[1] -= other.coords[1];
return *this;
};
// prolly only for commutative product, i don't think it's needed because can't be chained
template<class T>
Vector2D<T>& Vector2D<T>::operator*=(const T& scalar)
{
coords[0] *= scalar;
coords[1] *= scalar;
return *this;
};
template<class T>
Vector2D<T> operator+(const Vector2D<T>& left, const Vector2D<T>& right)
{
return Vector2D<T>(left.coords[0] + right.coords[0], left.coords[1] + right.coords[1]);
};
template<class T>
Vector2D<T> operator-(const Vector2D<T>& left, const Vector2D<T>& right)
{
return Vector2D<T>(left.coords[0] - right.coords[0], left.coords[1] - right.coords[1]);
};
template<class T>
Vector2D<T> operator*(const T& scalar, const Vector2D<T>& vector)
{
return Vector2D<T>(scalar * vector.coords[0], scalar * vector.coords[1]);
};
template<class T>
ostream& operator<<(ostream& os, const Vector2D<T>& vector)
{
return os << "(" << vector.coords[0] << ", " << vector.coords[1] << ")" ;
};
template<class T>
istream& operator>>(istream& is, Vector2D<T>& vector)
{
is >> vector.coords[0];
is >> vector.coords[1];
return is;
};
template<class T>
T euclidean_norm(const Vector2D<T>& v)
{
return dot_product(v,v);
}
template<class T>
T dot_product(const Vector2D<T>& vl, const Vector2D<T>& vr)
{
return
(vl.coords[0] * conj(vr.coords[0])) +
(vl.coords[1] * conj(vr.coords[1]));
}
template<class T>
Vector2D<T> normalized(const Vector2D<T>& v, function<T(const Vector2D<T> &)> norm_function)
{
T norm = norm_function(v);
return (static_cast<T>(1)/norm) * v;
}
int main(int, char**)
{
using namespace std::complex_literals;
typedef complex<double> C;
Vector2D<C> first(1. + 2i, 3. + 1i);
Vector2D<C> second(-1. - 2i, 3. + 1i);
C fugg(-1.0 - 0.5i);
cout << euclidean_norm<C>(fugg * first + second) << endl;
cout << dot_product(first, second) << endl;
cout << normalized(first+second, euclidean_norm<C>) << endl;
return 0;
}

View File

@ -1,21 +0,0 @@
#include <array>
#ifndef WHOOO
#define WHOOO
template<class T>
struct Vector2D
{
std::array<T,2> coords;
Vector2D(): coords(std::array<T,2>{0.0,0.0}) {}
Vector2D(const T& x, const T& y): coords(std::array<T,2>({x,y})) {}
Vector2D(const std::array<T,2>& xy): coords(xy) {}
Vector2D<T>& operator+=(const Vector2D<T>&);
Vector2D<T>& operator-=(const Vector2D<T>&);
// prolly only for commutative product
Vector2D<T>& operator*=(const T&);
};
#endif

BIN
modfizlab/1.png Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 68 KiB

172
modfizlab/15anna.tex Normal file
View File

@ -0,0 +1,172 @@
\documentclass[12pt,a4paper]{article}
\usepackage[utf8]{inputenc}
\usepackage[magyar]{babel}
\usepackage[T1]{fontenc}
\usepackage{graphicx}
\usepackage{epstopdf}
\usepackage{caption}
\usepackage{subcaption}
\usepackage{xcolor}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{multirow}
\usepackage{hyperref}
\usepackage[noadjust]{cite}
\usepackage[all]{hypcap}
\usepackage[hypcap]{caption}
\usepackage{braket}
\usepackage{url}
\usepackage{float}
\usepackage{listings}
\setlength{\textwidth}{17cm}
\setlength{\textheight}{24cm}
\setlength{\topmargin}{-2cm}
\setlength{\footskip}{1cm}
\setlength{\evensidemargin}{-0.3cm}
\setlength{\oddsidemargin}{-0.3cm}
\setlength{\parindent}{0cm}
\hypersetup{
colorlinks=true,
linkcolor=blue,
linktoc=all,
filecolor=magenta,
urlcolor=cyan,
pdfstartview=FitB,
}
\urlstyle{same}
\title{\Huge{\textbf{15. Kvantumradír}}}
\author{\textsl{Barna Zsombor, Rédl Anna}}
\date{Mérés időpontja: 2021.02.23.\\
\textsl{Keddi csoport}\\}
\begin{document}
\begin{titlepage}
\pagenumbering{gobble}
\maketitle
\end{titlepage}
\section{Bevezetés}
\pagenumbering{arabic}
\setcounter{page}{1}
\section{A mérési összeállítás}
\section{Rövid elméleti összefoglaló}
\section{Számolási feladatok}
\subsection{A lézernyalábok által bezárt szög számítása Mach-Zender elrendezésben}
Ha a négy tükör nem tökéletesen párhuzamos, a lencsén áthaladva a nyalábok a fókuszsíkban nem egyetlen pontban metszik egymást, hanem kettő pontban ( $S'$ és $S''$).
Ha feltételezzük, hogy a nyalábok $\alpha$ szöge nagyon kicsi, akkor az $S'$ és $S''$ pontok, valamint a lencse fősíkjának az $S'S''$ szakasz felezőjével szemközti pontja által definiált egyenlő szárú háromszög szárai $\approx f$, ahol $f$ a lencse fókusztávolsága.
Ebből a cosinus tétel alapján az $S'S''$ szakasz $d$ hosszúsága az alábbiak szerint becsülhető:
\begin{equation}
d^2\approx 2f^2(1-\cos\alpha)
\end{equation}
$\cos{\alpha}$ -t másodrendig sorbafejtve:
\begin{equation*}
\cos{\alpha}\approx 1-\frac{\alpha^2}{2}+\mathcal{O}(\alpha^4)
\end{equation*}
Ezt az előző kifejezésbe behelyettesítve:
\begin{equation*}
d=f\alpha
\end{equation*}
A $d$ szakasz hosszára kaphatunk egy másik egyenletet is az alaphán, hogy az $S'$ és $S''$ pontokból gömbhullámok indulnak ki. A $\theta_m$ irányok, az $m$ rendek, a $d$ réstávolság és a $\lambda$ hullámhossz között az alábbi összefüggés írható fel:
\begin{equation}
d\sin\theta_m=m\lambda
\end{equation}
Ha $\theta_m$ kicsi és a lencse és az ernyő $L$ távolságára igaz, hogy $L\gg f$, akkor élhetünk az alábbi közelítéssekkel a $\theta_m$ irányra és az adott interferenciacsíknak a középső csíktól való $l_m$ távolságára vonatkozóan:
\begin{equation*}
\theta \approx \sin\theta\approx\tan\theta\approx\frac{l_m}{L}
\end{equation*}
Ez alapján az egyenlet a következő formát ölti:
\begin{equation*}
d \approx \frac{m\lambda L}{l_m}
\end{equation*}
Ebbe behelyettesítve a $d$-re kapott előző kifejezést, a nyalábok párhuzamosságának $\alpha$ hibájára a következő kifejezés adódik:
\begin{equation}\label{eq: alpha}
\alpha\approx \frac{m\lambda L}{fl_m}
\end{equation}
\section{Mérési adatok}
\subsection{Jelölt és jelöletlen útvonal}
\begin{figure}[h]
\centering
\begin{minipage}{.5\textwidth}
\centering
\includegraphics[width=.9\linewidth]{1.png}
\captionof{figure}{Jelöletlen útvonal}
\label{fig:test1}
\end{minipage}%
\begin{minipage}{.5\textwidth}
\centering
\includegraphics[width=.9\linewidth]{2.png}
\captionof{figure}{Jelölt útvonal}
\label{fig:test2}
\end{minipage}
\end{figure}
\subsection{A kvantumradírozás szemléltetése}
\begin{figure}[h]
\centering
\includegraphics[width=8cm]{3.png}
\captionof{figure}{Bekapcsolt radírozó polárszűrő}
\end{figure}
\pagebreak
\subsection{Egyik ernyőn nincs jel}
\begin{figure}[h]
\centering
\begin{minipage}{.5\textwidth}
\centering
\includegraphics[width=.9\linewidth]{4.png}
\captionof{figure}{Másik ernyőn nincs interferencia}
\label{fig:test3}
\end{minipage}%
\begin{minipage}{.5\textwidth}
\centering
\includegraphics[width=.9\linewidth]{5.png}
\captionof{figure}{Másik ernyőn van interferencia}
\label{fig:test4}
\end{minipage}
\end{figure}
\section{A mérés kiértékelése}
\begin{thebibliography}{1}
\bibitem{leiras}
Koltai János: Kvantumradír
(\href{http://wigner.elte.hu/koltai/labor/parts/modern15.pdf}{Mérési leírás})
\end{thebibliography}
\end{document}

View File

@ -1,40 +0,0 @@
$T = 26°C$
$p = 1001 hPa$
$p2 = 999 hPa$
15 egységenként mérünk
#+NAME: example-table
| tér nélkül[s] | térrel[s] | tér[V] |
| 4.70 | 3.62 | 508 |
| 12.56 | 6.60 | 508 |
| 19.77 | 14.90 | 508 |
| 8.26 | 12.43 | 506 |
| 8.37 | 7.21 | 506 |
| 11.03 | 6.47 | 506 |
| 9.67 | 3.67 | 506 |
| 12.87 | 3.89 | 507 |
| 4.61 | 35.76 | 507 |
| 13.86 | 10.02 | 506 |
| 7.33 | 11.83 | 506 |
| 6.69 | 6.79 | 506 |
| 9.43 | 9.41 | 506 |
| 8.59 | 4.92 | 506 |
| 20.13 | 17.81 | 506 |
| 11.66 | 7.93 | 506 |
| 7.27 | 10.93 | 506 |
| 5.69 | 9.46 | 505 |
| 11.52 | 6.59 | 505 |
| 23.49 | 11.75 | 505 |
| 10.61 | 7.69 | 505 |
| 6.49 | 8.25 | 505 |
| 7.37 | 22.08 | 506 |
| 7.93 | 8.86 | 505 |
| 26.77 | 11.42 | 507 |
| 21.59 | 13.86 | 505 |
#+begin_src emacs-lisp :var
#+end_src

BIN
modfizlab/2.png Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 68 KiB

BIN
modfizlab/3.png Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 65 KiB

BIN
modfizlab/4.png Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 86 KiB

BIN
modfizlab/5.png Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 86 KiB

View File

@ -1,41 +0,0 @@
* A hidrogén és az alkálifémek spektumának történelmi jelentősége
Az esszém célja az, hogy megmutassam, hogy hogyan vezetett a spektroszkópia, és ezen belül az első főcsoport elemeinek a színképe a kvantummechanika kialakulásához.
** Rövid összefoglaló a történelmi előzményekről
A spektroszkópia tudománya rendkívül termékenynek bizonyult a történelem során. Forradalmasította a kémia anyagok vizsgálatát, hasznos segédeszköze lett a csillagászatnak, elválaszthatatlan a kvantummechanika felfedezésétől, és máig egy rendkívül elterjedt vizsgálati módszer. A története rendkívül korán elkezdődött, már ókori római írók (pl. Plinius) is említettek tárgyakat, amelyekben ma a prizmát véljük felfedezni. A spektrum mai nevét közismerten Newton adta, utána Fraunhofer tett nagy előrelépéseket: felvette a Nap spektrumát és észrevett benne vonalakat amelyek hiányoznak, ő készítette el az első diffrakciós rácsot, továbbá ő fedezte fel, hogy a hevített anyagok csakis bizonyos hullámhossztartományokban bocsájtanak ki fényt. Tiszteletből róla neveztük el a Fraunhofer-diffrakciót.
Fraunhofert követően először prizmákat használtak a folytonos színkép felvételére, de a tudomány hamar áttért a rácsok használatára, mert azokkal pontosabb munkát lehetett végezni. A spektroszkópia népszerűségének a növekedése magával vonta azt is, hogy 19. században megnőtt az igény az egyre kisebb rácsállandójú diffrakciós rácsokra. Ebben a versenyben egy neves magyar fizikus is képviseltette magát: Jedlik Ányos. Jedlik a rácsai révén Európa-szerte híressé vált.
A hidrogén színképének a felvétele először Ångström nevéhez köthető, ő vetette fel először, hogy anyagok ugyanazokban a sávokban nyelvnek el fényt, mint amelyekben kibocsájtanak a vonalas spektrumba.
A hevített kémai elemek sprektumának a szisztematikus vizsgálatát Robert Bunsen nevéhez szoktuk kötni. Ő és kollégája, Gustav Kirchoff úttörők voltak a spektroszkópia terén, és egy remek példa az interdiszciplinaritás hasznosságára. Ketten igazolták Ångströmnek az emissziós és abszorpciós vonalakra vonatkozó állítását. Kirchoff fizikus volt az eredmények alapján ő fogalmazta a spektroszkópia három törvényét:
- hevített szilárd anyagok és folyadok, illetve nagynyomású gázok folytonos spektrumban világítanak
- alacsony nyomású gázok vonalas színképet produkálnak
- ha egy folytonos színképet alacsony nyomású gázon átnézve szemlélünk, egy elnyelési vonalakat fogounk látni a folytonos spektrumban
Bunsen egy, az elődeinél jobb hevítőeszközt alkotott (a Bunsen-égőt), amellyel könnyebben tudták anyagok spektrumát vizsgálni. Bunsen és Kirchoff voltak az első, akik felvették a lítium, nátrium és a kálium spektrumát, illetve a módszer segítségével felfedezték a céziumot és a rubídiumot. Ezen kívül ketten bebizonyították, hogy a spektroszkópia alkalmas olyan esetekben is, amikor az adott elem csak nyomokban fedezhető fel a vizsgált anyagban.
Ők voltak az elsők abban is, hogy felvették egy olyan csillag színképét az éjjeli égboltról. Ezt ők maguk is annyira nagy előrelépésnek élték meg, hogy Bunsen a görög mitológiára hivatkozva így fogalmazott: "megloptuk az isteneket". Az Kirchoff és Bunsen az eredményeiket közösen publikálták.
** A kvantummechanikához vezető lépések
Johann Balmer volt az, akinek eszébe jutott, hogy a hidrogén színképében szereplő hullámhosszokat felírja egy képlettel, amiben az egyik paraméter a természetes számok közül kerül ki. A képlete jóslatai nem csak arra voltak hasznosak, hogy a mérési eredmények pontatlanságára rámutassanak, de aképlet néhány, még nem felfedezett vonalat is megjósolt. Sajnos ennek a képletnek a magyarázó ereje nem volt túl nagy, de Rydberg egy átrendezéssel olyan formára alakította a képletet, amelyet könnyebb volt általánosítani, illetve az értelmezése a mai kvantummechanika eredményeivel triviális.[fn:1] A képletében szereplő konstanst t Rydberg-állandónak neveztük el. Rydberg képlete jól modellezte a hidrogén és az alkálifémek spektrumát[fn:2], az általánosítását Walter Ritz készítette el. Ez akkor még szintúgy egy tisztán empirikus képlet volt.
Ezen korábbi eredmények alapján Bohr megalkotta az atommodelljét és megmutatta az első négy kémiai elem (H, He, Li, Be) esetében azt, hogy hogy kell a modellt alkamazni. Ezt modellt tanítják ma középiskolákban.
[fn:1] A trivialitás ténye szerintem jól szemlélteti, hogy mennyire fontos volt ez az eredmény a kvantummechanika kialakításában. Ha nem ehhez "illesztették" volna a modellt, akkor előfordulhat, hogy pl. perturbáció útján, valamely mellékfejezetben tárgyalnánk ma a témakört.
[fn:2] Az alkálifémek spektrumának levezetéséhez már kellenek a relativisztikus korrekciók és a spin-pálya kölcsönhatás figyelembe vétele, tehát a Dirac-egyenletből indulunk ki, de ezek elsősorban a szétvált vonalak modellezéséhez szükségesek: maguk a dublette a megjósolt hely közelében helyezkednek el.
** Spektroszkópia ma
** Tárgyalás
A történelmi előzmények rámutatnak, hogy azért az első főcsoport elemeinek színképe volt úttörő jelentőségű, mert az egyszerű elektronszerkezetük miatt egyrészt a színkép struktúrája egyszerű lesz, másrészt a kvantummechanika egyenletei viszonylag könnyen jóslatokat produkálnak, amelyeket utána ellenőrizni tudunk. Az új fizikai elméletek első igazoló kísérletei rendre ezt a két tulajdonságot ötvözik: ha az mérési hiba nagy, illetve a kép zavaros, akkor bármilyen egyszerű is a modell, a kísérletben nehézkesen lesz használható. Ha azonban maga az rendszer bonyolult, amit mérünk, akkor a még gyerekcipőben járó elmélet nincs eléggé kidolgozva ahhoz, hogy hamar igazolható jóslatokat tegyünk. A kísérlet megválasztása tehát a tudomány inkrementális tudásképző modelljébe illeszkedik, és kulcsfontosságú például akkor, ha egy fizikus demonstrálni kívánja laikusok számára, hogy a kvantummechanika nem a semmiből jött vagy néhány zseni fejében pattant ki és köze sincs a valósághoz, illetve nem puszta értelmiségi hóbort, hanem egy történelmileg mélyen beágyazott, kísrletekkel jól alátámasztott, a magyarázatokhoz szükséges modell. Szeretem idézni az alábbi (*alátámasztatlanul* Einstein-nek tulajdonított https://quoteinvestigator.com/2011/05/13/einstein-simple/) megfogalmazást Occam borotvapengéjéről: "Everything should be made as simple as possible, but not simpler." A kísérlet az új modellek szükségességét támasztja alá.
Bunsen-Kirchoff https://books.google.com/books?id=3KwRAAAAYAAJ
Ritz https://gallica.bnf.fr/ark:/12148/bpt6k33824.image.f185

View File

@ -1,130 +0,0 @@
% Created 2023-05-01 h 14:22
% Intended LaTeX compiler: pdflatex
\documentclass[11pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage[magyar]{babel}
\usepackage{graphicx}
\usepackage{longtable}
\usepackage{wrapfig}
\usepackage{rotating}
\usepackage[normalem]{ulem}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{capt-of}
\usepackage{hyperref}
\author{Barna Zsombor}
\date{\today}
\title{A hidrogén és az alkálifémek spektruma}
\hypersetup{
pdfauthor={Barna Zsombor},
pdftitle={},
pdfkeywords={},
pdfsubject={},
pdflang={Magyar}}
\begin{document}
\maketitle
Az esszém célja az, hogy megmutassam, hogy hogyan vezetett a spektroszkópia, és ezen belül az első főcsoport elemeinek a színképe a kvantummechanika kialakulásához, és ezzel megmagyarázzam azt, hogy a ``Modern fizika laboratórium'' kísérlete a látszólag ``ártatlan'' neve ellenére rendkívül nagy erejű bizonyíték a kvantummechanika szükségességéhez. Ezt a Simonyi~Károly: A~fizika~kultúrtörténete című könyvében felvázolt, a kvantummechanikához vezető utak ``Balmer-Rydberg-Ritz-Bohr'' ágának kibontásával kivánom elérni.
\section{Rövid összefoglaló a történelmi előzményekről}
\label{sec:org1ecbaae}
A spektroszkópia tudománya rendkívül termékenynek bizonyult a történelem során. Forradalmasította a kémia anyagok vizsgálatát, hasznos segédeszköze lett a csillagászatnak, elválaszthatatlan a kvantummechanika felfedezésétől, és máig egy rendkívül elterjedt vizsgálati módszer. A története rendkívül korán elkezdődött, már ókori római írók (pl. Plinius\cite[p.~438-439]{plinius}) is említettek tárgyakat, amelyekben ma a prizmát véljük felfedezni. A spektrum mai nevét közismerten Newton adta, majd utána Fraunhofer tett nagy előrelépéseket: felvette a Nap spektrumát és észrevett benne vonalakat amelyek hiányoznak, ő készítette el az első diffrakciós rácsot, továbbá ő fedezte fel, hogy a hevített anyagok csakis bizonyos hullámhossztartományokban bocsájtanak ki fényt. Tiszteletből róla neveztük el a Fraunhofer-diffrakciót.
Fraunhofert követően először prizmákat használtak a folytonos színkép felvételére, de a tudomány hamar áttért a rácsok használatára, mert azokkal pontosabb munkát lehetett végezni. A spektroszkópia népszerűségének a növekedése magával vonta azt is, hogy 19. században megnőtt az igény az egyre kisebb rácsállandójú diffrakciós rácsokra. Ebben a versenyben egy neves magyar fizikus is képviseltette magát: Jedlik~Ányos\cite{jedlik}\footnote{Erről először Groma István előadásán hallottam, ezt ellenőriztem le}. Jedlik a rácsai révén Európa-szerte híressé vált.
A hidrogén színképének a felvétele először Anders~Jones~Ångström nevéhez köthető, továbbá ő vetette fel először, hogy anyagok ugyanazokban a sávokban nyelvnek el fényt, mint amelyekben kibocsájtanak a vonalas spektrumban.
A hevített kémai elemek spektrumának a szisztematikus vizsgálatát Robert~Bunsen és Gustav~Kirchoff neveihez szoktuk kötni. Bunsent ma kémikusnak tartjuk, még Kirchoffot fizikusnak, a spektroszkópiát pedig a ``fizikai kémia'' címszó alatt találni a könyvtárakban. Úttörők voltak a spektroszkópia terén, és egy remek példa az interdiszciplinaritás hasznosságára. Igazolták Ångströmnek az emissziós és abszorpciós vonalakra vonatkozó állítását. Kirchoff az eredményeik alapján fogalmazta meg a spektroszkópia három törvényét:
\begin{itemize}
\item hevített szilárd anyagok és folyadok, illetve nagynyomású gázok folytonos spektrumban világítanak
\item alacsony nyomású gázok vonalas színképet produkálnak
\item ha egy folytonos színképet alacsony nyomású gázon átnézve szemlélünk, elnyelési vonalakat fogunk látni a folytonos spektrumban
\end{itemize}
Bunsen egy, az elődeinél jobb hevítőeszközt alkotott (a Bunsen-égőt), amelyet a kutatásaikban használtak is. Bunsen és Kirchoff voltak az első, akik felvették a lítium, nátrium és a kálium spektrumát, illetve a módszer segítségével felfedezték a céziumot és a rubídiumot. Ezen kívül bebizonyították, hogy a spektroszkópia alkalmas olyan esetekben is, amikor az adott elem csak nyomokban fedezhető fel a vizsgált anyagban. Az Kirchoff és Bunsen az eredményeiket két további kollégájukkal közösen publikálták.
Ők voltak az elsők abban is, hogy felvették egy csillag színképét az éjjeli égboltról (tehát nem a Nap színképét). Ezt ők maguk is annyira nagy előrelépésnek élték meg, hogy Bunsen a görög mitológiára hivatkozva így fogalmazott: ``megloptuk az isteneket''\cite{simonyi}.
\section{A kvantummechanikához vezető lépések}
\label{sec:orgdc5d95b}
Johann~Balmer volt az, akinek először feltűnt, hogy fel lehet írni a hidrogén színképében szereplő hullámhosszokat egy olyan képlettel, amelyben bizonyos paraméterek természetes számok -- nem pedig ``folytonos'', valós számok. A képlete a következő volt:
\begin{equation}
\lambda = B\left( \frac{n^2}{n^2-m^2}\right)
\end{equation}
ahol $\lambda$ a kibocsátott vagy elnyelt hullámhossz, $B$ egy empirikus konstans, $m=2$ és $n$ a term száma (a vonal sorszáma), ahol $n > m$. A képlete jóslatai nem csak arra voltak hasznosak, hogy a mérési eredmények pontatlanságára rámutassanak, de a képlet néhány, még nem felfedezett vonalat is megjósolt.
Sajnos ennek a képletnek a magyarázó ereje nem volt elegendő, azonban néhány lépéssel más alakra hozható. Johannes~Rydberg egy átrendezéssel olyan formára alakította a képletet, amelyet könnyebb volt általánosítani, illetve az értelmezése könnyebb -- a mai kvantummechanika eredményeivel majdnem triviális.\footnote{Az egyszerűség ténye szerintem jól szemlélteti, hogy mennyire fontos volt ez az eredmény a kvantummechanika kialakításában: ha nem ehhez "illesztették" volna a modellt, akkor előfordulhatna, hogy pl. perturbáció útján vezettük volna le, vagy valamely mellékfejezetben tárgyalnánk ma a témakört.}
\begin{equation}
\frac{1}{\lambda} = R\left( \frac{1}{n_1^2} - \frac{1}{n_2^2} \right)
\end{equation}
A képletében szereplő $R$ konstanst tiszteletből Rydberg-állandónak neveztük el, illetve $n_1$ és $n_2$ azok, amelyeket ma az elektronpályák főkvantumszámának nevezünk. $n_1$ a kiindulási héj száma, $n_2$ annak a héjnak a száma, ahova az elektron érkezik. A képlet egyszerűen kiterjeszthető olyan ionokra is, ahol 1 elektron van az atommag körül. Rydberg képlete jól modellezte a hidrogén és az alkálifémek spektrumát\footnote{Az alkálifémek spektrumának egy, modern igényeket kielégítő levezetéséhez már kellenek a relativisztikus korrekciók és a spin-pálya kölcsönhatás -- azaz a finomfelhasadás és a hiperfinom struktúra -- figyelembe vétele -- tehát a Dirac-egyenletből indulunk ki --, de ezek elsősorban a felhasad színképvonalak modellezéséhez szükségesek: maguk a dublettek a megjósolt helyek közelében helyezkednek el.\cite{függ}}, az általánosítását -- azaz korrekcióját a többi kémiai elemre -- Walther~Ritz készítette el (ez akkor még szintúgy egy tisztán empirikus képlet volt):
\begin{equation}
\tilde{\nu} = A - \frac{N}{(n + \alpha + \beta (A - \nu))^2}
\end{equation}
Itt $\tilde{\nu}$ a korrigált hullámszám, $\nu$ a hullámszám ($\frac{1}{\lambda}$), $A$ a sorozat végtelenben vett határértéke, $R$ a Rydberg-állandó, $\alpha$ és $\beta$ empirikus konstansok, $n$ pedig a term száma.\footnote{A forrásban történelmi okokból más betűket használt a szerző.}
Ezen korábbi eredmények alapján Niels~Bohr megalkotta az atommodelljét és megmutatta az első négy kémiai elem (H, He, Li, Be) esetében azt, hogy hogy kell a modellt alkamazni. Ezt modellt tanítják ma a középiskolákban. A Bohr-modell az előbb vázolt képleteket nem csupán magába foglalja, de korrekciókat is ad hozzá:
\begin{equation}
R_M = R_{\infty} \left(\frac{1}{1 + \frac{m_e}{M}} \right)
\end{equation}
ahol $m_e$ az elektron tömege, $M$ az atommag tömege, és $R_{\infty}$ a végtelen tömegű atommaghoz tartozó Rydberg-állandó:
\begin{equation}
R_{\infty} = \frac{1}{hc} \cdot \frac{m_ee^4}{2\hbar^2(4\pi\epsilon_0)^2}
\end{equation} és
\begin{equation}
\alpha = \frac{e^2}{\hbar c(4\pi\epsilon_0)}
\end{equation}
Ahol az $h$ a Planck-állandó $c$ a fénysebesség, $\epsilon_0$ a vákuum permittivitása és $e$ az elektron töltése. Az újonnan bevezetett $\alpha$-t a finomszerkezeti állandónak neveztük el.
Mint látható, a fentebbi képletek a klasszikus mechanikából már ismert redukált tömeg fogalmát alkalmazzák.
\section{Konklúzió}
\label{sec:org18f9e85}
A történelmi előzmények rámutatnak, hogy azért az első főcsoport elemeinek színképe volt úttörő jelentőségű, mert az egyszerű elektronszerkezetük miatt egyrészt a színkép struktúrája egyszerű lesz, másrészt a kvantummechanika egyenletei viszonylag könnyen jóslatokat produkálnak, amelyeket utána ellenőrizni tudunk.
Az új fizikai elméletek első igazoló kísérletei rendre ezt a két tulajdonságot ötvözik: ha az mérési hiba nagy lenne, illetve a kép zavaros, akkor bármilyen egyszerű is a modell, a kísérletben nehézkesen lesz használható. Ha azonban maga az rendszer lenne bonyolult, amit mérünk, akkor a még gyerekcipőben járó elmélet nincs eléggé kidolgozva ahhoz, hogy hamar igazolható jóslatokat tegyünk. A kísérlet megválasztása tehát a tudomány inkrementális-építkező (``óriások vállán állunk'') tudásképző modelljébe illeszkedik, és kulcsfontosságú például akkor, ha egy fizikus demonstrálni kívánja laikusok számára, hogy a kvantummechanika nem a semmiből jött vagy néhány zseni fejében pattant ki és köze sincs a valósághoz, illetve nem puszta értelmiségi hóbort, hanem egy történelmileg mélyen beágyazott, kísérletekkel jól alátámasztott, a magyarázatokhoz szükséges modell.
Szeretem idézni az alábbi ( \href{https://quoteinvestigator.com/2011/05/13/einstein-simple/}{alátámasztatlanul Einstein-nek tulajdonított} ) megfogalmazást Occam borotvapengéjéről: ``Everything should be made as simple as possible, but not simpler.'' A laboratóriumi kísérlet az új modellek szükségességét támasztja alá.
\begin{thebibliography}{9}
\bibitem{plinius}
The Natural History of Pliny,
\url{https://babel.hathitrust.org/cgi/pt?id=mdp.39015020434141&view=1up&seq=457}
\bibitem{Bunsen-Kirchoff}
Pierre Prevost, Balfour Stewart, Gustav Kirchhoff, Robert Bunsen,
The Laws of Radiation and Absorption
\url{https://books.google.com/books?id=3KwRAAAAYAAJ}
American Book Company,
1901
\bibitem{ritz}
Walther Ritz,
New Law of Series Spectra
\url{https://gallica.bnf.fr/ark:/12148/bpt6k33824.image.f185}
Astrophysical Journal, Vol XXVIII, No. 3,
October 1908
\bibitem{függ}
Kürti Jenő,
Modern fizika laboratórium, A függelék: Atomspektroszkópia
\url{http://wigner.elte.hu/koltai/labor/parts/modernA.pdf}
\bibitem{simonyi}
Simonyi Károly,
A fizika kultúrtörténete
\url{https://mersz.hu/hivatkozas/m432afk_281}
Akadémiai Kiadó,
2021
\bibitem{jedlik}
\url{https://web.archive.org/web/20100819022858/http://jedliktarsasag.hu/mszh/jedlik.html}
\end{thebibliography}
\end{document}

View File

@ -1,840 +0,0 @@
0 0.0
1 0.0
2 0.0
3 0.0
4 0.0
5 0.0
6 0.0
7 0.0
8 0.0
9 0.0
10 0.0
11 0.0
12 0.0
13 0.0
14 0.0
15 0.0
16 0.0
17 0.0
18 0.0
19 0.0
20 0.0
21 0.0
22 0.0
23 0.0
24 0.0
25 0.0
26 0.0
27 0.0
28 0.0
29 0.0
30 0.0
31 0.0
32 0.0
33 0.0
34 0.0
35 0.0
36 0.0
37 0.0
38 0.0
39 0.0
40 0.0
41 0.0
42 0.0
43 0.0
44 0.0
45 0.0
46 0.0
47 0.0
48 0.0
49 0.0
50 0.0
51 0.0
52 0.0
53 0.0
54 0.0
55 0.0
56 0.0
57 0.0
58 0.0
59 0.0
60 0.0
61 0.0
62 0.0
63 0.0
64 0.0
65 0.0
66 0.0
67 0.0
68 0.0
69 0.0
70 0.0
71 0.0
72 0.0
73 0.0
74 0.0
75 0.0
76 0.0
77 0.0
78 0.0
79 0.0
80 0.0
81 0.0
82 0.0
83 0.0
84 0.0
85 0.0
86 0.0
87 0.0
88 0.0
89 0.0
90 0.0
91 0.0
92 0.0
93 0.0
94 0.0
95 0.0
96 0.0
97 0.0
98 0.0
99 0.0
100 0.0
101 0.0
102 0.0
103 0.0
104 0.0
105 0.0
106 0.0
107 0.0
108 0.0
109 0.0
110 0.0
111 0.0
112 0.0
113 0.0
114 0.0
115 0.0
116 0.0
117 0.0
118 0.0
119 0.0
120 0.0
121 0.0
122 0.0
123 0.0
124 0.0
125 0.0
126 0.0
127 0.0
128 0.0
129 0.0
130 0.0
131 0.0
132 0.0
133 0.0
134 0.0
135 0.0
136 0.0
137 0.0
138 0.0
139 0.0
140 0.0
141 0.0
142 0.0
143 0.0
144 0.0
145 0.0
146 0.0
147 0.0
148 0.0
149 0.0
150 0.0
151 0.0
152 0.0
153 0.0
154 0.0
155 0.0
156 0.0
157 0.0
158 0.0
159 0.0
160 0.0
161 0.0
162 0.0
163 0.0
164 0.0
165 0.0
166 0.0
167 0.0
168 0.0
169 0.0
170 0.0
171 0.0
172 0.0
173 0.0
174 0.0
175 0.0
176 0.0
177 0.0
178 0.0
179 0.0
180 0.0
181 0.0
182 0.0
183 0.0
184 0.0
185 0.0
186 0.0
187 0.0
188 0.0
189 0.0
190 0.0
191 0.0
192 0.0
193 0.0
194 0.0
195 0.0
196 0.0
197 0.0
198 0.0
199 0.0
200 0.0
201 0.0
202 0.0
203 0.0
204 0.0
205 0.0
206 0.0
207 0.0
208 0.0
209 0.0
210 0.0
211 0.0
212 0.0
213 0.0
214 0.0
215 0.0
216 0.0
217 0.0
218 0.0
219 0.0
220 0.0
221 0.0
222 0.0
223 0.0
224 0.0
225 0.0
226 0.0
227 0.0
228 0.0
229 0.0
230 0.0
231 0.0
232 0.0
233 0.0
234 0.0
235 0.0
236 0.0
237 0.0
238 0.0
239 0.0
240 0.0
241 0.0
242 0.0
243 0.0
244 0.0
245 0.0
246 0.0
247 0.0
248 0.0
249 2.0
250 0.0
251 0.0
252 4.0
253 0.0
254 0.0
255 45.0
256 0.0
257 0.0
258 0.0
259 75.0
260 0.0
261 0.0
262 511.0
263 0.0
264 0.0
265 135.0
266 0.0
267 0.0
268 0.0
269 115.0
270 0.0
271 0.0
272 193.0
273 0.0
274 0.0
275 93.0
276 0.0
277 0.0
278 177.0
279 0.0
280 0.0
281 0.0
282 75.0
283 0.0
284 0.0
285 178.0
286 0.0
287 0.0
288 98.0
289 0.0
290 0.0
291 0.0
292 87.0
293 0.0
294 0.0
295 172.0
296 0.0
297 0.0
298 93.0
299 0.0
300 0.0
301 174.0
302 0.0
303 0.0
304 0.0
305 107.0
306 0.0
307 0.0
308 221.0
309 0.0
310 0.0
311 133.0
312 0.0
313 0.0
314 0.0
315 137.0
316 0.0
317 0.0
318 271.0
319 0.0
320 0.0
321 136.0
322 0.0
323 0.0
324 290.0
325 0.0
326 0.0
327 0.0
328 168.0
329 0.0
330 0.0
331 463.0
332 0.0
333 0.0
334 204.0
335 0.0
336 0.0
337 234.0
338 0.0
339 0.0
340 0.0
341 467.0
342 0.0
343 0.0
344 221.0
345 0.0
346 0.0
347 470.0
348 0.0
349 0.0
350 0.0
351 250.0
352 0.0
353 0.0
354 486.0
355 0.0
356 0.0
357 255.0
358 0.0
359 0.0
360 216.0
361 0.0
362 0.0
363 0.0
364 497.0
365 0.0
366 0.0
367 257.0
368 0.0
369 0.0
370 540.0
371 0.0
372 0.0
373 0.0
374 287.0
375 0.0
376 0.0
377 271.0
378 0.0
379 0.0
380 567.0
381 0.0
382 0.0
383 314.0
384 0.0
385 0.0
386 0.0
387 663.0
388 0.0
389 0.0
390 314.0
391 0.0
392 0.0
393 686.0
394 0.0
395 0.0
396 0.0
397 375.0
398 0.0
399 0.0
400 412.0
401 0.0
402 0.0
403 786.0
404 0.0
405 0.0
406 396.0
407 0.0
408 0.0
409 0.0
410 716.0
411 0.0
412 0.0
413 380.0
414 0.0
415 0.0
416 730.0
417 0.0
418 0.0
419 0.0
420 348.0
421 0.0
422 0.0
423 368.0
424 0.0
425 0.0
426 746.0
427 0.0
428 0.0
429 367.0
430 0.0
431 0.0
432 0.0
433 692.0
434 0.0
435 0.0
436 357.0
437 0.0
438 0.0
439 803.0
440 0.0
441 0.0
442 356.0
443 0.0
444 0.0
445 0.0
446 387.0
447 0.0
448 0.0
449 766.0
450 0.0
451 0.0
452 458.0
453 0.0
454 0.0
455 0.0
456 790.0
457 0.0
458 0.0
459 448.0
460 0.0
461 0.0
462 876.0
463 0.0
464 0.0
465 440.0
466 0.0
467 0.0
468 0.0
469 443.0
470 0.0
471 0.0
472 928.0
473 0.0
474 0.0
475 510.0
476 0.0
477 0.0
478 0.0
479 1045.0
480 0.0
481 0.0
482 568.0
483 0.0
484 0.0
485 1096.0
486 0.0
487 0.0
488 550.0
489 0.0
490 0.0
491 0.0
492 581.0
493 0.0
494 0.0
495 1166.0
496 0.0
497 0.0
498 634.0
499 0.0
500 0.0
501 0.0
502 1260.0
503 0.0
504 0.0
505 656.0
506 0.0
507 0.0
508 1320.0
509 0.0
510 0.0
511 716.0
512 0.0
513 0.0
514 0.0
515 651.0
516 0.0
517 0.0
518 1325.0
519 0.0
520 0.0
521 686.0
522 0.0
523 0.0
524 0.0
525 1374.0
526 0.0
527 0.0
528 652.0
529 0.0
530 0.0
531 1483.0
532 0.0
533 0.0
534 715.0
535 0.0
536 0.0
537 0.0
538 720.0
539 0.0
540 0.0
541 1480.0
542 0.0
543 0.0
544 773.0
545 0.0
546 0.0
547 1662.0
548 0.0
549 0.0
550 0.0
551 830.0
552 0.0
553 0.0
554 1808.0
555 0.0
556 0.0
557 861.0
558 0.0
559 0.0
560 0.0
561 886.0
562 0.0
563 0.0
564 1719.0
565 0.0
566 0.0
567 848.0
568 0.0
569 0.0
570 1630.0
571 0.0
572 0.0
573 0.0
574 822.0
575 0.0
576 0.0
577 1586.0
578 0.0
579 0.0
580 802.0
581 0.0
582 0.0
583 0.0
584 781.0
585 0.0
586 0.0
587 1567.0
588 0.0
589 0.0
590 888.0
591 0.0
592 0.0
593 1604.0
594 0.0
595 0.0
596 0.0
597 831.0
598 0.0
599 0.0
600 1570.0
601 0.0
602 0.0
603 738.0
604 0.0
605 0.0
606 0.0
607 855.0
608 0.0
609 0.0
610 1561.0
611 0.0
612 0.0
613 815.0
614 0.0
615 0.0
616 1706.0
617 0.0
618 0.0
619 0.0
620 830.0
621 0.0
622 0.0
623 841.0
624 0.0
625 0.0
626 1645.0
627 0.0
628 0.0
629 0.0
630 871.0
631 0.0
632 0.0
633 1723.0
634 0.0
635 0.0
636 842.0
637 0.0
638 0.0
639 1816.0
640 0.0
641 0.0
642 0.0
643 967.0
644 0.0
645 0.0
646 1011.0
647 0.0
648 0.0
649 2087.0
650 0.0
651 0.0
652 973.0
653 0.0
654 0.0
655 0.0
656 2019.0
657 0.0
658 0.0
659 969.0
660 0.0
661 0.0
662 2009.0
663 0.0
664 0.0
665 0.0
666 1003.0
667 0.0
668 0.0
669 994.0
670 0.0
671 0.0
672 2011.0
673 0.0
674 0.0
675 976.0
676 0.0
677 0.0
678 0.0
679 1956.0
680 0.0
681 0.0
682 1055.0
683 0.0
684 0.0
685 2113.0
686 0.0
687 0.0
688 0.0
689 1050.0
690 0.0
691 0.0
692 1110.0
693 0.0
694 0.0
695 2237.0
696 0.0
697 0.0
698 1160.0
699 0.0
700 0.0
701 0.0
702 2286.0
703 0.0
704 0.0
705 1137.0
706 0.0
707 0.0
708 2359.0
709 0.0
710 0.0
711 0.0
712 1246.0
713 0.0
714 0.0
715 1227.0
716 0.0
717 0.0
718 2523.0
719 0.0
720 0.0
721 1248.0
722 0.0
723 0.0
724 0.0
725 2604.0
726 0.0
727 0.0
728 1258.0
729 0.0
730 0.0
731 2761.0
732 0.0
733 0.0
734 0.0
735 1403.0
736 0.0
737 0.0
738 1380.0
739 0.0
740 0.0
741 2885.0
742 0.0
743 0.0
744 1521.0
745 0.0
746 0.0
747 0.0
748 3055.0
749 0.0
750 0.0
751 1642.0
752 0.0
753 0.0
754 3287.0
755 0.0
756 0.0
757 1696.0
758 0.0
759 0.0
760 0.0
761 1774.0
762 0.0
763 0.0
764 3447.0
765 0.0
766 0.0
767 1866.0
768 0.0
769 0.0
770 0.0
771 3657.0
772 0.0
773 0.0
774 1957.0
775 0.0
776 0.0
777 4087.0
778 0.0
779 0.0
780 2053.0
781 0.0
782 0.0
783 0.0
784 2089.0
785 0.0
786 0.0
787 4214.0
788 0.0
789 0.0
790 2256.0
791 0.0
792 0.0
793 0.0
794 4403.0
795 0.0
796 0.0
797 2247.0
798 0.0
799 0.0
800 4514.0
801 0.0
802 0.0
803 2348.0
804 0.0
805 0.0
806 0.0
807 2427.0
808 0.0
809 0.0
810 5158.0
811 0.0
812 0.0
813 2796.0
814 0.0
815 0.0
816 0.0
817 5996.0
818 0.0
819 0.0
820 3432.0
821 0.0
822 0.0
823 8402.0
824 0.0
825 0.0
826 5198.0
827 0.0
828 0.0
829 0.0
830 6574.0
831 0.0
832 0.0
833 22671.0
834 0.0
835 0.0
836 416677.0
837 0.0
838 0.0
839 0.0

Binary file not shown.

After

Width:  |  Height:  |  Size: 311 KiB

BIN
modfizlab/ernyo.png Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 290 KiB

BIN
modfizlab/histogram.png Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 30 KiB

View File

@ -26,17 +26,6 @@ Mi a mérésünkben egy efféle elrendezésben vizsgáltuk a "kvantumradírt". A
Egy [[https://en.wikipedia.org/wiki/Mach%E2%80%93Zehnder_interferometer][Mach-Zehnder interferométerrel]] dolgoztunk. Ennek minkét nyalábútvonalába egy-egy polárszűrőt tettünk, és ezekkel a fénynyalábokat különbözőképp polarizáltuk. A megvalósításunk virtuális volt, amelyet a laborvezető bocsájtott rendekezésünkre. A fényképezést egy screenshot készítésével tudtuk helyettesíteni.
A mérés elvégzéséhez a laborvezető egy szoftvert adott rendelkezésünkre, mely a következőképp nézett ki:
\begin{figure}[h]
\centering
\includegraphics[width=0.9\linewidth]{sw.png}
\caption{A szoftver}
\label{fig:test1}
\end{figure}
Itt találhatók angol nyelven felcímkézett gombok, amelyekkel válhatunk a klasszikus és a kvantumos eset közt, engedélyezhetjük és állíthatjuk a polárszűrőket, illetve leolvashatjuk az ún. distinguishability paramétert ( erről később beszámolunk ).
* Számolási feladatok
** Kvantumradír-kísérlet a Mach-Zehnder interferométerben
@ -206,9 +195,11 @@ I\sim 1+\cos(2\phi)\cos\alpha
\end{equation}
adódik, tehát a kontraszt most is
\begin{equation}
V_2=cos (2\phi)
V_2=\cos (2\phi).
\end{equation}
* Mérési elrendezés
* Mérési adatok
** 1. Jelölt és jelöletlen útvonal
@ -225,8 +216,8 @@ V_2=cos (2\phi)
\caption{Jelölt útvonal}
\label{fig:test2}
\end{figure}
\pagebreak
\pagebreak
** 2. Kvantumradírozás szemléltetése
\begin{figure}[h]
@ -235,7 +226,6 @@ V_2=cos (2\phi)
\caption{Bekapcsolt radírozó polárszűrő}
\end{figure}
\pagebreak
** 3. Az egyik ernyőn nincs jel
\begin{figure}[h]
@ -253,88 +243,36 @@ V_2=cos (2\phi)
\end{figure}
\pagebreak
** Visibility paraméter és distinguishability összefüggése
** 4. mérési feladat
|----------+---------+---------+-------|
| $\alpha$ | $\beta$ | $D_1$ | $V_1$ |
|----------+---------+---------+-------|
| 75 | 120 | 0.81650 | TBD |
|----------+---------+---------+-------|
\begin{figure}[h]
\centering
\includegraphics[width=.9\linewidth]{ernyo.png}
\caption{Az ernyő képe}
\label{fig:ernyo1}
\end{figure}
A számításhoz először a képernyőképet GIMP segítśegével dolgoztam fel, egy perspektívatranszformációval négyzetbe transzformáltam az ernyő képét. Ezután szürkeárnyalatossát tettem a képet, és körülvágtam.
\begin{figure}[h]
\centering
\includegraphics[width=.9\linewidth]{ernyo-transformed.png}
\caption{Az ernyő képe a tranformációk után}
\label{fig:ernyo2}
\end{figure}
** 5. mérési feladat
Ezután betöltöttem a képet a python nyelv parancsértelmezőjébe, negatívba tettem, és hisztogramot csináltam 15-ös thresholddal (Ha a pixel fényessége 15 felett volt, akkor a az oszlop megfelelő gyűjtőjében a számlálót növeltem.)
|----------+---------+---------|
| $\alpha$ | $\beta$ | $D_1$ |
|----------+---------+---------|
| 45 | 135 | 1.0 |
| 50 | 130 | 0.98481 |
| 55 | 125 | 0.93969 |
| 60 | 120 | 0.86603 |
| 65 | 115 | 0.76604 |
| 70 | 110 | 0.64279 |
| 75 | 105 | 0.50000 |
| 80 | 100 | 0.34202 |
| 85 | 95 | 0.17365 |
| 90.5 | 90.5 | 0.000 |
|----------+---------+---------|
\begin{figure}[h]
\centering
\includegraphics[width=.9\linewidth]{histogram.png}
\caption{A hisztogram, csakis szemléltetésnek használtuk}
\label{fig:ernyo1}
\end{figure}
A hisztogram adatai közt megkerestem a nekünk kellő minimum-és maximumértékeket, amelyek $I_{min}=155\pm12.5$, illetve $I_{max}=510\pm22.6$-nak adódtak (a hibát a poisson-eloszlás szórásával, $\sqrt{\mu}$-vel becsültem).
\begin{equation}
V_1=\frac{I_{max}-I_{min}}{I_{max}+I_{min}}=0.534\pm0.045
\end{equation}
A hiba számítására a következő képleteket használtam:
\begin{equation}
\frac{\Delta V_1}{\Delta I_{max}}=|
\frac{1}{I_{max}+I_{min}}-\frac{I_{max}}{(I_{max}+I_{min})^2}+\frac{I_{min}}{(I_{max}+I_{min})^2}|
\end{equation}
\begin{equation}
\frac {\Delta V_1}{\Delta I_{min}} = |
-\frac{1}{I_{max}+I_{min}}-\frac{I_{max}}{(I_{max}+I_{min})^2}+\frac{I_{min}}{(I_{max}+I_{min})^2}|
\end{equation}
|----------+---------+---------+-----------------|
| $\alpha$ | $\beta$ | $D_1$ | $V_1$ |
|----------+---------+---------+-----------------|
| 75 | 120 | 0.81650 | $0.534\pm0.045$ |
|----------+---------+---------+-----------------|
Ennyire kevés adat alapján az adat alapján nem még nem megállapítható az összefüggés, bár annyi már sejthető, hogy a kettő négyzetösszege mindig 1 körül lesz.
\pagebreak
** Kontraszt paraméterek kvantumradír nélkül
#+CAPTION: /hmmm/
|---------+---------|
| $2\phi$ | $D_1$ |
|---------+---------|
| 90 | 1.0 |
| 80 | 0.98481 |
| 70 | 0.93969 |
| 60 | 0.86603 |
| 50 | 0.76604 |
| 40 | 0.64279 |
| 30 | 0.50000 |
| 20 | 0.34202 |
| 10 | 0.17365 |
| 0 | 0.000 |
|---------+---------|
A táblázatból is jól látszik, hogy a $D_1=cos(\pi/2-2\phi)$, azaz 90°-kal eltolt fázisú a visibility paraméterhez képest, amit a számolás eredményeképp kaptunk.
Minthogy ez egy virtuális labor által adott eredmény és szemmel is látható az elméleti eredmény helyessége, továbbá a hibaszámítás ill. az illesztés hibája gyakorlatilag 0 kell legyen ( leolvasási hiba egy kerekítés eredménye ($\pm5*10^{-6}$), illetve a lebegőpontos számok pontatlansága jelen esetben elhanyagolható), jelen és követkeő esetekben eltekintettem a koszinuszfüggvény illetszésétől.
Megj.: D1-re a szoftver hibát jelez, ha mindkét polarizátor 90 fokra van állítva.
\pagebreak
** Kontraszt paraméterek kvantumradírral
Megj.: D1-re a szoftver ***-ot jelez, ha mindkét polarizátor 90 fokra van állítva
** 6. mérési feladat
|----------+---------|
| $\alpha$ | $D_1$ |
| $\gamma$ | $D_1$ |
|----------+---------|
| 45 | 1.00000 |
| 54 | 0.95106 |
@ -348,10 +286,3 @@ Megj.: D1-re a szoftver hibát jelez, ha mindkét polarizátor 90 fokra van áll
| 126 | 0.95106 |
| 135 | 1.00000 |
|----------+---------|
Ismét igazolódott az a számolás: a periódus 90°
\pagebreak
* Diszkusszió
A kísérlet során egy szimulációban meg tudtuk figyelni kétréses kísérletben az interferenciakép eltűnését útvonaljelölés esetén, a kvantumradírozást, megtaláltuk a forrás polarizációs szögét, megtaláltuk a szoftverben levő distiguishability paraméter és a visibility paraméter közti összefüggést, illetve kimuttatuk szögfüggéses számolások helyességét. Sajnos a $D_1$ és $V_1$ közti összfüggést nem sikerült elsőre igazolni, azonban további mérésekkel endönthetőnek bizonyult a probléma.

333
modfizlab/lab15.tex Normal file
View File

@ -0,0 +1,333 @@
% Created 2021-03-09 Tue 21:58
% Intended LaTeX compiler: pdflatex
\documentclass[12pt,a4paper]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{graphicx}
\usepackage{grffile}
\usepackage{longtable}
\usepackage{wrapfig}
\usepackage{rotating}
\usepackage[normalem]{ulem}
\usepackage{amsmath}
\usepackage{textcomp}
\usepackage{amssymb}
\usepackage{capt-of}
\usepackage{hyperref}
\usepackage[magyar, english]{babel}
\hypersetup{ colorlinks=true,linkcolor=blue, linktoc=all, filecolor=magenta, urlcolor=cyan, pdfstartview=FitB,}
\usepackage{multirow}
\usepackage{braket}
\urlstyle{same}
\author{Rédl Anna, Barna Zsombor}
\date{Mérés időpontja: 2021. 02. 23.}
\title{15. Kvantumradír}
\hypersetup{
pdfauthor={Rédl Anna, Barna Zsombor},
pdftitle={15. Kvantumradír},
pdfkeywords={},
pdfsubject={},
pdfcreator={Emacs 27.1 (Org mode 9.1.8)},
pdflang={English}}
\begin{document}
\maketitle
\begin{titlepage}
\pagenumbering{gobble}
\maketitle
\end{titlepage}
\pagebreak
\section{Bevezetés}
\label{sec:org5341b67}
Ugyan Richard Feynman szerint a kétréses interferenciakísérlet magyarázata a kvantummechanika egyetlen rejtélye, azonban a jelenség bizonyos esetei klasszikus keretek közt is értelmezhetők.
Mi a mérésünkben egy efféle elrendezésben vizsgáltuk a "kvantumradírt". A kvantumradírnak nevezzük azt a jelenséget ami akkor lép fel, ha egy interferométerben az ún. "útvonaljelölés" ( a mi esetünkben a megkülönböztető tulajdonság a fotonok polarizációja ) eltünteti az interferenciát, de még a mielőtt a detektorba érnének a részecskesugár, megszüntetjük ("kiradírozzuk") az útvonaljelölést -- azaz kitöröljük a korábban szerzett útvonaljelölést, és az interferenciakép helyreáll.
\section{Mérési összeállítás}
\label{sec:orgc577efa}
Egy \href{https://en.wikipedia.org/wiki/Mach\%E2\%80\%93Zehnder\_interferometer}{Mach-Zehnder-interferométerrel} dolgozunk. Ennek minkét nyalábútvonalába egy-egy polárszűrőt tettünk, és ezekkel a fénynyalábokat különbözőképp polarizáltuk. A megvalósításunk virtuális volt, amelyet a laborvezető bocsájtott rendekezésünkre. A fényképezést egy screenshot készítésével tudtuk helyettesíteni.
\section{Számolási feladatok}
\label{sec:org0765e8b}
\begin{enumerate}
\item A lézernyalábok által bezárt szög számítása Mach-Zender elrendezésben
\label{sec:org393a8fe}
Ha a négy tükör nem tökéletesen párhuzamos, a lencsén áthaladva a nyalábok a fókuszsíkban nem egyetlen pontban metszik egymást, hanem kettő pontban ( \(S'\) és \(S''\)).
Ha feltételezzük, hogy a nyalábok \(\alpha\) szöge nagyon kicsi, akkor az \(S'\) és \(S''\) pontok, valamint a lencse fősíkjának az \(S'S''\) szakasz felezőjével szemközti pontja által definiált egyenlő szárú háromszög szárai \(\approx f\), ahol \(f\) a lencse fókusztávolsága.
Ebből a koszinusztétel alapján az \(S'S''\) szakasz \(d\) hosszúsága az alábbiak szerint becsülhető:
\begin{equation}
d^2\approx 2f^2(1-\cos\alpha)
\end{equation}
\(\cos{\alpha}\) -t másodrendig sorbafejtve:
\begin{equation*}
\cos{\alpha}\approx 1-\frac{\alpha^2}{2}+\mathcal{O}(\alpha^4)
\end{equation*}
Ezt az előző kifejezésbe behelyettesítve:
\begin{equation*}
d=f\alpha
\end{equation*}
A \(d\) szakasz hosszára kaphatunk egy másik egyenletet is az alaphán, hogy az \(S'\) és \(S''\) pontokból gömbhullámok indulnak ki. A \(\theta_m\) irányok, az \(m\) rendek, a \(d\) réstávolság és a \(\lambda\) hullámhossz között az alábbi összefüggés írható fel:
\begin{equation}
d\sin\theta_m=m\lambda
\end{equation}
Ha \(\theta_m\) kicsi és a lencse és az ernyő \(L\) távolságára igaz, hogy \(L\gg f\), akkor élhetünk az alábbi közelítéssekkel a \(\theta_m\) irányra és az adott interferenciacsíknak a középső csíktól való \(l_m\) távolságára vonatkozóan:
\begin{equation*}
\theta \approx \sin\theta\approx\tan\theta\approx\frac{l_m}{L}
\end{equation*}
Ez alapján az egyenlet a következő formát ölti:
\begin{equation*}
d \approx \frac{m\lambda L}{l_m}
\end{equation*}
Ebbe behelyettesítve a \$d\$-re kapott előző kifejezést, a nyalábok párhuzamosságának \(\alpha\) hibájára a következő kifejezés adódik:
\begin{equation}\label{eq: alpha}
\alpha\approx \frac{m\lambda L}{fl_m}
\end{equation}
\item Láthatósági paraméter
\label{sec:orga23221c}
\begin{enumerate}
\item Első eset
\label{sec:org6c6787f}
Vegyünk két, lineárisan polarizált állapotot \(2\cdot\phi\) szögeltéréssel és \(\alpha\) fáziskülönbséggel. Ezek Jones-vektorai a következők:
\begin{align}
\mathbf{E}_1=
\left(
\begin{array}{c}
\cos\phi \\
\sin\phi
\end{array}
\right)
&\hspace{24pt}
\mathbf{E}_2=
\left(
\begin{array}{c}
\cos\phi \\
-\sin\phi
\end{array}
\right)
\cdot
e^{i\alpha}
\end{align}
A szuperpozíció elve szerint összeadhatjuk a fentieket, hogy az eredő térerősséget megkapjuk:
\begin{equation}
\mathbf{E} =\left(
\begin{array}{c}
\cos\phi(1+e^{i\alpha}) \\
\sin\phi(1-e^{i\alpha})
\end{array}
\right)
\end{equation}
Az intenzitás a térerősségnégyzet abszolút értékével arányos, ezért felírhatjuk így is:
\begin{equation}
I\sim \mathbf{E}\mathbf{E}^{*}=\left(
\begin{array}{c c}
\cos\phi(1+e^{-i\alpha}); &
\sin\phi (1-e^{-i\alpha})
\end{array}
\right)\cdot
\left(
\begin{array}{c}
\cos\phi(1+e^{i\alpha}) \\
\sin\phi (1-e^{i\alpha})
\end{array}
\right)
\end{equation}
A beszorzás után a következő marad:
\begin{equation}
I\sim 1+\cos(2\phi)\cos\alpha
\end{equation}
Ha ezt behelyzettesítjük a kontraszt paraméter definíciójába, akkor a következőt kapjuk:
\begin{equation}
V_1=\frac{I_\text{max}-I_\text{min}}{I_\text{max}+I_\text{min}}=\cos(2\phi),
\end{equation}
Ez a fenti eredmény lesz a kontraszt paraméter képlete az első elrendezésünkben.
\item Második eset
\label{sec:org8519298}
A két polarizátort ±45°-os állásba tesszük és a harmadik polárszűrőt (a kvantumradírt) \(\phi\) szögű irányba állítjuk. A ±45°-os állású Jones-vektorok:
\begin{align}
\mathbf{E}_1 =\frac{1}{\sqrt{2}}\left(
\begin{array}{c}
1 \\
1
\end{array}
\right)&\hspace{24pt}
\mathbf{E}_2 =\frac{1}{\sqrt{2}}\left(
\begin{array}{c}
1 \\
-1
\end{array}
\right)e^{i\alpha}.
\end{align}
A polarizátor Jones-mátrixa pedig legyen egy általános polarizátor mátrixa:
\begin{equation}
\mathbf{P}=\left(\begin{array}{c c}
\cos^2\phi& \sin\phi\cos\phi \\
\sin\phi\cos\phi & \sin^2\phi
\end{array}
\right).
\end{equation}
Az eredő térerősséget ismét szuperpozícióból kapjuk:
\begin{equation}
\mathbf{E}=\mathbf{P}(\mathbf{E}_1+\mathbf{E}_2)=\frac{1}{\sqrt{2}}(\cos\phi (1+e^{i\alpha})+\sin\phi (1-e^{i\alpha}))
\cdot
\left(
\begin{array}{c}
\cos\phi \\
\sin\phi
\end{array}
\right).
\end{equation}
Ennek kiszámítva az abszolút értékének négyzetét, az intenzitásra ebben az esetben is
\begin{equation}
I\sim 1+\cos(2\phi)\cos\alpha
\end{equation}
adódik, tehát a kontraszt most is
\begin{equation}
V_2=\cos (2\phi).
\end{equation}
\end{enumerate}
\end{enumerate}
\section{Mérési elrendezés}
\label{sec:orgbb2e722}
\section{Mérési adatok}
\label{sec:org4965fd4}
\begin{enumerate}
\item 1. Jelölt és jelöletlen útvonal
\label{sec:orge785d7c}
\begin{figure}[h]
\centering
\includegraphics[width=0.9\linewidth]{1.png}
\caption{Jelöletlen útvonal}
\label{fig:test1}
\end{figure}
\begin{figure}[h]
\centering
\includegraphics[width=0.9\linewidth]{2.png}
\caption{Jelölt útvonal}
\label{fig:test2}
\end{figure}
\pagebreak
\item 2. Kvantumradírozás szemléltetése
\label{sec:org17ab964}
\begin{figure}[h]
\centering
\includegraphics[width=8cm]{3.png}
\caption{Bekapcsolt radírozó polárszűrő}
\end{figure}
\pagebreak
\item 3. Az egyik ernyőn nincs jel
\label{sec:org5a73f61}
\begin{figure}[h]
\centering
\includegraphics[width=.9\linewidth]{4.png}
\caption{Másik ernyőn nincs interferencia}
\label{fig:test3}
\end{figure}
\begin{figure}[h]
\centering
\includegraphics[width=.9\linewidth]{5.png}
\caption{Másik ernyőn van interferencia}
\label{fig:test4}
\end{figure}
\pagebreak
\item 4. mérési feladat
\label{sec:org3ac45a7}
\begin{center}
\begin{tabular}{rrrl}
\hline
\(\alpha\) & \(\beta\) & \(D_1\) & \(V_1\)\\
\hline
75 & 120 & 0.81650 & TBD\\
\hline
\end{tabular}
\end{center}
\item 5. mérési feladat
\label{sec:org5d37545}
\begin{center}
\begin{tabular}{rrr}
\hline
\(\alpha\) & \(\beta\) & \(D_1\)\\
\hline
45 & 135 & 1.0\\
50 & 130 & 0.98481\\
55 & 125 & 0.93969\\
60 & 120 & 0.86603\\
65 & 115 & 0.76604\\
70 & 110 & 0.64279\\
75 & 105 & 0.50000\\
80 & 100 & 0.34202\\
85 & 95 & 0.17365\\
90.5 & 90.5 & 0.000\\
\hline
\end{tabular}
\end{center}
Megj.: D1-re a szoftver \textbf{*}-ot jelez, ha mindkét polarizátor 90 fokra van állítva
\item 6. mérési feladat
\label{sec:org14d32b0}
\begin{center}
\begin{tabular}{rr}
\hline
\(\gamma\) & \(D_1\)\\
\hline
45 & 1.00000\\
54 & 0.95106\\
63 & 0.80902\\
72 & 0.58779\\
81 & 0.30902\\
90 & 0.00000\\
99 & 0.30902\\
108 & 0.58779\\
117 & 0.80902\\
126 & 0.95106\\
135 & 1.00000\\
\hline
\end{tabular}
\end{center}
\end{enumerate}
\end{document}

View File

@ -1,248 +1,28 @@
#+OPTIONS: tex:t
#+AUTHOR: Rédl Anna, Barna Zsombor
#+DATE: Mérés időpontja: 2021. 03. 09.
#+TITLE: 13. Molekulamodellezés
#+LATEX_CLASS: article
#+LATEX_CLASS_OPTIONS: [12pt,a4paper]
#+LATEX_HEADER: \usepackage[utf8]{inputenc}
#+LATEX_HEADER: \usepackage[magyar]{babel}
#+LATEX_HEADER: \hypersetup{ colorlinks=true,linkcolor=blue, linktoc=all, filecolor=magenta, urlcolor=cyan, pdfstartview=FitB,}
#+LATEX_HEADER: \usepackage{multirow}
#+LATEX_HEADER: \usepackage{braket}
#+LATEX_HEADER: \urlstyle{same}
#+OPTIONS: toc:nil
\newpage
* A mérés célja
A molekulamodellezés kísérlet célja többatomos rendszerek (molekulák) modellezése. Ez a feladat nem magától érthetődő, mivel a Schrödinger-egyenlet egzaktul a legtöbb molekulára nem oldható meg. Ezért molekulák modellezéséhez különböző közelítések használatára van szűkség. A mérés során Hartee-Fock közelítést használunk.
\newpage
* A mérési összeállítás és a mérés menete
A mérés során a molekulákat az Avogadro programmal modelleztük. Elsőnek kétatomos molekulákat ($H_2$, $N_2$, $O_2$ és $F_2$) modelleztünk, szingulett és triplett spinállapotban, és dokumentáltuk a kötési távolságot, és az optimális energiát.
Második feladatként a $H_2$ molekula elektronszerkezetét vizsgáltuk. Kirajzoltuk a HOMO és az első négy LUMO pályát, hogy aztán a LCAO megfontolások alapján elemezzük azokat.
Harmadik feladatként a lineáris víz problémáját vizsgáltuk rezgési analízis használatával. Az Avogadro programban a kötésszög rögzítésével hoztunk létre lineáris vízmolekulát, és feljegyeztük a rezgési frekvenciákat és a rezgési módusok infraaktivitását. Elkészítettük a valódi vízmolekulát is, és összehasonlítottuk energiáját a lineáris vízmolekuláéval.
Negyedik feladatként az izotópeffektust vizsgáltuk a benzol és a deuterizált benzol infravörös spektrumának kiszámolásával. Ehhez az előző feladathoz hasonló módon kimentettük a rezgési frekvenciákat és a rezgési módusok infraaktivitását.
Utolsó feladatként választottunk további két $C_6H_6$ molekulát a ChemSpider adatbázisból (divubykacetylene, dipropargyl), majd elmentettük mindkettőnek a keletkezési energiáját összehasonlítás céljából.
\newpage
* Mérési adatok
A mérés során végig B3LYP módszert használtuk $6-311+G$ bázison. A továbbiakban az $E$ mindenhol a keletkezési energiát jelöli
\newpage
** Elemmolekulák
| anyag | E[kJ/mol] | kötéshossz[Å] |
|-----------------+-------------+---------------|
| $H_2$ szinglett | -3088.294 | 0.743 |
| $H_2$ triplett | -2628.693 | 5.801 |
| $N_2$ szinglett | -287747.996 | 1.106 |
| $N_2$ triplett | -287046.505 | 1.284 |
| $O_2$ szinglett | -394764.823 | 1.216 |
| $O_2$ triplett | -394929.329 | 1.214 |
| $F_2$ szinglett | -524133.219 | 1.404 |
| $F_2$ triplett | -523980.926 | 1.719 |
Az 5.801Å-ös kötéstávolság rendkívül nagy. A kötés túl gyenge lenne, s így a molekula instabil lenne. A triplett-hidrogén a gyakorlatban ezért nem fordul elő.
\newpage
** Hidrogén HOMO és LUMO pályái
#+Caption: HOMO pálya
[[./h2e_h.png]]
\newpage
#+Caption: LUMO-1 pálya
[[./h2e_l.png]]
\newpage
#+Caption: LUMO-2 pálya
[[./h2e_l1.png]]
\newpage
#+Caption: LUMO-3 pálya
[[./h2e_l2.png]]
\newpage
** Rezgési analízis
$I$ Az infraaktivitás, és \nu a hullámszám.
*** Lineáris víz
| $\nu[cm^{-1}]$ | $I[km/mol]$ |
|----------------+-------------|
| -1733.29 | 564.319 |
| -1733.29 | 564.319 |
| 4073.74 | 0.000 |
| 4463.31 | 604.517 |
| E[kJ/mol] | kötéshossz[Å] | szög[°] |
|-------------+---------------+---------|
| -200597.195 | 0.9365 | 180.0 |
A $\nu=4073$ -es módusrezgés feltehetően a forgatási szimmetria a molekulával párhuzamos tengely körül, $I=0$ miatt. Ebbe az irányba "könnyű forgatni". Az $\nu=-1733$ -es hullámszámú módus degenerált, a másik két merőleges tengelyirányú forgatásokhoz tartozik, hiszen a molekula azokban az irányokban ugyanúgy viselkedik. A $\nu=4464$ -s módus a molekula rugózásához tartozik.
*** Normális víz
| $\nu[cm^{-1}]$ | $I[km/mol]$ |
|----------------+-------------|
| -1713.25 | 75.739 |
| -3729.65 | 1.698 |
| 3851.36 | 19.319 |
| E[kJ/mol] | kötéshossz[Å] | szög[°] |
|-------------+---------------+---------|
| -200745.930 | 0.969 | 103.6 |
A kötésszög és a kötéshossz az irodalmi érték[fn:1] közelében van. A szög alig tér el, 1% alatti a hiba.
\newpage
** Izotópeffektus
*** Normál
#+tblname: benzene
| $\nu_H[cm^{-1}]$ | $I_H[km/mol]$ |
|------------------+---------------|
| 414.76 | 0.000 |
| 415.57 | 0.000 |
| 621.97 | 0.000 |
| 622.01 | 0.000 |
| 717.70 | 0.000 |
| 864.37 | 0.000 |
| 865.15 | 0.000 |
| 969.11 | 0.000 |
| 969.62 | 0.000 |
| 1011.07 | 0.000 |
| 1019.93 | 0.000 |
| 1020.47 | 0.000 |
| 1185.79 | 0.000 |
| 1207.88 | 0.000 |
| 1208.04 | 0.000 |
| 1356.55 | 0.000 |
| 1387.34 | 0.000 |
| 1655.99 | 0.000 |
| 1656.13 | 0.000 |
| 3175.15 | 0.000 |
| 3184.75 | 0.000 |
| 3184.82 | 0.000 |
| 3211.23 | 0.001 |
| 1069.47 | 3.249 |
| 1069.06 | 3.263 |
| 3200.64 | 51.998 |
| 3200.50 | 52.002 |
| 1531.73 | 6.555 |
| 1531.38 | 6.562 |
| 694.56 | 77.548 |
#+call: plot(inputs=benzene, fname="fig1.png") :noexport:
\newpage
#+Caption: Benzol
[[./fig1.png]]
Ez már egy nagyobbacska molekula, így a modellünk pontatlanabb. A módusokat most nem a hullámszámmal jellemzem, hiszen azok alapján nem olyan egyszerű csoportosítani A $3.249 \frac{km}{mol}$ -es páros, az $52.002\frac{km}{mol}$ -es páros, illetve a $6.555\frac{km}{mol}$ -ös párosok degenerált, rotációs módusnak tűnnek. A nullásakkal most nem foglalkozom, ellenben a $77\frac{km}{mol}$ -es módus feltehetően egy rugózása a molekulának.
\newpage
*** Deutérium-benzol
Itt a hidrogénatomokat a program segítségével lecseréltük deutériumatomokra.
#+tblname: deuterium
| $\nu_D[cm^{-1}]$ | $I_D[km/mol]$ |
|------------------+---------------|
| 361.53 | 0.000 |
| 362.23 | 0.000 |
| 592.82 | 0.000 |
| 592.85 | 0.000 |
| 615.29 | 0.000 |
| 672.55 | 0.000 |
| 673.04 | 0.000 |
| 786.48 | 0.000 |
| 786.88 | 0.000 |
| 834.25 | 0.000 |
| 882.37 | 0.000 |
| 882.37 | 0.000 |
| 882.52 | 0.000 |
| 973.11 | 0.000 |
| 979.51 | 0.000 |
| 1079.31 | 0.000 |
| 1346.88 | 0.000 |
| 1612.36 | 0.000 |
| 1612.52 | 0.000 |
| 2338.71 | 0.000 |
| 2350.64 | 0.000 |
| 2350.68 | 0.000 |
| 2382.10 | 0.000 |
| 844.82 | 0.001 |
| 1377.46 | 1.341 |
| 1377.26 | 1.343 |
| 2370.10 | 29.070 |
| 2369.99 | 29.071 |
| 834.12 | 3.359 |
| 833.72 | 3.368 |
| 509.97 | 41.805 |
#+call: plot(inputs=deuterium, fname="fig2.png") :noexport:
#+Caption: Hydrogen
| multiplicity | E[kJ/mol] | dist[Å] |
| 1 | -3088.294 | 0.743 |
| 3 | -2628.693 | 5.801 |
#+Caption: Deutérium-benzol
[[./fig2.png]]
\newpage
** $C_6H_6$
A molekulákat az alakjuk alapján elneveztem.
- "Kurbli" szénlánca: $C\equiv C-C-C-C \equiv C$
- "Bari" szénlánca: $C=C-C\equiv C-C=C$
Nitrogen
| multiplicity | E[kJ/mol] | dist[Å] |
| 1 | -287747.996 | 1.106 |
| 3 | -287046.505 | 1.284 |
Keletkezési energiák [ $\frac{kJ}{mol}$ ] -ban:
| "Kurbli" | "Bari" |
| -609794.419 | -609889.479 |
Oxygen
| multiplicity | E[kJ/mol] | dist[Å] |
| 1 | -394764.823 | 1.216 |
| 3 | -394929.329 | 1.214 |
A két energia egymáshoz nagyon közel van, de ez nem csoda. Mindkettő hasonló struktúrájú.
* Diszkusszió
A mérés során végig az irodalmi értékek közelében levő eredményeket kaptunk, azonban ebben az is szerepet játszik, hogy a legtöbb irodalom hasonlót, vagy egyenesen ugyanazt a modellt használta a számításokhoz.
\newpage
** Feldolgozó program leírása
Először elkészítem az egyes tag-függvényeket legyártó kontruktor-függvényt. Egy closure-re egyszerű.
#+name: constructor
#+begin_src python
Γ=10
def get_distr(A, nu_0):
def distr(nu):
return A / ((nu - nu_0)**2 + (Γ/2)**2)
return distr
#+end_src
Fluorine
| multiplicity | E[kJ/mol] | dist[Å] |
| 1 | -524133.219 | 1.404 |
| 3 | -523980.926 | 1.719 |
Szűrök, és megfeleltetem a módusokat a disztribúcióknak.
A 0-s tagokat elhagyom:
#+name: funcmap
#+begin_src python :noweb yes
<<constructor>>
func_list = []
for row in inputs:
nu = row[0]
A = row[1]
if A > 0.001:
func_list.append(get_distr(A, nu))
#+end_src
Elkészítem az eredő függvényt, ami az egyes tagokat hívja meg, majd összeadja az eredményeket.
#+name: resultant
#+begin_src python :noweb yes
<<funcmap>>
def resulting(nu):
return sum([f(nu) for f in func_list])
#return resulting(694)
#+end_src
Plottolok.
#+name: plot
#+begin_src python :var inputs=benzene :var fname="" :noweb yes
<<resultant>>
import matplotlib.pyplot as plt
import numpy as np
nulist = np.linspace(450, 3500, 3000)
y = np.array([resulting(nu) for nu in nulist])
plt.plot(nulist, y, label='Az összegfüggvény')
plt.title('Konvolvált eredmény')
plt.xlim([500, 4000])
plt.ylim([0, 5.0])
plt.xlabel('ν[1/cm]')
plt.ylabel('I[km/mol]')
plt.legend()
plt.savefig(fname)
#+end_src
* Irodalom :noexport:
[fn:1] l=0.958Å és $\alpha=104.5$ °, forrás: [[http://www.idc-online.com/technical_references/pdfs/chemical_engineering/Water_molecule_structure.pdf]]
Hydrogen, electron stuff
| wat | q |
| | |

View File

@ -1,85 +0,0 @@
#+Title: Millikan-szimulátorok
Az esszém célja bemutatni néhány, könnyen hozzáférhető szimulátort, amelyekkel a kísérletet felszerelés hiányában is jól meg lehet ismerni.
* Bevezetés. Történelmi kontextus.
Az elektron létezésének első bizonyítékait ugyan már korábban is ismerte az emberiség, azonban a képtelenül nagy töltés-tömeg arány miatt, ami a protonénál (a kémia már ismert a hidrogén pozitív töltésű ionját) 1800-szor nagyob, de J.J. Thomson kísérletei nyomán elfogadottá vált, hogy az elektron egy részecske, minden vizsgált anyagban jelen van, és a mérés igazolta a korábbi becsléseket is a töltés-tömeg arányra.
Az elemi töltésre létezésére utaló mérések és elméleti eredmények már voltak: Max Planck megjósolta az elemi töltést a Faraday-állandó, illetve a Loschmidt-féle szám, továbbá a feketetest-sugárzás elméletének segítségével. Planck jóslata $e = 4.69 \cdot 10^{-10}$ /elektrosztatikus töltésegység/ volt, míg Richarz-é $1.29\cdot 10^{-10}$ , és J.J. Thomson-é $6.5\cdot 10^{-10}$ .[fn:1]. Természetesen ekkor még nem lehetett tudni, hogy az elektron töltése lesz az, amit mi ma elemi töltésnek nevezünk, ezt csak utólag tudjuk.
Az első ilyen mérést Robert Millikan és Harvey Fletcher végezték, magyar kultúrterületen Millikan után neveztük el a kísérletet, angolul pedig /oil drop experiment/ ismert a fogalom. A magyar elnevezés ellenére Fletchernek is fontos szerepe volt. Kísérlet több szempontból is nagy jelentőséggel rendelkezett:
- egy nagyon pontos mérést adott az elektron töltéséről (így a töltés-tömeg arányból rögtön kaptunk egy eredményt a tömegére, s így megismerhettük a mozgását uraló számokat is)
- egy iskolapéldáját adta annak, hogy a megerősítési torzítás (/confirmation bias/) miféle módokon hat ki a korábban már mért eredmények újramérésére[fn:2]
Az első eredmények pontatlanságának több oka is volt, pl. a kísérletben Millikanék még nem alkalmazták azokat a korrekciókat, amiket mi a Modern Fizika Laboratórium mérésben felhasználtunk az eredményeink pontosabbá tételére.
* [[http://www.college-physics.com/book/electric-field/oil-drop-experiment/#the-experiment][A college physics szimulátora]]
#+CAPTION: A College Physics szimulátorának a kezelőfelülete
#+NAME: fig:cphy
[[~/college-physics.png]]
A mérés fizikai elméletéhez a weblap bőséges információt ad, levezetésekkel.
A kezelőfelület már majdnem spártai, de nem teljesen egyértelmű a használata, és a leírás ebben nem ad támpontot. A "run" gomb benyomásával indíthatjuk el a szimulációt. A "spray" gomb csakis akkor csinál valamit, ha a szoftver "fut", azaz "running" állapotban van. Ekkor is csak némi késéssel szór be véletlenszerű mennyiségű cseppet, *ami lehet 0 is*. Ha sorozatosan klikkelünk a "spray" gombra, akkor előbb-utóbb inaktívvá válik, tehát meg van adva egy felső limit az olajcseppek számához. Ez a kísérlet elvégzése szemponjából jelentéktelen, mert a szoftver csakis olyan cseppeket generál, amelyeket könnyűszerrel meg lehet mérni. Érdekes /bug/, hogy ha megnyomjuk a "pause" gombot a "weiter" felirat jelenik meg, ami német nyelvű, és azt jelenti ebben a kontextusban, hogy "folytasd". Ez egy fordítási hiányosság, hiszen az angolban a "continue" lenne a megfelelő szó. Ez alapján valószínűnek tartom, hogy ez egy német szoftver adaptációja.
A felület úgy van kialakítva, és egyedül ezt a módszert is ajánlja, hogy próbáljuk meg addig variálni a feszültséget,
A kiértékeléshez egyetlen támpontot kapunk, de az jobban emlékeztet a "ránézésre" mérésekre, mint a mi laborunkban leírt módszerhez. A képet a wikipediáról vették, és ezt le is hivatkozták.
#+CAPTION: A kiértékelést szemléltető kép
#+NAME: fig:ladung
[[~/millikan-ladung.png]]
* [[https://ophysics.com/em2.html][Az OPhysics szimulátora]]
#+CAPTION: Az OPhysics szimulátorának a kezelőfelülete
#+NAME: fig:cphy
[[~/ophysics.png]]
Ez a szoftver inkább középiskolai demonstrációs célt szolgál, mint egy egyetemi laboratóriumi mérés szimulációját felszerelés hiányában. A használt képletek megválasztása alapján kijelenthető, hogy próbáltak az ennek a szintnek megfelelő anyagra hagyatkozni. A mérési leírása alapos, és a képletek mellé sok szöveges magyarázatot nyújt.
A kezelőfelület kőbunkó egyszerűségű, és magától értetődő is. A szoftver már induláskor legenerált egy részecskét, amelynek a helyébe gombbal újat generálhatunk. A szimulációt elindítani a "start" gombbal tudjuk, és megállíthatjuk a "pause" gombbal, illetve a bal alsó sarokban található egyezményes szimbólumú gombbal is, ami egy "toggle" típusú gomb, tehát kattintás után funkciót vált és "start"-ként működik. A gomb eltűnik, ha a "pause" gombbal állítjuk meg a szimulációt. Tehát csakis akkor láthatjuk megállított állapot esetén, ha magát ezt a gombot használtuk a szüneteltetéshez.
A feszültséghez ugyan nem írtak mértékegységet, de a leírás alapján kitalálható, hogy voltban mutatja a csúszka. Sajnos ennek a felbontása nem elég sűrű, tehát nem lehet minden csepp esetén nullázni a végsebességet. Minthogy a feladatleírás mégis ezt kéri tőlünk, ezért feltehetően addig érdemes generálni az újabb cseppeket, amíg nem kapunk egy, a célnak megfelelő alanyt.
* [[https://www.thephysicsaviary.com/Physics/Programs/Labs/MillikanOilDropLab/][A "Millikan Oil Drop Lab" szimulátor]]
#+CAPTION: A Millikan Oil Drop Lab kezdőlapja
#+NAME: fig:cphy
[[~/oil-drop-lab-1.png]]
Ez a szoftver az előzőhöz hasonlóan középiskolai órákhoz ajánlható. A mérési leírás elegendő a megértéshez és a feladat elvégzéséhez.
A program más módon működik, mint az eddigiek. Az olajcseppeket a pumpára kattintással szórhatjuk be, és minden esetben lesz egy olyan csepp, amelyik stabil. A "mikroszkópra" kattintva -- azért nevezem mikroszkópnak, mert a valóságban cseppek igen aprók és közönséges nagyítóval nem mérhetők le pontosan -- a csepp mérhető méretűre lesz kinagyítva, és a felhasználó feladata ez alapján a tömeget kiszámolni. Az átmérő véletlenszerűen változik minden egyes befecskendezésnél. A mikroszkóp felületénél kapjuk meg az elektromos mezőre vonatkozó adatokat is.
#+CAPTION: A Millikan Oil Drop Lab "mikroszkópja"
#+NAME: fig:cphy
[[~/oil-drop-lab-2.png]]
A tipp, miszerint az osztály összegezheti a mérési eredményeit, iskolai környezetben pedagógiailag hasznos tud lenni, azonban a mi kontextusunkban -- egyetemi laboratóriumban, amiatt a kettős elvárás miatt, hogy az oktatók értékeljék az egyes a diákok munkáját, illetve teljesítményét, illetve a diákok teljesítsenek, továbbá amiatt, hogy a diákokat alkalmassá tegyük a mérés egyedüli elvégzésére is -- nehezen alkalmazható. Egy ilyen módszertan inkább egy előadás keretei közt tud hasznosulni, amennyiben az oktató szeretne a /frontális oktatás/ módszertanától eltérni.
* [[https://vlab.amrita.edu/?sub=1&brch=195&sim=357&cnt=1][A Virtual Amrita Laboratories szimulátora]] [fn:4]
#+CAPTION: A Millikan Oil Drop Lab "mikroszkópja"
#+NAME: fig:cphy
[[~/amritsar.png]]
Ez a szoftver meglehetősen összetett az előzőkhöz képest, egy egyetemi laboratóriumi mérés összecsomagolva egybe. Különöző fülekre elhelyezve tartalmaz egy elméleti leírást, egy használati leírást, egy "beugrót", a szimulátort, mérési feladatokat, irodalmi hivatkozásokat, és visszajelzési lehetőséget is.
A szimulátorban külön súgófunkció van, azonban a "start" gombot nem találjuk meg, amíg nem tesszük teljes képernyőre a szimulációt. Innentől viszont egyértelmű a kezelés. A "voltage ON" gomb azt jelzi, hogy jelenleg épp be van kapcsolva a feszültség. A gombbal kapcsolhatjuk *ki* a cseppekre a feszültséget, és ekkor a gomb neve átvált "voltage OFF"-ra. A feszültséget csakis kikapcsolt állapotban engedi a program állítani.
A csöppeket a program automatikusan "betereli" a nagyító elé, ezek közül válaszhatunk.
A szoftver kétféle folyadékot is kínál cseppekhez: olívaolajat és glicerint, továbbá ionizálhatjuk a közeget röntgensugárzással is. Ez megakadályozza, hogy az olajcseppek leessenek, illetve visszalöki a más hulló cseppeket a tartályba.
[fn:1] Max Planck: Válogatott tanulmányok. Gondolat, 1982. 356-357
Az /elektrosztatikus töltésegység/ fogalom a [[https://en.wikipedia.org/wiki/Statcoulomb][CGS-mértékegység]], de nem ugyanazt a fogalmat fejezik ki, így a váltás nem egyértelmű. Az $1C$ -ot az egyszerűség kedvéért közelíthetjük $3.00 \cdot 10^9 statC$ -bal.
Megjegyzés: az esszé írásakor Thomson eredményére a [[https://en.wikipedia.org/w/index.php?title=Electron&oldid=1145699944][Wikipedia]] $6.8\cdot 10^{-10}$ -t írt. Se a Planck-szöveg, se a Wikipedia nem tekinthető elsődleges forrásnak, de Planck közvetlenül Thomson-t hivatkozza: J.J. Thomson: Phil. Mag. 46, 528 (1898). Az esszé témájából fakadóan mélyebb forráskritika nem szükségeltetik, ezért a kíváncsi Olvasóra bízom.
[fn:2] Richard P. Feynman: "Tréfál, Feynman úr?" Park, 2001. 334
[fn:3] By Christian Hill - https://scipython.com/blog/measurements-of-the-electron-charge-over-time/, CC BY-SA 4.0, https://commons.wikimedia.org/w/index.php?curid=114882411
[fn:4] vlab.amrita.edu,. (2011). Millikan's oil drop experiment. Retrieved 1 April 2023, from vlab.amrita.edu/?sub=1&brch=195&sim=357&cnt=1

View File

@ -1,10 +0,0 @@
90 1.0
80 0.98481
70 0.93969
60 0.86603
50 0.76604
40 0.64279
30 0.50000
20 0.34202
10 0.17365
0 0.000

View File

@ -1,34 +0,0 @@
#+OPTIONS: tex:t
#+AUTHOR: Rédl Anna, Barna Zsombor
#+DATE: 2023. 07. 25.
#+TITLE: Összegyűjtött
#+LATEX_CLASS: article
#+LATEX_CLASS_OPTIONS: [12pt,a4paper]
#+LATEX_HEADER: \usepackage[utf8]{inputenc}
#+LATEX_HEADER: \usepackage[magyar]{babel}
#+LATEX_HEADER: \hypersetup{ colorlinks=true,linkcolor=blue, linktoc=all, filecolor=magenta, urlcolor=cyan, pdfstartview=FitB,}
#+LATEX_HEADER: \usepackage{multirow}
#+LATEX_HEADER: \usepackage{braket}
#+LATEX_HEADER: \urlstyle{same}
#+OPTIONS: toc:nil
#+tblname: deuterium
| $\nu_D[cm^{-1}]$ | $I_D[km/mol]$ |
|------------------+---------------|
| 361.53 | 0.000 |
| 362.23 | 0.000 |
| 592.82 | 0.000 |
| 592.85 | 0.000 |
| 615.29 | 0.000 |
#+call: plot(inputs=deuterium, fname="fig2.png") :noexport:
#+name: plot
#+begin_src python :var inputs=benzene :var fname="" :noweb yes
func_list = []
for row in inputs:
nu = row[0]
A = row[1]
if A > 0.001:
func_list.append(get_distr(A, nu))
#+end_src

View File

@ -1,36 +0,0 @@
##################################################
# Author: József Stéger
# Date: 02 13 2017
# Usage (tested in linux):
# make
# Summary:
# 1. compiles objects in ./build/
# 2. creates lib archive in .
##################################################
UTILS=datafit fft interpolate matrix odeint optimize random vector
BUILD=build
OBJ_LIBS=$(patsubst %,$(BUILD)/%.o,$(UTILS))
A_LIB=libcpl.a
.PHONY: all distclean clean
all: $(A_LIB)
clean:
rm -rf $(A_LIB) $(BUILD)
distclean: clean
rm -rf $(PDIR)
$(OBJ_LIBS):
@make $(BUILD)
g++ -Wall -c -I. -O2 $(patsubst $(BUILD)/%.o,./%.cpp,$@) -o $@
$(A_LIB): $(OBJ_LIBS)
ar -srcl $@ $?
$(BUILD):
@mkdir -p $@

Binary file not shown.

View File

@ -1,178 +0,0 @@
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <iostream>
#include <limits>
#include <vector>
using namespace std;
#include "vector.hpp"
#include "datafit.hpp"
namespace cpl {
inline double SQR(double a) {
return a * a;
}
static void error(const char *msg) {
cerr << msg << endl;
exit(EXIT_FAILURE);
}
double gammln(const double xx)
{
int j;
double x,y,tmp,ser;
static const double cof[6]=
{76.18009172947146,-86.50532032941677,
24.01409824083091,-1.231739572450155,0.1208650973866179e-2,
-0.5395239384953e-5};
y=x=xx;
tmp=x+5.5;
tmp -= (x+0.5)*log(tmp);
ser=1.000000000190015;
for (j=0;j<6;j++) ser += cof[j]/++y;
return -tmp+log(2.5066282746310005*ser/x);
}
void gser(double& gamser, const double a, const double x, double& gln)
{
const int ITMAX=100;
const double EPS=numeric_limits<double>::epsilon();
int n;
double sum,del,ap;
gln=gammln(a);
if (x <= 0.0) {
if (x < 0.0) error("x less than 0 in routine gser");
gamser=0.0;
return;
} else {
ap=a;
del=sum=1.0/a;
for (n=0;n<ITMAX;n++) {
++ap;
del *= x/ap;
sum += del;
if (abs(del) < abs(sum)*EPS) {
gamser=sum*exp(-x+a*log(x)-gln);
return;
}
}
error("a too large, ITMAX too small in routine gser");
return;
}
}
void gcf(double& gammcf, const double a, const double x, double& gln)
{
const int ITMAX=100;
const double EPS=numeric_limits<double>::epsilon();
const double FPMIN=numeric_limits<double>::min()/EPS;
int i;
double an,b,c,d,del,h;
gln=gammln(a);
b=x+1.0-a;
c=1.0/FPMIN;
d=1.0/b;
h=d;
for (i=1;i<=ITMAX;i++) {
an = -i*(i-a);
b += 2.0;
d=an*d+b;
if (abs(d) < FPMIN) d=FPMIN;
c=b+an/c;
if (abs(c) < FPMIN) c=FPMIN;
d=1.0/d;
del=d*c;
h *= del;
if (abs(del-1.0) <= EPS) break;
}
if (i > ITMAX) error("a too large, ITMAX too small in gcf");
gammcf=exp(-x+a*log(x)-gln)*h;
}
double gammq(const double a, const double x)
{
double gamser,gammcf,gln;
if (x < 0.0 || a <= 0.0)
error("Invalid arguments in routine gammq");
if (x < a+1.0) {
gser(gamser,a,x,gln);
return 1.0-gamser;
} else {
gcf(gammcf,a,x,gln);
return gammcf;
}
}
void fit(
const Vector& x,
const Vector& y,
const Vector& sig,
const bool mwt,
double& a,
double& b,
double& siga,
double& sigb,
double& chi2,
double& q)
{
int i;
double wt,t,sxoss,sx=0.0,sy=0.0,st2=0.0,ss,sigdat;
int ndata=x.size();
b=0.0;
if (mwt) {
ss=0.0;
for (i=0;i<ndata;i++) {
wt=1.0/SQR(sig[i]);
ss += wt;
sx += x[i]*wt;
sy += y[i]*wt;
}
} else {
for (i=0;i<ndata;i++) {
sx += x[i];
sy += y[i];
}
ss=ndata;
}
sxoss=sx/ss;
if (mwt) {
for (i=0;i<ndata;i++) {
t=(x[i]-sxoss)/sig[i];
st2 += t*t;
b += t*y[i]/sig[i];
}
} else {
for (i=0;i<ndata;i++) {
t=x[i]-sxoss;
st2 += t*t;
b += t*y[i];
}
}
b /= st2;
a=(sy-sx*b)/ss;
siga=sqrt((1.0+sx*sx/(ss*st2))/ss);
sigb=sqrt(1.0/st2);
chi2=0.0;
q=1.0;
if (!mwt) {
for (i=0;i<ndata;i++)
chi2 += SQR(y[i]-a-b*x[i]);
sigdat=sqrt(chi2/(ndata-2));
siga *= sigdat;
sigb *= sigdat;
} else {
for (i=0;i<ndata;i++)
chi2 += SQR((y[i]-a-b*x[i])/sig[i]);
if (ndata>2) q=gammq(0.5*(ndata-2),0.5*chi2);
}
}
} /* end namespace cpl */

View File

@ -1,23 +0,0 @@
#ifndef CPL_DATAFIT_HPP
#define CPL_DATAFIT_HPP
namespace cpl {
class Vector;
extern void fit( // Least Squares fit to y(x) = a + b x
const Vector& x, // Input: Vector of x values
const Vector& y, // Input: Vector of y values
const Vector& sig, // Input: Vector of individual standard deviations
const bool mwt, // Input: if false sig's assumed to be unavailable
double& a, // Output: intercept a
double& b, // Output: slope b
double& siga, // Output: uncertainty in a
double& sigb, // Output: uncertainty in b
double& chi2, // Output: chi-square
double& q // Output: goodness of fit
);
} /* end namespace cpl */
#endif /* CPL_DATAFIT_HPP */

View File

@ -1,214 +0,0 @@
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <iostream>
#include <vector>
using namespace std;
#include "fft.hpp"
#include "vector.hpp"
namespace cpl {
ComplexVector::ComplexVector(int dim) {
v = new std::complex<double> [this->dim = dim];
for (int i = 0; i < dim; i++) v[i] = 0.0;
}
ComplexVector::ComplexVector(const ComplexVector& cv) {
v = new std::complex<double> [dim = cv.dim];
for (int i = 0; i < dim; i++) v[i] = cv.v[i];
}
ComplexVector& ComplexVector::operator = (const ComplexVector& cv) {
if (this != &cv) {
if (dim != cv.dim) {
delete [] v;
v = new std::complex<double> [dim = cv.dim];
}
for (int i = 0; i < dim; i++) v[i] = cv[i];
}
return *this;
}
// FFT implementation
void FFT::transform(ComplexVector& data) {
N = data.dimension();
f = &data;
bitReverse();
for (int n = 1; n < N; n *= 2)
DanielsonLanczos(n);
for (int i = 0; i < N; ++i)
(*f)[i] /= std::sqrt(double(N));
}
void FFT::inverseTransform(ComplexVector& data) {
inverse = true;
transform(data);
inverse = false;
}
void FFT::bitReverse() {
int j = 1;
for (int i = 1; i < N; ++i) {
if (i < j) {
std::complex<double> temp = (*f)[i-1];
(*f)[i-1] = (*f)[j-1];
(*f)[j-1] = temp;
}
int k = N / 2;
while ( k < j ) {
j -= k;
k /= 2;
}
j += k;
}
}
void FFT::DanielsonLanczos(int n) {
const double pi = 4 * atan(1.0);
std::complex<double> W(0, pi / n);
W = inverse ? std::exp(W) : std::exp(-W);
std::complex<double> W_j(1, 0);
for (int j = 0; j < n; ++j) {
for (int i = j; i < N; i += 2 * n) {
std::complex<double> temp = W_j * (*f)[n+i];
(*f)[n+i] = (*f)[i] - temp;
(*f)[i] += temp;
}
W_j *= W;
}
}
Vector FFT::power(ComplexVector& data) {
Vector P(1 + N / 2);
P[0] = std::norm(data[0]) / double(N);
for (int i = 1; i < N / 2; i++)
P[i] = (std::norm(data[i]) + std::norm(data[N-i])) / double(N);
P[N/2] = std::norm(data[N/2]) / double(N);
return P;
}
static void error(const std::string str) {
std::cerr << "cpl::ComplexMatrix error: " << str << std::endl;
std::exit(EXIT_FAILURE);
}
static void error(const std::string str, int i) {
std::ostringstream os;
os << str << " " << i;
error(os.str());
}
ComplexMatrix::ComplexMatrix(int rows, int cols, complex<double> d) {
if (rows <= 0)
error("bad number of rows", rows);
if (cols <= 0)
error("bad number of columns", cols);
this->rows = rows;
this->cols = cols;
m = new complex<double>* [rows];
for (int i = 0; i < rows; i++) {
m[i] = new complex<double> [cols];
for (int j = 0; j < cols; j++)
m[i][j] = d;
}
}
ComplexMatrix::ComplexMatrix(const ComplexMatrix& mat) {
rows = mat.rows;
cols = mat.cols;
m = new complex<double>* [rows];
for (int i = 0; i < rows; i++) {
m[i] = new complex<double> [cols];
for (int j = 0; j < cols; j++)
m[i][j] = mat[i][j];
}
}
ComplexMatrix::~ComplexMatrix() {
for (int i = 0; i < rows; i++)
delete [] m[i];
delete [] m;
}
const complex<double>* ComplexMatrix::operator[](const int row) const {
if (row < 0 || row >= rows)
error("bad row index", row);
return m[row];
}
complex<double>* ComplexMatrix::operator[](const int row) {
if (row < 0 || row >= rows)
error("bad row index", row);
return m[row];
}
ComplexMatrix& ComplexMatrix::operator=(const ComplexMatrix& mat) {
if (this != &mat) {
if (rows != mat.rows || cols != mat.cols) {
for (int i = 0; i < rows; i++)
delete [] m[i];
delete [] m;
rows = mat.rows;
cols = mat.cols;
m = new complex<double>* [rows];
for (int i = 0; i < rows; i++)
m[i] = new complex<double> [cols];
}
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
m[i][j] = mat[i][j];
}
return *this;
}
void FFT2::transform(ComplexMatrix& data) {
FFT fft;
int rows = data.numRows();
int cols = data.numCols();
ComplexVector r(cols), c(rows);
// FFT rows of data
for (int j = 0; j < rows; j++) {
for (int k = 0; k < cols; k++)
r[k] = data[j][k];
fft.transform(r);
for (int k = 0; k < cols; k++)
data[j][k] = r[k];
}
// FFT columns of data
for (int k = 0; k < cols; k++) {
for (int j = 0; j < rows; j++)
c[j] = data[j][k];
fft.transform(c);
for (int j = 0; j < rows; j++)
data[j][k] = c[j];
}
}
void FFT2::inverseTransform(ComplexMatrix& data) {
FFT fft;
int rows = data.numRows();
int cols = data.numCols();
ComplexVector r(cols), c(rows);
// FFT rows of data
for (int j = 0; j < rows; j++) {
for (int k = 0; k < cols; k++)
r[k] = data[j][k];
fft.inverseTransform(r);
for (int k = 0; k < cols; k++)
data[j][k] = r[k];
}
// FFT columns of data
for (int k = 0; k < cols; k++) {
for (int j = 0; j < rows; j++)
c[j] = data[j][k];
fft.inverseTransform(c);
for (int j = 0; j < rows; j++)
data[j][k] = c[j];
}
}
} /* end namespace cpl */

View File

@ -1,75 +0,0 @@
#ifndef CPL_FFT_HPP
#define CPL_FFT_HPP
#include <complex>
#include <iostream>
namespace cpl {
class Vector;
class ComplexVector {
public:
ComplexVector(int dim = 1);
ComplexVector(const ComplexVector& cv);
~ComplexVector() { delete [] v; }
int dimension() const { return dim; }
const std::complex<double> operator[](const int i) const { return v[i]; }
std::complex<double>& operator[](const int i) { return v[i]; }
ComplexVector& operator = (const ComplexVector& cv);
private:
int dim;
std::complex<double> *v;
};
class FFT {
public:
FFT() { N = 0; f = 0; inverse = false; }
void transform(ComplexVector& data);
void inverseTransform(ComplexVector& data);
Vector power(ComplexVector& data);
private:
int N;
ComplexVector *f;
bool inverse;
void bitReverse();
void DanielsonLanczos(int n);
};
class ComplexMatrix {
public:
ComplexMatrix(int rows=1, int cols=1, complex<double> d=0);
ComplexMatrix(const ComplexMatrix& m);
~ComplexMatrix();
int numRows() const { return rows; }
int numCols() const { return cols; }
const std::complex<double>* operator[](const int row) const;
std::complex<double>* operator[](const int row);
std::complex<double>& operator()(int row, int col);
ComplexMatrix& operator=(const ComplexMatrix& m);
private:
int rows;
int cols;
std::complex<double> **m;
};
class FFT2 {
public:
FFT2() { }
void transform(ComplexMatrix& data);
void inverseTransform(ComplexMatrix& data);
private:
};
} /* end namespace cpl */
#endif /* CPL_FFT_HPP */

View File

@ -1,59 +0,0 @@
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <iostream>
#include <vector>
using namespace std;
#include "interpolate.hpp"
namespace cpl {
// polynomial interpolation routine adapted from Numerical Recipes
// polynomial degree inferred from dimension of input vector
void Interpolator::polint( // uses Neville's algorithm
const vector<double>& xa, // vector of abscissa values
const vector<double>& ya, // vector of ordinate values
const double x, // point x
double& y, // interpolated y(x)
double& dy) // estimate of error in interpolation
{
int i,m,ns=0;
double den,dif,dift,ho,hp,w;
int n=xa.size();
vector<double> c(n),d(n);
dif=abs(x-xa[0]);
for (i=0;i<n;i++) {
if ((dift=abs(x-xa[i])) < dif) {
ns=i;
dif=dift;
}
c[i]=ya[i];
d[i]=ya[i];
}
y=ya[ns--];
for (m=1;m<n;m++) {
for (i=0;i<n-m;i++) {
ho=xa[i]-x;
hp=xa[i+m]-x;
w=c[i+1]-d[i];
if ((den=ho-hp) == 0.0) { // nrerror("Error in routine polint");
cerr << "Error in Interpolator::polint\n"
<< "Neville's algorithm breaks down when two "
<< "input x values are equal within roundoff\n"
<< "x, y values" << endl;
for (int j = 0; j < n; j++)
cerr << xa[j] << '\t' << ya[j] << '\n';
exit(EXIT_FAILURE);
}
den=w/den;
d[i]=hp*den;
c[i]=ho*den;
}
y += (dy=(2*(ns+1) < (n-m) ? c[ns+1] : d[ns--]));
}
}
} /* end namespace cpl */

View File

@ -1,19 +0,0 @@
#ifndef CPL_INTERPOLATE_HPP
#define CPL_INTERPOLATE_HPP
namespace cpl {
class Interpolator {
public:
static void polint( // adapted from Numerical Recipes
const vector<double>& xa, // vector of x values - input
const vector<double>& ya, // vector of y values - input
const double x, // x at which to find y - input
double& y, // value of y - output
double& dy); // estimate or error in y
};
} /* end namespace cpl */
#endif /* CPL_INTERPOLATE_HPP */

View File

@ -1,895 +0,0 @@
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <iostream>
#include <limits>
#include <string>
#include <sstream>
#include "matrix.hpp"
#include "vector.hpp"
namespace cpl {
static void error(const std::string str) {
std::cerr << "cpl::Matrix error: " << str << std::endl;
std::exit(EXIT_FAILURE);
}
static void error(const std::string str, int i) {
std::ostringstream os;
os << str << " " << i;
error(os.str());
}
Matrix::Matrix(int rows, int cols, double d) {
if (rows <= 0)
error("bad number of rows", rows);
if (cols <= 0)
error("bad number of columns", cols);
this->rows = rows;
this->cols = cols;
m = new double* [rows];
for (int i = 0; i < rows; i++) {
m[i] = new double [cols];
for (int j = 0; j < cols; j++)
m[i][j] = d;
}
}
Matrix::Matrix(const Matrix& mat) {
rows = mat.rows;
cols = mat.cols;
m = new double* [rows];
for (int i = 0; i < rows; i++) {
m[i] = new double [cols];
for (int j = 0; j < cols; j++)
m[i][j] = mat[i][j];
}
}
Matrix::~Matrix() {
for (int i = 0; i < rows; i++)
delete [] m[i];
delete [] m;
}
const double* Matrix::operator[](const int row) const {
if (row < 0 || row >= rows)
error("bad row index", row);
return m[row];
}
double* Matrix::operator[](const int row) {
if (row < 0 || row >= rows)
error("bad row index", row);
return m[row];
}
Matrix& Matrix::operator=(const Matrix& mat) {
if (this != &mat) {
if (rows != mat.rows || cols != mat.cols) {
for (int i = 0; i < rows; i++)
delete [] m[i];
delete [] m;
rows = mat.rows;
cols = mat.cols;
m = new double* [rows];
for (int i = 0; i < rows; i++)
m[i] = new double [cols];
}
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
m[i][j] = mat[i][j];
}
return *this;
}
Matrix& Matrix::operator=(const double d) {
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
m[i][j] = d;
return *this;
}
Matrix& Matrix::operator+=(const Matrix& mat) {
if (this != &mat) {
if (rows != mat.rows || cols != mat.cols)
error("matrix dimension mismatch");
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
m[i][j] += mat[i][j];
}
return *this;
}
Matrix& Matrix::operator-=(const Matrix& mat) {
if (this != &mat) {
if (rows != mat.rows || cols != mat.cols)
error("matrix dimension mismatch");
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
m[i][j] -= mat[i][j];
}
return *this;
}
Matrix& Matrix::operator*=(const double d) {
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
m[i][j] *= d;
return *this;
}
Matrix& Matrix::operator/=(const double d) {
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
m[i][j] /= d;
return *this;
}
Matrix operator - (const Matrix& mat) {
int rows = mat.numRows();
int cols = mat.numCols();
Matrix temp(rows, cols);
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
temp[i][j] = -mat[i][j];
return temp;
}
Matrix operator * (const Matrix& mat, double d) {
int rows = mat.numRows();
int cols = mat.numCols();
Matrix temp(rows, cols);
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
temp[i][j] = mat[i][j] * d;
return temp;
}
Matrix operator * (double d, const Matrix& mat) {
int rows = mat.numRows();
int cols = mat.numCols();
Matrix temp(rows, cols);
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
temp[i][j] = d * mat[i][j];
return temp;
}
Matrix operator / (const Matrix& mat, double d) {
int rows = mat.numRows();
int cols = mat.numCols();
Matrix temp(rows, cols);
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
temp[i][j] = mat[i][j] / d;
return temp;
}
Matrix operator + (const Matrix& m1, const Matrix& m2) {
int rows = m1.numRows();
int cols = m1.numCols();
if (rows != m2.numRows() || cols != m2.numCols())
error("matrix dimension mismatch");
Matrix temp(rows, cols);
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
temp[i][j] = m1[i][j] + m2[i][j];
return temp;
}
Matrix operator - (const Matrix& m1, const Matrix& m2) {
int rows = m1.numRows();
int cols = m1.numCols();
if (rows != m2.numRows() || cols != m2.numCols())
error("matrix dimension mismatch");
Matrix temp(rows, cols);
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
temp[i][j] = m1[i][j] - m2[i][j];
return temp;
}
Matrix operator * (const Matrix& m1, const Matrix& m2) {
int n = m1.numCols();
if (n != m2.numRows())
error("matrix dimension mismatch");
int rows = m1.numRows();
int cols = m2.numCols();
Matrix temp(rows, cols);
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
for (int k = 0; k < n; k++)
temp[i][j] += m1[i][k] * m2[k][j];
return temp;
}
Matrix Matrix::transpose() {
Matrix temp(cols, rows);
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
temp[j][i] = m[i][j];
return temp;
}
std::ostream& operator<<(std::ostream& os, const Matrix& mat) {
for (int i = 0; i < mat.rows; i++) {
for (int j = 0; j < mat.cols; j++) {
os << '\t' << mat.m[i][j];
}
os << '\n';
}
return os;
}
void solveGaussJordan(Matrix& A, Matrix& B) {
int n = A.numRows();
if (A.numCols() != n)
error("matrix dimension mismatch");
if (B.numRows() != n)
error("matrix dimension mismatch");
int m = B.numCols();
int *indxc = new int [n];
int *indxr = new int [n];
int *ipiv = new int [n];
for (int j = 0; j < n; j++)
ipiv[j] = 0;
for (int i = 0; i < n; i++) {
double big = 0;
int irow, icol;
for (int j = 0; j < n; j++) {
if (ipiv[j] != 1) {
for (int k = 0; k < n; k++) {
if (ipiv[k] == 0) {
double Ajk = A[j][k];
if (std::abs(Ajk) >= big) {
big = std::abs(Ajk);
irow = j;
icol = k;
}
}
}
}
}
++ipiv[icol];
if (irow != icol) {
for (int j = 0; j < n; j++)
std::swap(A[irow][j], A[icol][j]);
for (int j = 0; j < m; j++)
std::swap(B[irow][j], B[icol][j]);
}
indxr[i] = irow;
indxc[i] = icol;
if (A[icol][icol] == 0.0)
error("solveGaussJordan singular matrix");
double pivinv = 1 / A[icol][icol];
A[icol][icol] = 1;
for (int j = 0; j < n; j++)
A[icol][j] *= pivinv;
for (int j = 0; j < m; j++)
B[icol][j] *= pivinv;
for (int j = 0; j < n; j++) {
if (j != icol) {
double mult = A[j][icol];
for (int k = 0; k < n; k++)
A[j][k] -= A[icol][k] * mult;
for (int k = 0; k < m; k++)
B[j][k] -= B[icol][k] * mult;
}
}
}
for (int i = n - 1; i >= 0; i--) {
if (indxr[i] != indxc[i]) {
for (int j = 0; j < n; j++)
std::swap(A[j][indxr[i]], A[j][indxc[i]]);
}
}
delete [] indxc;
delete [] indxr;
delete [] ipiv;
}
void ludcmp(Matrix& A, int *indx, double& d) {
const double TINY=1.0e-20;
int i,imax,j,k;
double big,dum,sum,temp;
int n=A.numRows();
Vector vv(n);
d=1.0;
for (i=0;i<n;i++) {
big=0.0;
for (j=0;j<n;j++)
if ((temp=std::abs(A[i][j])) > big) big=temp;
if (big == 0.0) error("Singular matrix in routine ludcmp");
vv[i]=1.0/big;
}
for (j=0;j<n;j++) {
for (i=0;i<j;i++) {
sum=A[i][j];
for (k=0;k<i;k++) sum -= A[i][k]*A[k][j];
A[i][j]=sum;
}
big=0.0;
for (i=j;i<n;i++) {
sum=A[i][j];
for (k=0;k<j;k++) sum -= A[i][k]*A[k][j];
A[i][j]=sum;
if ((dum=vv[i]*std::abs(sum)) >= big) {
big=dum;
imax=i;
}
}
if (j != imax) {
for (k=0;k<n;k++) {
dum=A[imax][k];
A[imax][k]=A[j][k];
A[j][k]=dum;
}
d = -d;
vv[imax]=vv[j];
}
indx[j]=imax;
if (A[j][j] == 0.0) A[j][j]=TINY;
if (j != n-1) {
dum=1.0/(A[j][j]);
for (i=j+1;i<n;i++) A[i][j] *= dum;
}
}
}
void lubksb(const Matrix& A, int *indx, Matrix& B) {
int i,ii=0,ip,j;
double sum;
int n=A.numRows();
for (int k = 0; k < B.numCols(); k++) {
for (i=0;i<n;i++) {
ip=indx[i];
sum=B[ip][k];
B[ip][k]=B[i][k];
if (ii != 0)
for (j=ii-1;j<i;j++) sum -= A[i][j]*B[j][k];
else if (sum != 0.0)
ii=i+1;
B[i][k]=sum;
}
for (i=n-1;i>=0;i--) {
sum=B[i][k];
for (j=i+1;j<n;j++) sum -= A[i][j]*B[j][k];
B[i][k]=sum/A[i][i];
}
}
}
void solveLUDecompose(Matrix& A, Matrix& B) {
int n = A.numRows();
if (A.numCols() != n)
error("matrix dimension mismatch");
if (B.numRows() != n)
error("matrix dimension mismatch");
int m = B.numCols();
int *indx = new int [n];
double d;
ludcmp(A, indx, d);
lubksb(A, indx, B);
delete [] indx;
}
void reduceHouseholder(Matrix& A, Vector& d, Vector& e) {
int n = d.dimension();
for (int i = n - 1; i > 0; i--) {
int l = i - 1;
double h = 0;
double scale = 0;
if (l > 0) {
for (int k = 0; k <= l; k++)
scale += std::abs(A[i][k]);
if (scale == 0.0)
e[i] = A[i][l];
else {
for (int k = 0; k <= l; k++) {
A[i][k] /= scale;
h += A[i][k] * A[i][k];
}
double f = A[i][l];
double g = (f >= 0.0 ? -std::sqrt(h) : std::sqrt(h));
e[i] = scale * g;
h -= f * g;
A[i][l] = f - g;
f = 0.0;
for (int j = 0; j <= l; j++) {
A[j][i] = A[i][j] / h;
g = 0.0;
for (int k = 0; k <= j; k++)
g += A[j][k] * A[i][k];
for (int k = j + 1; k <= l; k++)
g += A[k][j] * A[i][k];
e[j] = g / h;
f += e[j] * A[i][j];
}
double hh = f / (h + h);
for (int j = 0; j <= l; j++) {
f = A[i][j];
e[j] = g = e[j] - hh * f;
for (int k = 0; k <= j; k++)
A[j][k] -= f * e[k] + g * A[i][k];
}
}
} else
e[i] = A[i][l];
d[i] = h;
}
d[0] = 0.0;
e[0]=0.0;
for (int i = 0; i < n; i++) {
if (d[i] != 0.0) {
for (int j = 0; j < i; j++) {
double g = 0;
for (int k = 0; k < i; k++)
g += A[i][k] * A[k][j];
for (int k = 0; k < i; k++)
A[k][j] -= g * A[k][i];
}
}
d[i] = A[i][i];
A[i][i] = 1.0;
for (int j = 0; j < i; j++)
A[j][i] = A[i][j] = 0.0;
}
}
static double pythag(double a, double b) {
double absa = std::abs(a);
double absb = std::abs(b);
if (absa > absb) {
double ratio = absb / absa;
return absa * std::sqrt(1 + ratio * ratio);
} else {
if (absb == 0.0)
return 0.0;
else {
double ratio = absa / absb;
return absb * std::sqrt(1 + ratio * ratio);
}
}
}
inline double sign(double a, double b) {
if (b >= 0.0) {
if (a >= 0.0)
return a;
else
return -a;
} else {
if (a >= 0.0)
return -a;
else
return a;
}
}
void solveTQLI(Vector& d, Vector& e, Matrix& z) {
int n = d.dimension();
for (int i = 1; i < n; i++)
e[i-1] = e[i];
e[n-1] = 0.0;
for (int l = 0; l < n; l++) {
int iter = 0, m;
do {
for (m = l ; m < n-1; m++) {
double dd = std::abs(d[m]) + std::abs(d[m+1]);
if ((std::abs(e[m]) + dd) == dd)
break;
}
if (m != l) {
if (iter++ == 30)
error("Too many iterations in solveTQLI");
double g = (d[l+1] - d[l]) / (2.0 * e[l]);
double r = pythag(g, 1.0);
g = d[m] - d[l] + e[l] / (g + sign(r, g));
double s = 1.0;
double c = 1.0;
double p = 0.0;
int i;
for (i = m-1; i >= l; i--) {
double f = s * e[i];
double b = c * e[i];
e[i+1] = r = pythag(f, g);
if (r == 0.0) {
d[i+1] -= p;
e[m] = 0.0;
break;
}
s = f / r;
c = g / r;
g = d[i+1] - p;
r = (d[i] - g) * s + 2.0 * c * b;
d[i+1] = g + (p = s * r);
g = c * r - b;
for (int k = 0; k < n; k++) {
f = z[k][i+1];
z[k][i+1] = s * z[k][i] + c * f;
z[k][i] = c * z[k][i] - s * f;
}
}
if (r == 0.0 && i >= l)
continue;
d[l] -= p;
e[l] = g;
e[m] = 0.0;
}
} while (m != l);
}
}
void sortEigen(Vector& d, Matrix& z) {
// sorts eigenvalues and eigenvector in descending order
int n = d.dimension();
if (z.numRows() != n || z.numCols() != n)
error("Bad vector, matrix dimensions in sortEigen");
for (int i = 0; i < n - 1; i++) {
int k = i;
double p = d[k];
for (int j = i; j < n; j++)
if (d[j] >= p)
p = d[k = j];
if (k != i) {
d[k] = d[i];
d[i] = p;
for (int j = 0; j < n; j++) {
p = z[j][i];
z[j][i] = z[j][k];
z[j][k] = p;
}
}
}
}
void singularValueDecompose(Matrix& a, Vector& w, Matrix& v) {
bool flag;
int i,its,j,jj,k,l,nm;
double anorm,c,f,g,h,s,scale,x,y,z;
int m=a.numRows();
int n=a.numCols();
Vector rv1(n);
g=scale=anorm=0.0;
for (i=0;i<n;i++) {
l=i+2;
rv1[i]=scale*g;
g=s=scale=0.0;
if (i < m) {
for (k=i;k<m;k++) scale += std::abs(a[k][i]);
if (scale != 0.0) {
for (k=i;k<m;k++) {
a[k][i] /= scale;
s += a[k][i]*a[k][i];
}
f=a[i][i];
g = -sign(std::sqrt(s),f);
h=f*g-s;
a[i][i]=f-g;
for (j=l-1;j<n;j++) {
for (s=0.0,k=i;k<m;k++) s += a[k][i]*a[k][j];
f=s/h;
for (k=i;k<m;k++) a[k][j] += f*a[k][i];
}
for (k=i;k<m;k++) a[k][i] *= scale;
}
}
w[i]=scale *g;
g=s=scale=0.0;
if (i+1 <= m && i != n) {
for (k=l-1;k<n;k++) scale += std::abs(a[i][k]);
if (scale != 0.0) {
for (k=l-1;k<n;k++) {
a[i][k] /= scale;
s += a[i][k]*a[i][k];
}
f=a[i][l-1];
g = -sign(std::sqrt(s),f);
h=f*g-s;
a[i][l-1]=f-g;
for (k=l-1;k<n;k++) rv1[k]=a[i][k]/h;
for (j=l-1;j<m;j++) {
for (s=0.0,k=l-1;k<n;k++) s += a[j][k]*a[i][k];
for (k=l-1;k<n;k++) a[j][k] += s*rv1[k];
}
for (k=l-1;k<n;k++) a[i][k] *= scale;
}
}
anorm=std::max(anorm,(std::abs(w[i])+std::abs(rv1[i])));
}
for (i=n-1;i>=0;i--) {
if (i < n-1) {
if (g != 0.0) {
for (j=l;j<n;j++)
v[j][i]=(a[i][j]/a[i][l])/g;
for (j=l;j<n;j++) {
for (s=0.0,k=l;k<n;k++) s += a[i][k]*v[k][j];
for (k=l;k<n;k++) v[k][j] += s*v[k][i];
}
}
for (j=l;j<n;j++) v[i][j]=v[j][i]=0.0;
}
v[i][i]=1.0;
g=rv1[i];
l=i;
}
for (i=std::min(m,n)-1;i>=0;i--) {
l=i+1;
g=w[i];
for (j=l;j<n;j++) a[i][j]=0.0;
if (g != 0.0) {
g=1.0/g;
for (j=l;j<n;j++) {
for (s=0.0,k=l;k<m;k++) s += a[k][i]*a[k][j];
f=(s/a[i][i])*g;
for (k=i;k<m;k++) a[k][j] += f*a[k][i];
}
for (j=i;j<m;j++) a[j][i] *= g;
} else for (j=i;j<m;j++) a[j][i]=0.0;
++a[i][i];
}
for (k=n-1;k>=0;k--) {
for (its=0;its<30;its++) {
flag=true;
for (l=k;l>=0;l--) {
nm=l-1;
if ((std::abs(rv1[l])+anorm) == anorm) {
flag=false;
break;
}
if ((std::abs(w[nm])+anorm) == anorm) break;
}
if (flag) {
c=0.0;
s=1.0;
for (i=l-1;i<k+1;i++) {
f=s*rv1[i];
rv1[i]=c*rv1[i];
if ((std::abs(f)+anorm) == anorm) break;
g=w[i];
h=pythag(f,g);
w[i]=h;
h=1.0/h;
c=g*h;
s = -f*h;
for (j=0;j<m;j++) {
y=a[j][nm];
z=a[j][i];
a[j][nm]=y*c+z*s;
a[j][i]=z*c-y*s;
}
}
}
z=w[k];
if (l == k) {
if (z < 0.0) {
w[k] = -z;
for (j=0;j<n;j++) v[j][k] = -v[j][k];
}
break;
}
if (its == 30) error("singularValueDecomposition: no convergence in 30 iterations");
x=w[l];
nm=k-1;
y=w[nm];
g=rv1[nm];
h=rv1[k];
f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
g=pythag(f,1.0);
f=((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x;
c=s=1.0;
for (j=l;j<=nm;j++) {
i=j+1;
g=rv1[i];
y=w[i];
h=s*g;
g=c*g;
z=pythag(f,h);
rv1[j]=z;
c=f/z;
s=h/z;
f=x*c+g*s;
g = g*c-x*s;
h=y*s;
y *= c;
for (jj=0;jj<n;jj++) {
x=v[jj][j];
z=v[jj][i];
v[jj][j]=x*c+z*s;
v[jj][i]=z*c-x*s;
}
z=pythag(f,h);
w[j]=z;
if (z != 0.0) {
z=1.0/z;
c=f*z;
s=h*z;
}
f=c*g+s*y;
x=c*y-s*g;
for (jj=0;jj<m;jj++) {
y=a[jj][j];
z=a[jj][i];
a[jj][j]=y*c+z*s;
a[jj][i]=z*c-y*s;
}
}
rv1[l]=0.0;
rv1[k]=f;
w[k]=x;
}
}
}
Vector solveEigenSymmetric(Matrix& A) {
int n = A.numRows();
Vector d(n), e(n);
reduceHouseholder(A, d, e);
solveTQLI(d, e, A);
sortEigen(d, A);
return d;
}
void lineSearch(Vector& xold, double fold, Vector& g, Vector& p, Vector& x,
double& f, double stpmax, bool& check, double (*func)(Vector&))
{
const double ALF=1.0e-4, TOLX=std::numeric_limits<double>::epsilon();
int i;
double a,alam,alam2=0.0,alamin,b,disc,f2=0.0;
double rhs1,rhs2,slope,sum,temp,test,tmplam;
int n = xold.dimension();
check = false;
sum=0.0;
for (i=0;i<n;i++) sum += p[i]*p[i];
sum=std::sqrt(sum);
if (sum > stpmax)
for (i=0;i<n;i++) p[i] *= stpmax/sum;
slope = 0.0;
for (i=0;i<n;i++)
slope += g[i]*p[i];
if (slope >= 0.0) error("Roundoff problem in lineSearch");
test=0.0;
for (i=0;i<n;i++) {
temp=std::abs(p[i])/std::max(std::abs(xold[i]),1.0);
if (temp > test) test=temp;
}
alamin=TOLX/test;
alam=1.0;
for (;;) {
for (i=0;i<n;i++) x[i]=xold[i]+alam*p[i];
f=func(x);
if (alam < alamin) {
for (i=0;i<n;i++) x[i]=xold[i];
check = true;
return;
} else if (f <= fold+ALF*alam*slope) return;
else {
if (alam == 1.0)
tmplam = -slope/(2.0*(f-fold-slope));
else {
rhs1 = f-fold-alam*slope;
rhs2=f2-fold-alam2*slope;
a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
if (a == 0.0) tmplam = -slope/(2.0*b);
else {
disc=b*b-3.0*a*slope;
if (disc < 0.0) tmplam=0.5*alam;
else if (b <= 0.0) tmplam=(-b+std::sqrt(disc))/(3.0*a);
else tmplam=-slope/(b+std::sqrt(disc));
}
if (tmplam > 0.5*alam)
tmplam=0.5*alam;
}
}
alam2=alam;
f2 = f;
alam=std::max(tmplam,0.1*alam);
}
}
void minimizeBFGS(Vector& p, const double gtol, int& iter, double& fret,
double (*func)(Vector&), void (*dfunc)(Vector&, Vector&)) {
const int ITMAX = 200;
const double EPS = std::numeric_limits<double>::epsilon();
const double TOLX = 4*EPS, STPMX = 100.0;
bool check;
int i,its,j;
double den,fac,fad,fae,fp,stpmax,sum=0.0,sumdg,sumxi,temp,test;
int n = p.dimension();
Vector dg(n),g(n),hdg(n),pnew(n),xi(n);
Matrix hessin(n, n);
fp=func(p);
dfunc(p,g);
for (i=0;i<n;i++) {
for (j=0;j<n;j++) hessin[i][j]=0.0;
hessin[i][i]=1.0;
xi[i] = -g[i];
sum += p[i]*p[i];
}
stpmax=STPMX*std::max(std::sqrt(sum),double(n));
for (its=0;its<ITMAX;its++) {
iter=its;
lineSearch(p,fp,g,xi,pnew,fret,stpmax,check,func);
fp = fret;
for (i=0;i<n;i++) {
xi[i]=pnew[i]-p[i];
p[i]=pnew[i];
}
test=0.0;
for (i=0;i<n;i++) {
temp=std::abs(xi[i])/std::max(std::abs(p[i]),1.0);
if (temp > test) test=temp;
}
if (test < TOLX)
return;
for (i=0;i<n;i++) dg[i]=g[i];
dfunc(p,g);
test=0.0;
den=std::max(fret,1.0);
for (i=0;i<n;i++) {
temp=std::abs(g[i])*std::max(std::abs(p[i]),1.0)/den;
if (temp > test) test=temp;
}
if (test < gtol)
return;
for (i=0;i<n;i++) dg[i]=g[i]-dg[i];
for (i=0;i<n;i++) {
hdg[i]=0.0;
for (j=0;j<n;j++) hdg[i] += hessin[i][j]*dg[j];
}
fac=fae=sumdg=sumxi=0.0;
for (i=0;i<n;i++) {
fac += dg[i]*xi[i];
fae += dg[i]*hdg[i];
sumdg += dg[i]*dg[i];
sumxi += xi[i]*xi[i];
}
if (fac > std::sqrt(EPS*sumdg*sumxi)) {
fac=1.0/fac;
fad=1.0/fae;
for (i=0;i<n;i++) dg[i]=fac*xi[i]-fad*hdg[i];
for (i=0;i<n;i++) {
for (j=i;j<n;j++) {
hessin[i][j] += fac*xi[i]*xi[j]
-fad*hdg[i]*hdg[j]+fae*dg[i]*dg[j];
hessin[j][i]=hessin[i][j];
}
}
}
for (i=0;i<n;i++) {
xi[i]=0.0;
for (j=0;j<n;j++) xi[i] -= hessin[i][j]*g[j];
}
}
error("Too many iterations in minimizeBFGS");
}
} /* end namespace cpl */

View File

@ -1,86 +0,0 @@
#ifndef CPL_MATRIX_HPP
#define CPL_MATRIX_HPP
namespace cpl {
class Vector;
class Matrix {
public:
Matrix(int rows=1, int cols=1, double d=0);
Matrix(const Matrix& m);
~Matrix();
int numRows() const { return rows; }
int numCols() const { return cols; }
const double* operator[](const int row) const;
double* operator[](const int row);
double& operator()(int row, int col);
Matrix& operator=(const Matrix& m);
Matrix& operator=(const double d);
Matrix& operator += (const Matrix& m);
Matrix& operator -= (const Matrix& m);
Matrix& operator *= (double d);
Matrix& operator /= (double d);
Matrix transpose();
friend std::ostream& operator<<(std::ostream& os, const Matrix& mat);
private:
int rows;
int cols;
double **m;
};
inline Matrix operator + (const Matrix& m) {
return m;
}
extern Matrix operator - (const Matrix& m);
extern Matrix operator * (const Matrix& m, double d);
extern Matrix operator * (double d, const Matrix& m);
extern Matrix operator / (const Matrix& m, double d);
extern Matrix operator + (const Matrix& m1, const Matrix& m2);
extern Matrix operator - (const Matrix& m1, const Matrix& m2);
extern Matrix operator * (const Matrix& m1, const Matrix& m2);
extern void solveGaussJordan(Matrix& A, Matrix& B);
extern void solveLUDecompose(Matrix& A, Matrix& B);
extern void reduceHouseholder(Matrix& A, Vector& d, Vector& e);
extern void solveTQLI(Vector& d, Vector&e, Matrix& z);
extern void sortEigen(Vector& d, Matrix& z);
extern Vector solveEigenSymmetric(Matrix& A);
extern void singularValueDecompose(Matrix& a, Vector& w, Matrix& v);
extern void minimizeBFGS(Vector& p, const double gtol, int& iter, double& fret,
double (*func)(Vector&), void (*dfunc)(Vector&, Vector&));
} /* end namespace cpl */
#endif /* CPL_MATRIX_HPP */

View File

@ -1,136 +0,0 @@
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <iostream>
using namespace std;
#include "vector.hpp"
#include "odeint.hpp"
namespace cpl {
// Fourth order Runge-Kutta
void RK4Step(Vector& x, double tau,
Vector derivs(const Vector&))
{
Vector k1 = tau * derivs(x);
Vector k2 = tau * derivs(x + 0.5 * k1);
Vector k3 = tau * derivs(x + 0.5 * k2);
Vector k4 = tau * derivs(x + k3);
x += (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;
}
// Adaptive step size control using Runge-Kutta and step doubling
void adaptiveRK4Step(Vector& x, double& tau, double accuracy,
Vector derivs(const Vector&))
{
const double SAFETY = 0.9, PGROW = -0.2, PSHRINK = -0.25,
ERRCON = 1.89E-4, TINY = 1.0e-30;
int n = x.dimension();
Vector x_half(n), x_full(n), Delta(n);
Vector scale = derivs(x);
for (int i = 0; i < n; i++)
scale[i] = abs(x[i]) + abs(scale[i] * tau) + TINY;
double err_max;
while (true) {
// take two half steps
double tau_half = tau / 2;
x_half = x;
RK4Step(x_half, tau_half, derivs);
RK4Step(x_half, tau_half, derivs);
// take full step
x_full = x;
RK4Step(x_full, tau, derivs);
// estimate error
Delta = x_half - x_full;
err_max = 0;
for (int i = 0; i < n; i++)
err_max = max(err_max, abs(Delta[i]) / scale[i]);
err_max /= accuracy;
if (err_max <= 1.0)
break;
double tau_temp = SAFETY * tau * pow(err_max, PSHRINK);
if (tau >= 0.0)
tau = max(tau_temp, 0.1 * tau);
else
tau = min(tau_temp, 0.1 * tau);
if (abs(tau) == 0.0) {
cerr << "adaptiveRK4Step: step size underflow\naborting ..."
<< endl;
exit(EXIT_FAILURE);
}
}
tau *= (err_max > ERRCON ? SAFETY * pow(err_max, PGROW) : 5.0);
x = x_half + Delta / 15.0;
}
// Runge-Kutta-Cash-Karp including error estimate
static void rkck(Vector& x, double tau,
Vector derivs(const Vector&), Vector& x_err)
{
const double b21 = 1.0/5.0, b31 = 3.0/40.0, b41 = 3.0/10.0,
b51 = -11.0/54.0, b61 = 1631.0/55296.0, b32 = 9.0/40.0,
b42 = -9.0/10.0, b52 = 5.0/2.0, b62 = 175.0/512.0, b43 = 6.0/5.0,
b53 = -70.0/27.0, b63 = 575.0/13824.0, b54 = 35.0/27.0,
b64 = 44275.0/110592.0, b65 = 253.0/4096.0, c1 = 37.0/378.0,
c2 = 0.0, c3 = 250.0/621.0, c4 = 125.0/594.0, c5 = 0.0,
c6 = 512.0/1771.0, dc1 = c1 - 2825.0/27648.0, dc2 = c2 - 0.0,
dc3 = c3 - 18575.0/48384.0, dc4 = c4 - 13525.0/55296.0,
dc5 = c5 - 277.0/14336.0, dc6 = c6 - 1.0/4.0;
Vector
k1 = tau * derivs(x),
k2 = tau * derivs(x + b21*k1),
k3 = tau * derivs(x + b31*k1 + b32*k2),
k4 = tau * derivs(x + b41*k1 + b42*k2 + b43*k3),
k5 = tau * derivs(x + b51*k1 + b52*k2 + b53*k3 + b54*k4),
k6 = tau * derivs(x + b61*k1 + b62*k2 + b63*k3 + b64*k4 + b65*k5);
x += c1*k1 + c2*k2 + c3*k3 + c4*k4 + c5*k5 + c6*k6;
x_err = dc1*k1 + dc2*k2 + dc3*k3 + dc4*k4 + dc5*k5 + dc6*k6;
}
// Runge-Kutta-Cash-Karp step
void RKCKStep(Vector& x, double tau, Vector derivs(const Vector&))
{
Vector x_err(x.dimension());
rkck(x, tau, derivs, x_err);
}
// Adaptive step size control using Runge-Kutta-Cash-Karp
void adaptiveRKCKStep(Vector& x, double& tau, double accuracy,
Vector derivs(const Vector&))
{
const double SAFETY = 0.9, PGROW = -0.2, PSHRINK = -0.25,
ERRCON = 1.89E-4, TINY = 1.0e-30;
int n = x.dimension();
Vector x_err(n), x_temp(n);
Vector scale = derivs(x);
for (int i = 0; i < n; i++)
scale[i] = abs(x[i]) + abs(scale[i] * tau) + TINY;
double err_max;
while (true) {
// take Cash-Karp step including error estimate
x_temp = x;
rkck(x_temp, tau, derivs, x_err);
err_max = 0;
for (int i = 0; i < n; i++)
err_max = max(err_max, abs(x_err[i]) / scale[i]);
err_max /= accuracy;
if (err_max <= 1.0)
break;
double tau_temp = SAFETY * tau * pow(err_max, PSHRINK);
if (tau >= 0.0)
tau = max(tau_temp, 0.1 * tau);
else
tau = min(tau_temp, 0.1 * tau);
if (abs(tau) == 0.0) {
cerr << "adaptiveRKCKStep: step size underflow\naborting ..."
<< endl;
exit(EXIT_FAILURE);
}
}
tau *= (err_max > ERRCON ? SAFETY * pow(err_max, PGROW) : 5.0);
x = x_temp;
}
} /* end namespace cpl */

View File

@ -1,44 +0,0 @@
#ifndef CPL_ODEINT_HPP
#define CPL_ODEINT_HPP
#include <iostream>
namespace cpl {
class Vector;
// ODE integration routines adapted from Numerical Recipes
// take a single 4th order Runge-Kutta step
extern void RK4Step(
Vector& x, // extended solution vector
double tau, // step size
Vector derivs(const Vector& x) // extended derivative vector
);
// take one adaptive RK4 step using step doubling
extern void adaptiveRK4Step(
Vector& x, // extended solution vector
double& tau, // recommended step size
double accuracy, // desired solution accuracy
Vector derivs(const Vector&) // extended derivative vector
);
// take one Runge-Kutta-Cash-Karp step
extern void RKCKStep(
Vector& x, // extended solution vector
double tau, // step size
Vector derivs(const Vector& x) // derivative vector
);
// take one adaptive RKCK step using embedded error estimate
extern void adaptiveRKCKStep(
Vector& x, // extended solution vector
double& tau, // recommended step size
double accuracy, // desired solution accuracy
Vector derivs(const Vector&) // extended derivative vector
);
} /* end namespace cpl */
#endif /* CPL_ODEINT_HPP */

View File

@ -1,134 +0,0 @@
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <iostream>
using namespace std;
#include "optimize.hpp"
static int sign_of_func;
static double (*user_func)(double);
static double func(double x) {
return sign_of_func * (*user_func)(x);
}
inline void shft3(double &a, double &b, double &c, const double d) {
a=b; b=c; c=d;
}
inline double SIGN(const double& a, const double& b) {
return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
}
namespace cpl {
double Optimizer::findMaximum(const double a, const double b,
double f(double), double tol, double& xmin) {
sign_of_func = -1;
user_func = f;
double ax = a, bx = b, cx, fa, fb, fc;
mnbrak(ax, bx, cx, fa, fb, fc, func);
return -golden(ax, bx, cx, func, tol, xmin);
}
double Optimizer::findMinimum(const double a, const double b,
double f(double), double tol, double& xmin) {
double ax = a, bx = b, cx, fa, fb, fc;
mnbrak(ax, bx, cx, fa, fb, fc, f);
return golden(ax, bx, cx, f, tol, xmin);
}
void Optimizer::mnbrak(double &ax, double &bx, double &cx,
double &fa, double &fb, double &fc,
double func(const double)) {
const double GOLD=1.618034,GLIMIT=100.0,TINY=1.0e-20;
double ulim,u,r,q,fu;
fa=func(ax);
fb=func(bx);
if (fb > fa) {
swap(ax,bx);
swap(fb,fa);
}
cx=bx+GOLD*(bx-ax);
fc=func(cx);
while (fb > fc) {
r=(bx-ax)*(fb-fc);
q=(bx-cx)*(fb-fa);
u=bx-((bx-cx)*q-(bx-ax)*r)/
(2.0*SIGN(max(abs(q-r),TINY),q-r));
ulim=bx+GLIMIT*(cx-bx);
if ((bx-u)*(u-cx) > 0.0) {
fu=func(u);
if (fu < fc) {
ax=bx;
bx=u;
fa=fb;
fb=fu;
return;
} else if (fu > fb) {
cx=u;
fc=fu;
return;
}
u=cx+GOLD*(cx-bx);
fu=func(u);
} else if ((cx-u)*(u-ulim) > 0.0) {
fu=func(u);
if (fu < fc) {
shft3(bx,cx,u,cx+GOLD*(cx-bx));
shft3(fb,fc,fu,func(u));
}
} else if ((u-ulim)*(ulim-cx) >= 0.0) {
u=ulim;
fu=func(u);
} else {
u=cx+GOLD*(cx-bx);
fu=func(u);
}
shft3(ax,bx,cx,u);
shft3(fa,fb,fc,fu);
}
}
inline void shft2(double &a, double &b, const double c) {
a=b; b=c;
}
double Optimizer::golden(const double ax, const double bx, const double cx,
double f(const double), const double tol,
double &xmin) {
const double R=0.61803399,C=1.0-R;
double f1,f2,x0,x1,x2,x3;
x0=ax;
x3=cx;
if (abs(cx-bx) > abs(bx-ax)) {
x1=bx;
x2=bx+C*(cx-bx);
} else {
x2=bx;
x1=bx-C*(bx-ax);
}
f1=f(x1);
f2=f(x2);
while (abs(x3-x0) > tol*(abs(x1)+abs(x2))) {
if (f2 < f1) {
shft3(x0,x1,x2,R*x2+C*x3);
shft2(f1,f2,f(x2));
} else {
shft3(x3,x2,x1,R*x1+C*x0);
shft2(f2,f1,f(x1));
}
}
if (f1 < f2) {
xmin=x1;
return f1;
} else {
xmin=x2;
return f2;
}
}
} /* end namespace cpl */

View File

@ -1,47 +0,0 @@
#ifndef CPL_OPTIMIZE_HPP
#define CPL_OPTIMIZE_HPP
namespace cpl {
class Optimizer {
public:
double findMaximum( // returns f(x) at local maximum
const double xa, // guess for one point near desried maximum
const double xb, // guess for second point near maximum
double f(double), // name of function to be maximized
double accuracy, // desired accuracy in x
double& xMax); // value of x at local maximum
double findMinimum( // returns f(x) at local minimum
const double xa, // guess for one point near desried minimum
const double xb, // guess for second point near minimum
double f(double), // name of function to be minimized
double accuracy, // desired accuracy in x
double& xMin); // value of x at local minimum
double golden( // numerical recipes routine
const double,
const double,
const double,
double f(double),
const double,
double&
);
void mnbrak( // numerical recipes bracketing routine
double&,
double&,
double&,
double&,
double&,
double&,
double f(double)
);
};
} /* end namespace cpl */
#endif /* CPL_OPTIMIZE_HPP */

View File

@ -1,73 +0,0 @@
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <iostream>
#include <vector>
using namespace std;
#include "random.hpp"
namespace cpl {
int timeSeed() {
time_t t; // usually an unsigned long int
time(&t); // get seconds since Jan 1, 1900
tm* pt = gmtime(&t); // Universal (Greenwich Mean) Time
return pt->tm_year + 70 * (pt->tm_mon + 12 * (pt->tm_hour +
23 * (pt->tm_min + 59 * pt->tm_sec)));
}
double rand(int& seed, bool set) {
if (set)
srand(seed);
else
seed = std::rand();
return seed / (RAND_MAX + 1.0);
}
double ranf(int& seed, bool set) {
const int IA = 16807, IC = 2147483647, IQ = 127773, IR = 2836;
if (!set) {
int h = seed / IQ;
int l = seed % IQ;
seed = IA * l - IR * h;
}
if (seed <= 0)
seed += IC;
return seed / double(IC);
}
double qand(int& seed, bool set) {
static unsigned long idum;
const double TWO_POWER_32 = 4294967296.0;
if (set) {
idum = (unsigned long) seed;
return idum / TWO_POWER_32;
}
idum = 1664525L * idum + 1013904223L;
seed = int(idum);
return idum / TWO_POWER_32;
}
double gasdev(int& seed) {
static int iset = 0;
static double gset;
double fac, rsq, v1, v2;
if (iset == 0) {
do {
v1 = 2.0 * ranf(seed) - 1.0;
v2 = 2.0 * ranf(seed) - 1.0;
rsq = v1 * v1 + v2 * v2;
} while (rsq >= 1.0 || rsq == 0.0);
fac = sqrt(-2.0 * log(rsq)/rsq);
gset = v1 * fac;
iset = 1;
return v2 * fac;
} else {
iset = 0;
return gset;
}
}
} /* end namespace cpl */

View File

@ -1,23 +0,0 @@
#ifndef CPL_RANDOM_HPP
#define CPL_RANDOM_HPP
namespace cpl {
// Random integer seed based on the current time of day
extern int timeSeed();
// Standard C++ Library generator
extern double rand(int& seed, bool set=false);
// Park-Miller generator with Schrage's algorithm to prevent overflow
extern double ranf(int& seed, bool set=false);
// Quick and dirty generator according to Numerical Recipes
extern double qand(int& seed, bool set=false);
// Gaussian distributed random number based on Box-Muller algorithm
extern double gasdev(int& seed);
} /* end namespace cpl */
#endif /* CPL_RANDOM_HPP */

View File

@ -1,167 +0,0 @@
#include <algorithm>
#include <cmath>
#include <cstdarg>
#include <cstdlib>
#include <iostream>
using namespace std;
#include "vector.hpp"
namespace cpl {
Vector::Vector(int dim) {
v = new double [this->dim = dim];
for (int i = 0; i < dim; i++) v[i] = 0;
}
Vector::Vector(double v0, double v1) {
v = new double [dim = 2];
v[0] = v0;
v[1] = v1;
}
Vector::Vector(double v0, double v1, double v2) {
v = new double [dim = 3];
v[0] = v0;
v[1] = v1;
v[2] = v2;
}
Vector::Vector(const Vector& dv) {
v = new double [dim = dv.dim];
for (int i = 0; i < dim; i++) v[i] = dv.v[i];
}
void Vector::resize(const int dimension) {
double *old = new double[dim];
for (int i = 0; i < dim; i++) old[i] = v[i];
delete [] v;
v = new double [dimension];
for (int i = 0; i < dimension; i++)
v[i] = i < dim ? old[i] : 0;
dim = dimension;
delete [] old;
}
void Vector::push_back(const double d) {
resize(++dim);
v[dim-1] = d;
}
void Vector::set(const double d0, ...) {
va_list args;
v[0] = d0;
va_start(args, d0);
for (int i = 1; i < dim; i++)
v[i] = va_arg(args, double);
va_end(args);
}
Vector& Vector::operator = (const Vector& dv) {
if (this != &dv) {
if (dim != dv.dim) {
delete [] v;
v = new double [dim = dv.dim];
}
for (int i = 0; i < dim; i++) v[i] = dv[i];
}
return *this;
}
Vector& Vector::operator += (const Vector& dv) {
for (int i = 0; i < dim; i++) v[i] += dv[i];
return *this;
}
Vector& Vector::operator -= (const Vector& dv) {
for (int i = 0; i < dim; i++) v[i] -= dv[i];
return *this;
}
Vector& Vector::operator *= (double d) {
for (int i = 0; i < dim; i++) v[i] *= d;
return *this;
}
Vector& Vector::operator /= (double d) {
for (int i = 0; i < dim; i++) v[i] /= d;
return *this;
}
Vector operator - (const Vector& dv) {
int dim = dv.dimension();
Vector temp(dim);
for (int i = 0; i < dim; i++) temp[i] = -dv[i];
return temp;
}
Vector operator * (const Vector& dv, double d) {
int dim = dv.dimension();
Vector temp(dim);
for (int i = 0; i < dim; i++) temp[i] = dv[i] * d;
return temp;
}
Vector operator * (double d, const Vector& dv) {
int dim = dv.dimension();
Vector temp(dim);
for (int i = 0; i < dim; i++) temp[i] = dv[i] * d;
return temp;
}
Vector operator / (const Vector& dv, double d) {
int dim = dv.dimension();
Vector temp(dim);
for (int i = 0; i < dim; i++) temp[i] = dv[i] / d;
return temp;
}
Vector operator + (const Vector& v1, const Vector& v2) {
int dim = v1.dimension();
Vector temp(dim);
for (int i = 0; i < dim; i++) temp[i] = v1[i] + v2[i];
return temp;
}
Vector operator - (const Vector& v1, const Vector& v2) {
int dim = v1.dimension();
Vector temp(dim);
for (int i = 0; i < dim; i++) temp[i] = v1[i] - v2[i];
return temp;
}
double Vector::abs() {
return std::sqrt(norm());
}
double Vector::norm() {
double sum = 0;
for (int i = 0; i < dim; i++) sum += v[i] * v[i];
return sum;
}
double Vector::dot(const Vector& dv) {
double sum = 0;
for (int i = 0; i < dim; i++) sum += v[i] * dv[i];
return sum;
}
std::ostream& operator<<(std::ostream& os, const Vector& dv) {
for (int i = 0; i < dv.dim; i++) {
os << dv.v[i];
if (i < dv.dim-1)
os << '\t';
else
os << '\n';
}
return os;
}
double dot(const Vector& v1, const Vector& v2) {
double sum = 0;
for (int i = 0; i < v1.size(); i++)
sum += v1[i] * v2[i];
return sum;
}
} /* end namespace cpl */

View File

@ -1,124 +0,0 @@
#ifndef CPL_VECTOR_HPP
#define CPL_VECTOR_HPP
#include <iostream>
namespace cpl {
// Vector object with components of type double
class Vector {
public:
Vector( // constructor, creates a new vector
int dim = 1 // dimension, assumed = 1 if omitted
); // components initialized to 0
Vector( // constructor, creates a new 2-d vector
double v0, // value of first component
double v1 // value of second componet
);
Vector( // constructor, creates a new 3-d vector
double v0, // value of first component
double v1, // value of second component
double v2 // value of third component
);
Vector( // copy constructor, creates a new vector
const Vector& // vector to be copied
);
~Vector() { // destructor frees memory allocated by constructors
delete [] v; // inline code to delete array v
}
int dimension() const // returns number of components
{ return dim; } // inline code to get dimension
int size() const // same as dimension()
{ return dim; }
void resize( // resize the vector
const int // new number of components
); // old components preserved, any new zeroed
void push_back( // add a new component to the vector
const double // value to set the new component
);
void set( // set all the components of the vector
const double, // value of first component
... // second, third, ... , last component
); // all components must be specified
const double // value returned by
operator[]( // [ ] operator function
const int i // index of component to return
) const // value cannot be modified
{ return v[i]; } // inline code to get value
double& // reference returned by
operator[]( // [ ] operator function
const int i // index of component to return
) { return v[i]; } // inline code to get component
// operators to allow the vector to be assigned
Vector& operator = (const Vector&);
Vector& operator += (const Vector&);
Vector& operator -= (const Vector&);
// multiply or divide the vector by a scalar
Vector& operator *= (double);
Vector& operator /= (double);
double abs(); // return the length of the vector
double norm(); // return the norm of the vector
double dot( // return dot product of the vector
const Vector& // with this vector
);
// allow the << operator to be used to output vector components
friend std::ostream& operator<<(std::ostream& os, const Vector& dv);
private:
int dim; // store the number of components of the vector
double *v; // pointer to stored components values
};
// various arithmetic operations on vectors
// unary plus using inline code, just returns the same vector
inline Vector operator + (const Vector& vec) { return vec; }
// unary minus, flips the sign of each component
extern Vector operator - (const Vector&);
// multiplying and dividing a vector by a scalar
extern Vector operator * (const Vector&v, double);
extern Vector operator * (double, const Vector&);
extern Vector operator / (const Vector&, double);
// adding and subtracting two vectors
extern Vector operator + (const Vector&, const Vector&);
extern Vector operator - (const Vector&, const Vector&);
// dot product of two vectors
extern double dot(const Vector&, const Vector&);
} /* end namespace cpl */
#endif /* CPL_VECTOR_HPP */

View File

@ -1,30 +0,0 @@
PROG = sho
SRCS = sho.cpp
OBJS = ${SRCS:.cpp=.o}
HDRS =
DEPS = popt
CHEZ = /usr/lib/csv10.0.0/ta6le
CPPFLAGS += -D_FORTIFY_SOURCE=3 -DCHEZDIR=\"${CHEZ}\"
CXXFLAGS += -std=c++11 -Wall -Wextra -I${CHEZ} `pkg-config --cflags ${DEPS}`
LDLIBS = ${LDFLAGS} -lm `pkg-config --libs ${DEPS}`
all: ${PROG}
${OBJS}: Makefile
${PROG}: ${OBJS}
${CXX} -o $@ ${OBJS} ${LDLIBS} # ${CHEZ}/kernel.o
${SRCS}: ${HDRS}
lint:
clang-tidy ${SRCS} ${HDRS} -- ${CPPFLAGS} ${CXXFLAGS}
plot: plot1.png
plot1.png: test.txt plot1.gnuplot
gnuplot plot1.gnuplot
.PHONY: all lint

View File

@ -1,46 +0,0 @@
#include <scheme.h>
#include "actkbd.h"
/* is the device currently grabbed? */
static int grabbed_p(void) {
/* this value comes from C, so an extra assert won't hurt */
assert((grabbed == 0) || (grabbed == 1));
return grabbed;
}
void extension_init(void) {
FILE *stream;
assert(config != NULL);
Sscheme_init(0);
Sregister_boot_file(CHEZDIR "/petite.boot");
Sregister_boot_file(CHEZDIR "/scheme.boot");
Sbuild_heap(0, 0);
Sforeign_symbol("bind", (ptr) &set_cb);
Sforeign_symbol("grab", (ptr) &set_grab);
Sforeign_symbol("grabbed?", (ptr) &grabbed_p);
Sforeign_symbol("led", (ptr) &set_led);
/* test whether we can read the config */
stream = fopen(config, "r");
if (!stream) {
fprintf(stderr, "ERROR : %s not found or reading is not allowed.\n",
config);
exit(CONFERR);
}
fclose(stream);
Sscheme_program(config, 0, NULL);
}
/* a shortcut to reload stuff */
void extension_reload(void) {
Sscheme_program(config, 0, NULL);
}
void extension_cleanup_callback(CB fun) {
/* so the gc will be able to collect it */
Sunlock_object(fun);
}

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -1,19 +0,0 @@
set encoding utf8
set terminal png size 800,600 enhanced font 'IBM Plex Sans Text,14'
set output "plot1.png"
set xlabel 'idő [ms]'
set ylabel 'kitérés [m]'
set style line 100 lt 1 lc rgb "grey" lw 0.5
set grid ls 100
set style line 103 lw 4 lt rgb "#b8ee30"
set style line 104 pt 1 ps 1 lw 1 lt rgb "#8b1c62"
set title 'Harmonikus oszcillátor'
plot 'test.txt' title 'kitérés' with points ls 104

View File

@ -1,77 +0,0 @@
#include <cmath>
#include <iostream>
#include <functional>
#include <popt.h>
using namespace std;
double omega; // the natural frequency
double x, v; // position and velocity at time t
int periods; // number of periods to integrate
int steps; // number of time steps dt per period
void getInput(void); // for user to input parameters
void EulerCromer(double dt); // takes an Euler-Cromer step
void Euler(double dt); // takes an Euler step
double energy(void); // computes the energy
void sim(function<void(double)> integrator); // runs the simulation
typedef double v2d __attribute__ ((vector_size (2*8)));
enum {
OK=0,
USAGE=1,
HOSTFAIL=2,
};
void EulerCromer (double dt) {
double a = - omega * omega * x;
v += a * dt;
x += v * dt;
}
void Euler(double dt) {
double a = - omega * omega * x;
x += v * dt;
v += a * dt;
}
double energy(void) {
return 0.5 * (v * v + omega * omega * x * x);
}
void sim(function<void(double)> integrator) {
const double pi = 4 * atan(1.0);
double T = 2 * pi / omega;
double dt = T / steps;
double t = 0;
cout << t << '\t' << x << '\t' << v << '\t' << energy() << '\n';
for (int p = 1; p <= periods; p++) {
for (int s = 0; s < steps; s++) {
integrator(dt);
t += dt;
cout << t << '\t' << x << '\t' << v << '\t' << energy() << '\n';
}
// cout << "Period = " << p << "\tt = " << t
// cout << "\tx = " << x << "\tv = " << v
// cout << "\tenergy = " << energy() << endl;
}
}
int main (int argc, const char *argv[]) {
const struct poptOption options_table[] = {
{"omega", 'o', POPT_ARG_DOUBLE, &omega, 0, "Value of omega", NULL},
{"xnull", 'x', POPT_ARG_DOUBLE, &x, 0, "Value of x(0)", NULL},
{"vnull", 'v', POPT_ARG_DOUBLE, &v, 0, "Value of v(0)", NULL},
{"periods", 'p', POPT_ARG_INT, &periods, 0, "Number of periods", NULL},
{"steps", 's', POPT_ARG_INT, &steps, 0, "Steps per period", NULL},
POPT_AUTOHELP
{NULL, 0, 0, NULL, 0, NULL, NULL}};
poptContext opt_con = poptGetContext("sho", argc, argv, options_table, 0);
poptSetOtherOptionHelp(opt_con, "[OPTIONS]*");
poptGetNextOpt(opt_con);
sim(Euler);
return EXIT_SUCCESS;
}

Binary file not shown.

View File

@ -1,201 +0,0 @@
0 1 0.1 50.005
0.00628319 0.99668 -0.528319 49.8082
0.0125664 0.989426 -1.15455 49.6147
0.0188496 0.978266 -1.77623 49.4277
0.0251327 0.963243 -2.39089 49.2501
0.0314159 0.944418 -2.99611 49.0846
0.0376991 0.921865 -3.58951 48.934
0.0439823 0.895672 -4.16873 48.8006
0.0502655 0.865943 -4.7315 48.6864
0.0565487 0.832795 -5.27559 48.5933
0.0628319 0.79636 -5.79885 48.5228
0.069115 0.756781 -6.29922 48.4759
0.0753982 0.714214 -6.77472 48.4535
0.0816814 0.668828 -7.22347 48.4558
0.0879646 0.620801 -7.64371 48.4828
0.0942478 0.570323 -8.03377 48.5342
0.100531 0.517594 -8.39211 48.609
0.106814 0.462822 -8.71733 48.7061
0.113097 0.406222 -9.00813 48.824
0.119381 0.348018 -9.26336 48.9608
0.125664 0.288441 -9.48203 49.1144
0.131947 0.227725 -9.66326 49.2823
0.13823 0.16611 -9.80635 49.4618
0.144513 0.103839 -9.91072 49.6503
0.150796 0.0411583 -9.97596 49.8446
0.15708 -0.021685 -10.0018 50.0417
0.163363 -0.0844427 -9.9882 50.2386
0.169646 -0.146867 -9.93514 50.432
0.175929 -0.208712 -9.84286 50.619
0.182212 -0.269732 -9.71172 50.7965
0.188496 -0.329688 -9.54224 50.9619
0.194779 -0.388342 -9.3351 51.1125
0.201062 -0.445463 -9.09109 51.2458
0.207345 -0.500825 -8.8112 51.3599
0.213628 -0.554211 -8.49652 51.4529
0.219911 -0.605408 -8.1483 51.5233
0.226195 -0.654215 -7.76791 51.5701
0.232478 -0.70044 -7.35686 51.5925
0.238761 -0.743899 -6.91676 51.59
0.245044 -0.784421 -6.44935 51.5629
0.251327 -0.821847 -5.95649 51.5115
0.257611 -0.856028 -5.4401 51.4366
0.263894 -0.88683 -4.90225 51.3394
0.270177 -0.914131 -4.34503 51.2214
0.27646 -0.937822 -3.77067 51.0845
0.282743 -0.957812 -3.18142 50.9309
0.289027 -0.97402 -2.57961 50.7629
0.29531 -0.986383 -1.96761 50.5833
0.301593 -0.994852 -1.34785 50.3948
0.307876 -0.999393 -0.722766 50.2005
0.314159 -0.999989 -0.094829 50.0034
0.320442 -0.996637 0.533482 49.8065
0.326726 -0.98935 1.15969 49.6131
0.333009 -0.978158 1.78131 49.4262
0.339292 -0.963104 2.39591 49.2487
0.345575 -0.944248 3.00105 49.0833
0.351858 -0.921664 3.59433 48.9328
0.358142 -0.895442 4.17343 48.7996
0.364425 -0.865684 4.73606 48.6856
0.370708 -0.832509 5.27998 48.5927
0.376991 -0.796047 5.80306 48.5223
0.383274 -0.756443 6.30323 48.4757
0.389557 -0.713852 6.77852 48.4534
0.395841 -0.668443 7.22705 48.4559
0.402124 -0.620396 7.64704 48.4832
0.408407 -0.569899 8.03685 48.5347
0.41469 -0.517152 8.39493 48.6097
0.420973 -0.462363 8.71986 48.707
0.427257 -0.405749 9.01037 48.825
0.43354 -0.347534 9.26531 48.962
0.439823 -0.287946 9.48367 49.1157
0.446106 -0.227221 9.6646 49.2837
0.452389 -0.1656 9.80736 49.4634
0.458673 -0.103325 9.91141 49.6519
0.464956 -0.0406416 9.97633 49.8462
0.471239 0.022202 10.0019 50.0433
0.477522 0.0849579 9.98792 50.2402
0.483805 0.147379 9.93454 50.4336
0.490088 0.209217 9.84194 50.6205
0.496372 0.27023 9.71048 50.798
0.502655 0.330176 9.54069 50.9632
0.508938 0.388818 9.33324 51.1136
0.515221 0.445926 9.08894 51.2469
0.521504 0.501273 8.80875 51.3608
0.527788 0.554641 8.49379 51.4536
0.534071 0.605819 8.1453 51.5238
0.540354 0.654606 7.76465 51.5704
0.546637 0.700809 7.35335 51.5925
0.55292 0.744244 6.91302 51.5899
0.559203 0.784742 6.4454 51.5626
0.565487 0.822142 5.95233 51.511
0.57177 0.856296 5.43576 51.4359
0.578053 0.887069 4.89774 51.3385
0.584336 0.91434 4.34038 51.2204
0.590619 0.938002 3.76588 51.0833
0.596903 0.957961 3.17652 50.9296
0.603186 0.974137 2.57461 50.7615
0.609469 0.986469 1.96254 50.5818
0.615752 0.994905 1.34273 50.3933
0.622035 0.999414 0.717608 50.1989
0.628319 0.999977 0.0896581 50.0018
0.634602 0.996593 -0.538646 49.8049
0.640885 0.989274 -1.16482 49.6116
0.647168 0.97805 -1.7864 49.4247
0.653451 0.962964 -2.40093 49.2472
0.659734 0.944077 -3.00598 49.082
0.666018 0.921463 -3.59916 48.9317
0.672301 0.895211 -4.17813 48.7985
0.678584 0.865425 -4.74061 48.6847
0.684867 0.832222 -5.28437 48.592
0.69115 0.795734 -5.80727 48.5218
0.697434 0.756104 -6.30725 48.4754
0.703717 0.71349 -6.78232 48.4533
0.71 0.668059 -7.23062 48.456
0.716283 0.61999 -7.65037 48.4835
0.722566 0.569473 -8.03992 48.5352
0.728849 0.516709 -8.39774 48.6104
0.735133 0.461905 -8.72239 48.7079
0.741416 0.405277 -9.01262 48.8261
0.747699 0.347049 -9.26726 48.9632
0.753982 0.287451 -9.48532 49.117
0.760265 0.226718 -9.66593 49.2851
0.766549 0.16509 -9.80838 49.4649
0.772832 0.10281 -9.91211 49.6534
0.779115 0.0401249 -9.9767 49.8478
0.785398 -0.022719 -10.0019 50.045
0.791681 -0.0854732 -9.98764 50.2418
0.797965 -0.14789 -9.93394 50.4351
0.804248 -0.209723 -9.84102 50.622
0.810531 -0.270728 -9.70924 50.7994
0.816814 -0.330664 -9.53914 50.9645
0.823097 -0.389295 -9.33138 51.1148
0.82938 -0.446389 -9.08678 51.2479
0.835664 -0.50172 -8.8063 51.3616
0.841947 -0.555071 -8.49106 51.4543
0.84823 -0.606231 -8.1423 51.5243
0.854513 -0.654997 -7.76139 51.5707
0.860796 -0.701178 -7.34985 51.5926
0.86708 -0.74459 -6.90928 51.5898
0.873363 -0.785063 -6.44144 51.5623
0.879646 -0.822436 -5.94817 51.5104
0.885929 -0.856563 -5.43142 51.4352
0.892212 -0.887308 -4.89323 51.3376
0.898495 -0.91455 -4.33572 51.2193
0.904779 -0.938182 -3.76109 51.0821
0.911062 -0.958109 -3.17161 50.9282
0.917345 -0.974255 -2.56961 50.7601
0.923628 -0.986554 -1.95747 50.5803
0.929911 -0.994958 -1.3376 50.3917
0.936195 -0.999435 -0.71245 50.1973
0.942478 -0.999966 -0.0844871 50.0001
0.948761 -0.996549 0.54381 49.8033
0.955044 -0.989198 1.16996 49.61
0.961327 -0.977941 1.79149 49.4232
0.967611 -0.962824 2.40595 49.2458
0.973894 -0.943906 3.01091 49.0807
0.980177 -0.921262 3.60398 48.9305
0.98646 -0.89498 4.18283 48.7975
0.992743 -0.865166 4.74516 48.6839
0.999026 -0.831935 5.28876 48.5913
1.00531 -0.795421 5.81148 48.5214
1.01159 -0.755766 6.31126 48.4751
1.01788 -0.713127 6.78612 48.4533
1.02416 -0.667674 7.23419 48.4562
1.03044 -0.619584 7.6537 48.4838
1.03673 -0.569048 8.043 48.5357
1.04301 -0.516266 8.40054 48.6111
1.04929 -0.461446 8.72492 48.7088
1.05558 -0.404804 9.01486 48.8271
1.06186 -0.346564 9.2692 48.9644
1.06814 -0.286955 9.48696 49.1183
1.07442 -0.226214 9.66726 49.2866
1.08071 -0.16458 9.80939 49.4664
1.08699 -0.102296 9.9128 49.655
1.09327 -0.0396082 9.97707 49.8494
1.09956 0.023236 10.002 50.0466
1.10584 0.0859884 9.98736 50.2434
1.11212 0.148401 9.93333 50.4367
1.11841 0.210228 9.84009 50.6235
1.12469 0.271226 9.708 50.8008
1.13097 0.331152 9.53758 50.9658
1.13726 0.389771 9.32951 51.116
1.14354 0.446851 9.08461 51.2489
1.14982 0.502168 8.80385 51.3625
1.15611 0.555501 8.48833 51.4549
1.16239 0.606642 8.13929 51.5248
1.16867 0.655388 7.75813 51.571
1.17496 0.701546 7.34634 51.5927
1.18124 0.744935 6.90554 51.5897
1.18752 0.785383 6.43749 51.5619
1.19381 0.82273 5.94402 51.5099
1.20009 0.85683 5.42708 51.4345
1.20637 0.887546 4.88872 51.3367
1.21265 0.914759 4.33106 51.2182
1.21894 0.938361 3.7563 51.0809
1.22522 0.958258 3.16671 50.9269
1.2315 0.974372 2.56462 50.7586
1.23779 0.986639 1.9524 50.5788
1.24407 0.995011 1.33248 50.3901
1.25035 0.999455 0.707292 50.1957
1.25664 0.999954 0.0793161 49.9985

View File

@ -1,30 +0,0 @@
PROG = pen
SRCS = pendulum.cpp pendulum-gl.cpp
OBJS = ${SRCS:.cpp=.o}
HDRS =
DEPS = popt gl glu glut
CHEZ = /usr/lib/csv10.0.0/ta6le
CPPFLAGS += -D_FORTIFY_SOURCE=3 -DCHEZDIR=\"${CHEZ}\"
CXXFLAGS += -std=c++11 -Wall -Wextra -I${CHEZ} `pkg-config --cflags ${DEPS}`
LDLIBS = ${LDFLAGS} -lm `pkg-config --libs ${DEPS}`
all: ${PROG}
${OBJS}: Makefile
${PROG}: ${OBJS}
${CXX} -o $@ ${OBJS} ${LDLIBS} # ${CHEZ}/kernel.o
${SRCS}: ${HDRS}
lint:
clang-tidy ${SRCS} ${HDRS} -- ${CPPFLAGS} ${CXXFLAGS}
plot: plot1.png
plot1.png: test.txt plot1.gnuplot
gnuplot plot1.gnuplot
.PHONY: all lint

View File

@ -1,10 +0,0 @@
#include "pendulum.hpp"
void adaptiveRK4Step(State x, double dt, double accuracy, State f);
void adaptiveRK4Step(State x, double dt, double accuracy, State f)
{
}

View File

@ -1,246 +0,0 @@
#include <cmath>
#include <cstdlib>
#include <iostream>
#include <sstream>
#include <string>
using namespace std;
#ifdef __APPLE__ // Mac OS X uses different header
# include <GLUT/glut.h>
#else // Unix and Windows
# include <GL/gl.h>
# include <GL/glu.h>
# include <GL/glut.h>
#endif
#include "vector.hpp" // vectors with components of type double
#include "odeint.hpp" // ODE integration routines, Runge-Kutta ...
using namespace cpl;
const double pi = 4 * atan(1.0);
const double g = 9.8; // acceleration of gravity
double L = 1.0; // length of pendulum
double q = 0.5; // damping coefficient
double Omega_D = 2.0/3.0; // frequency of driving force
double F_D = 0.9; // amplitude of driving force
bool nonlinear; // linear if false
Vector f(const Vector& x) { // extended derivative vector
double t = x[0];
double theta = x[1];
double omega = x[2];
Vector f(3); // Vector with 3 components
f[0] = 1;
f[1] = omega;
if (nonlinear)
f[2] = - (g/L) * sin(theta) - q * omega + F_D * sin(Omega_D * t);
else
f[2] = - (g/L) * theta - q * omega + F_D * sin(Omega_D * t);
return f;
}
Vector x(3);
void getInput() {
cout << " Nonlinear damped driven pendulum\n"
<< " --------------------------------\n"
<< " Enter linear or nonlinear: ";
string response;
cin >> response;
nonlinear = (response[0] == 'n');
cout<< " Length of pendulum L: ";
cin >> L;
cout<< " Enter damping coefficient q: ";
cin >> q;
cout << " Enter driving frequency Omega_D: ";
cin >> Omega_D;
cout << " Enter driving amplitude F_D: ";
cin >> F_D;
cout << " Enter theta(0) and omega(0): ";
double theta, omega;
cin >> theta >> omega;
x[0] = 0;
x[1] = theta;
x[2] = omega;
}
double dt = 0.05;
double accuracy = 1e-6;
void step() {
adaptiveRK4Step(x, dt, accuracy, f);
}
double frames_per_second = 30; // for animation in real time
void animation_step() {
double start = x[0];
clock_t start_time = clock();
step();
double tau = 1.0 / frames_per_second;
while (x[0] - start < tau)
step();
while ((double(clock()) - start_time)/CLOCKS_PER_SEC < tau)
;
glutPostRedisplay();
}
void drawText(const string& str, double x, double y) {
glRasterPos2d(x, y);
int len = str.find('\0');
for (int i = 0; i < len; i++)
glutBitmapCharacter(GLUT_BITMAP_HELVETICA_12, str[i]);
}
void drawVariables() {
// write t, theta, omega values
glColor3ub(0, 0, 255);
ostringstream os;
os << "t = " << x[0] << ends;
drawText(os.str(), -1.1, -1.1);
os.seekp(0);
os << "theta = " << x[1] << ends;
drawText(os.str(), -1.1, 1.1);
os.seekp(0);
os << "omega = " << x[2] << ends;
drawText(os.str(), 0.3, 1.1);
}
void displayPendulum() {
glClear(GL_COLOR_BUFFER_BIT);
// draw the pendulum rod
double xp = sin(x[1]);
double yp = -cos(x[1]);
glColor3ub(0, 255, 0);
glBegin(GL_LINES);
glVertex2d(0, 0);
glVertex2d(xp, yp);
glEnd();
// draw the pendulum bob
glPushMatrix();
glTranslated(sin(x[1]), -cos(x[1]), 0);
glColor3ub(255, 0, 0);
const double r = 0.1;
glPolygonMode(GL_FRONT, GL_FILL);
glBegin(GL_POLYGON);
for (int i = 0; i < 12; i++) {
double theta = 2 * pi * i / 12.0;
glVertex2d(r * cos(theta), r * sin(theta));
}
glEnd();
glPopMatrix();
// write t, theta, and omega
drawVariables();
// we are using double buffering - write buffer to screen
glutSwapBuffers();
}
bool clearPhasePlot;
void displayPhasePlot() {
static double thetaOld, omegaOld;
if (clearPhasePlot) {
glClear(GL_COLOR_BUFFER_BIT);
clearPhasePlot = false;
thetaOld = 2 * pi, omegaOld = 2 * pi;
// draw axes
glColor3ub(0, 255, 0);
glBegin(GL_LINES);
glVertex2d(0, -1); glVertex2d(0, 1);
glVertex2d(-1, 0); glVertex2d(1, 0);
glEnd();
}
// draw the phase point
double theta = x[1];
while (theta >= pi) theta -= 2 * pi;
while (theta < -pi) theta += 2 * pi;
double omega = x[2];
glColor3ub(255, 0, 0);
if (abs(theta - thetaOld) < pi && abs(omega) < pi
&& abs(omega - omegaOld) < pi) {
glBegin(GL_LINES);
glVertex2d(thetaOld / pi, omegaOld / pi);
glVertex2d(theta / pi, omega / pi);
glEnd();
}
thetaOld = theta, omegaOld = omega;
glPopMatrix();
// write t, theta, and omega
glColor3ub(255, 255, 255);
glRectd(-1.2, 1, 1.2, 1.2);
glRectd(-1.2, -1, 1.2, -1.2);
drawVariables();
glutSwapBuffers();
}
void reshape(int w, int h) {
glViewport(0, 0, w, h);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
double aspectRatio = w/double(h);
double d = 1.2;
if (aspectRatio > 1)
glOrtho(-d*aspectRatio, d*aspectRatio, -d, d, -1.0, 1.0);
else
glOrtho(-d, d, -d/aspectRatio, d/aspectRatio, -1.0, 1.0);
glMatrixMode(GL_MODELVIEW);
}
bool running; // is the animation running?
bool phasePlot; // switch to phase plot if true
void mouse(int button, int state, int x, int y) {
switch (button) {
case GLUT_LEFT_BUTTON:
if (state == GLUT_DOWN) {
if (running) {
glutIdleFunc(NULL);
running = false;
} else {
glutIdleFunc(animation_step);
running = true;
}
}
break;
case GLUT_RIGHT_BUTTON:
if (state == GLUT_DOWN) {
if (phasePlot) {
glutDisplayFunc(displayPendulum);
phasePlot = false;
} else {
glutDisplayFunc(displayPhasePlot);
clearPhasePlot = phasePlot = true;
}
glutPostRedisplay();
}
break;
}
}
int main(int argc, char *argv[]) {
glutInit(&argc, argv);
getInput();
glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);
glutInitWindowSize(600, 600);
glutInitWindowPosition(100, 100);
glutCreateWindow(argv[0]);
glClearColor(1.0, 1.0, 1.0, 0.0);
glShadeModel(GL_FLAT);
glutDisplayFunc(displayPendulum);
glutReshapeFunc(reshape);
glutMouseFunc(mouse);
glutIdleFunc(NULL);
glutMainLoop();
}

View File

@ -1,76 +0,0 @@
#include <cmath>
#include <cstdlib>
#include <iostream>
#include <popt.h>
using namespace std;
#include "pendulum.hpp"
#include "odeint.hpp" // ODE integration routines, Runge-Kutta ...
using namespace cpl;
const double pi = 4 * atan(1.0);
const double g = 9.8; // acceleration of gravity
double L = 1.0; // length of pendulum
double q = 0.5; // damping coefficient
double Omega_D = 2.0/3.0; // frequency of driving force
double F_D = 0.9; // amplitude of driving force
bool nonlinear; // linear if false
State f(const State& s) { // extended derivative vector
double t = s.t;
double theta = s.theta;
double omega = s.omega;
State f; // Vector with 3 components
f.t = 1;
f.theta = omega;
f.omega = (- (g/L) * (nonlinear ? sin(theta) : theta ) -
q * omega + F_D * sin(Omega_D * t));
return f;
}
int main() {
const struct poptOption options_table[] = {
{"length", 'l', POPT_ARG_DOUBLE, &L, 0, "Length of pendulum", NULL},
{"damping", 'q', POPT_ARG_DOUBLE, &q, 0, "Damping coefficient", NULL},
{"omegaD", 'o', POPT_ARG_DOUBLE, &Omega_D, 0, "Value of Omega_D", NULL},
{"F_D", 'f', POPT_ARG_DOUBLE, &F_D, 0, "Value of F_D", NULL},
{"theta_0", 't', POPT_ARG_DOUBLE, &theta, 0, "Value of theta(0)", NULL},
{"omega_0", '0', POPT_ARG_DOUBLE, &omega, 0, "Value of omega(0)", NULL},
{"tmax", 'm', POPT_ARG_DOUBLE, &tMax, 0, "Integration time", NULL},
POPT_AUTOHELP
{NULL, 0, 0, NULL, 0, NULL, NULL}};
poptContext opt_con = poptGetContext("sho", argc, argv, options_table, 0);
poptSetOtherOptionHelp(opt_con, "[OPTIONS]*");
poptGetNextOpt(opt_con);
double dt = 0.05;
double accuracy = 1e-6;
double t = 0;
State x;
x.t = t;
x.theta = theta;
x.omega = omega;
while (t < tMax) {
adaptiveRK4Step(x, dt, accuracy, f);
t = x.t, theta = x.theta, omega = x;
if (nonlinear) {
while (theta >= pi) theta -= 2 * pi;
while (theta < -pi) theta += 2 * pi;
}
cout << t << '\t' << theta << '\t' << omega << '\n';
}
cerr << "Finished." << endl;
}

View File

@ -1,9 +0,0 @@
#ifndef PENDULUM_H
#define PENDULUM_H
struct State {
double t;
double theta;
double omega;
};
#endif

View File

@ -1,19 +0,0 @@
set encoding utf8
set terminal png size 800,600 enhanced font 'IBM Plex Sans Text,14'
set output "plot1.png"
set xlabel 'idő [ms]'
set ylabel 'kitérés [m]'
set style line 100 lt 1 lc rgb "grey" lw 0.5
set grid ls 100
set style line 103 lw 4 lt rgb "#b8ee30"
set style line 104 pt 1 ps 1 lw 1 lt rgb "#8b1c62"
set title 'Harmonikus oszcillátor'
plot 'test.txt' title 'kitérés' with points ls 104