Commits

kazeevn committed f60fc35

Finalized report.

Comments (0)

Files changed (8)

experiments/particle_size_peak_width/Makefile

-generate: FORCE clean
+generate: FORCE
 	python generator.py
 
 %.cfg: FORCE

experiments/particle_size_peak_width/example.cfg

 refractiveIndexStep          = 0.0002;
 outputFileName               = "results/ns.dat";
 
-numSteps		     = 30;
+numSteps		     = 40;

experiments/particle_size_peak_width/gauss.res

- results/r0.000316228.res 5259.2412516 25.698140693 12.4097947891
+ results/r0.0001.res 5259.93031459 17.1180053469 2.80428460993
+ results/r0.000316228.res 5259.25519823 25.6878108702 12.4047298026
  results/r0.000749894.res 5258.35015407 32.9797522383 29.9912317548
  results/r0.00177828.res 5257.34728265 47.2395585293 63.0021145047
+ results/r0.001.res 5257.82247952 36.7093257329 38.9134275676
  results/r0.00421697.res 5257.68720751 71.1049963587 115.881014449
  results/r0.01.res 5259.07201353 101.42841392 187.083063942
+ results/r4.53999e-05.res 5260.46551999 16.9844687403 1.14433485551

experiments/particle_size_peak_width/generator.py

 cfg = libconfig.ConfigLib()
 cfg.read('example.cfg')
 initial_r_log = np.log10(float(cfg['averageSphereRadius']))
-for radius in np.logspace(initial_r_log - 2, initial_r_log, num=5):
+for radius in np.logspace(np.log10(np.exp(-11)), np.log10(np.exp(-10)), num=5):
      cfg['averageSphereRadius'] = radius
      cfg['outputFileName'] = "results/r%g.dat" % radius
      cfg.write("len%g.cfg" % radius)

experiments/particle_size_peak_width/plot_column.py

 # Plots the third column.
 
 # TODO(kazeevn) acquaint it with exponential notation
-get_number = re.compile("\d+\.?\d*")
+get_number = re.compile("\d+(\.\d+)?(e-?\d+)?")
 for filename in sys.argv[1:2]:
     infile = open(filename)
     x = []
     for line in infile:
         # I'm terribly sorry if your filenames contain spaces
         the_split = line.split()
-        y.append(np.log(float(the_split[2])))
-        x.append(np.log(float(get_number.search(the_split[0]).group(0))))
+        y.append(np.log10(float(the_split[2])))
+	print float(get_number.search(the_split[0]).group(0))
+	print the_split[0]
+        x.append(np.log10(float(get_number.search(the_split[0]).group(0))))
     data = Serie(np.array(x), np.array(y))
 #ax.plot(x, y, '.', ms=10, label=filename)
 fig = data.plot(withfit=True)
-pl.xlabel("Average particle radius, ln(m)")
-pl.ylabel("Peak width, $\ln(\AA)$")
+pl.xlabel("Average particle radius, lg(m)")
+pl.ylabel("Peak width, $\lg(\AA)$")
 pl.legend(["$I=C\cdot l^{%.2f}$"%data.lfit_tg()[0], ])
 pl.show()
 
+data.save_plot("pics/plot.pdf")

report_spring_2013/report_spring_2013.pdf

Binary file modified.

report_spring_2013/report_spring_2013.tex

 
 \hypersetup{
   pdftitle={Mathematical modeling of Christiansen filter},
-  pdfauthor={Nikita Kazeev, Alexander Konovalov}
+  pdfauthor={Nikita Kazeev, Alexander Konovalov, Artem Baron}
 }
 
+\newcommand{\graph}[3]{ % Takes caption, label and image file
+\begin{landscape}
+\begin{figure}
+ \caption{#1}
+ \label{gr_#2}
+ \centering
+ \includegraphics{pics/#3}
+\end{figure}
+\end{landscape}
+}
+
+
 \title{Численное моделирование фильтра Христьянсона}
-\author{Никита Казеев \and Александр Коновалов \and МФТИ}
+\author{Никита Казеев \and Александр Коновалов \and Артём Баранов}
 \date{16 мая 2013}
 
+
 \begin{document}
 \maketitle
 \section{Введение}
-Фильтр Христьянсона -- это ... Работает так ...
+Если заполнить кювету жидкостью и твёрдыми частицами, такими, что
+показатели преломления совпадают для определённой длины
+волны, то получится спектральный фильтр -- для этой длины волны кювета
+будет изотропной средой, для других -- нет, и эти составляющие будут
+рассеиваться. \cite{clarke}
 \section{Описание метода}
-\subsection{Физика}
-Шарики, геометрическая оптика. Случайная генерация столкновений.
+Нашей задачей была разработка эффективного метода численного
+моделирования фильтра Христьянсона.
+\subsection{Физическая модель}
+Мы ограничились частицами с размером большим $100 \mu m$, что много
+больше длины волны $\lambda \approx 0.5 \mu m$ и использовали
+приближение геометрической оптики. В рамках этого предположения
+трассировались входящие в кювету лучи.
+\par Для упрощения расчёта, частицы считались шарами и вторичные
+отражения внутри частиц отбрасывались. \footnote{Во-первых, их мало,
+  во-вторых, если луч вышел в сторону, то скорее всего он всё равно
+  рассеется}.
+
+\par Вместо заданных положений шаров и расчёта пересечений, мы
+воспользовались неупорядоченностью среды и генерировали столкновения
+случайным образом, считая распределение расстояния до следующего
+столкновения Пуасоновым. Покажем это. Объём цилиндрической кюветы
+\begin{equation}
+  V_t = \pi  R_t^2 L_t,
+\end{equation}
+где $R_t$ и $L_t$ -- соответственно её радиус и длина. Всего сфер будет
+\begin{equation}
+  N_s = \frac{fV_t}{\frac{4}{3}\pi R_s^3},
+\end{equation}
+где $f$ -- плотность упаковки, $R_s$ -- радиус сферы. Тогда
+концентрация центров сфер
+\begin{equation}
+  n_s =\frac{N_s}{V_t}.
+\end{equation}
+
+Рассмотрим луч. Окружим его воображаемым цилиндром радиуса
+$R_s$. Предположим, что центры сфер равномерно распределены по
+объёму. \footnote{Вообще говоря, это очевидно неверно -- центры не
+  могут быть ближе, чем $R$ друг к другу.} Тогда столкновение луча со следующей сферой будет
+Пуассоновым процессом и среднее расстояние до него будет
+\begin{equation}
+  z_\text{mean} = \frac{1}{\pi R_s^2 n_s}.
+\end{equation}
+
+\par Для отдельного фотона важно лишь то, попал ли он на
+детектор. Таким образом, пропускаемость суть есть отношение числа
+успешных испытаний к общему:
+\begin{equation}
+  I = \frac{m}{N},
+\end{equation}
+где $N$ -- общее число лучей, $m$ -- число попавших на детектор.
+Необходимое число лучей может быть оценено по
+теореме Бернулли
+\begin{equation}
+  P\left(\left|\frac{m}{N}-p\right|\ge \epsilon\right)\le\frac{p(1-p)}{N\epsilon^2},
+\end{equation}
+где, $p$ -- вероятность
+прохождения, $\epsilon$ -- любое, по смыслу погрешность определения
+вероятности. В частности, при точности $0.01$ и допустимой вероятности
+отклонения $0.01$, в худшем случае, $p=0.5$, $N=2.5\cdot10^5$. Так как
+$p$ оценивается как $\frac{m}{N}$ -- стохастическая величина --
+количество испытаний было ограничено снизу 100.
 \subsection{Реализация}
-MPI/C++, сейчас параллелен по значениям n -- 30 нод, плюс по
- параметрам установки -- 30x5 нод. Потенциально -- неограниченно, каждый фотон.
+Написана программа на C++ с использованием MPI. Параллельность по
+значениям $n \approx 30$ узлов.  Потенциально -- по каждому фотону
+$\approx 30000$ узлов. Расчёт каждой из зависимостей ниже занял около
+5 минут на персональном компьютере.
+
 \section{Результаты}
-Графики - Гаусс, зависимости ширины/высоты пика. Сравнение с эксперементом.
+Для установки описанной в экспериментальной работе \cite{mc_alister}
+была рассчитана зависимость ширины пиков от различных
+параметров. Получена формула связывающая ширину полосы пропускания с
+параметрами установки:
+\begin{equation}
+  \Delta \lambda = C \frac{R_s^{0.35}}{L_t^{0.40}f^{0.33}}(1+2.6\cdot10^3\sigma_n),
+\end{equation}
+где $C$ -- константа, $\sigma_n$ -- стандартное отклонение показателя преломления. В
+работе \cite{prost} приведется следующая формула:
+\begin{equation}
+  \frac{\Delta\lambda}{\lambda} = C \sqrt{\frac{1}{f L_t R_s}}.
+\end{equation}
+Основное отличие -- отрицательная зависимость от радиуса частиц, что
+странно -- чем больше частицы, тем меньше актов преломления
+произойдёт, и, казалось бы, меньше света рассеется. Что и наблюдалось
+в экспериментальной работе \cite{bala}.
+\section{Дальнейшие перспективы}
+Больше всего времени у нас ушло на реализацию. Она прекрасна. В первую
+очередь, все \footnote{Кроме приближения геометрической оптики}
+недостатки модели можно устранить с небольшими
+трудозатратами. Следующим разумным шагом будет расчёт какого-то
+другого \footnote{У \cite{mc_alister} графики $n(\lambda)$ и спектры
+  пропускания приведены для разных условий -- положение максимума не
+  совпадает с пересечением  дисперсионных кривых.} эксперимента, сравнение результата и
+  соответствующая доработка моделей.
 
-\section{Дальнейшие перспективы}
-Есть метод, есть реализация, есть Issues, куда расти. Кто угодно может
-взять и продолжить.
+  \par Наша реализация доступна в виде исходных кодов по адресу
+  \url{https://bitbucket.org/alexknvl/poissonraytracer} под
+  BSD-подобной лицензией.
 
-\end{document}
+\begin{thebibliography}{9}
+\bibitem{clarke} \textit{A Theory for the Christiansen Filter},
+  R.~H.~Clarke, Applied Optics, May 1968, Vol. 7, No 5
+\bibitem{mc_alister} \textit{The Christiansen light filter: its
+    advantages and limitations}, E.~D.~McAlister, Smithsonian
+  miscellaneous collections, vol.~93, num.~7
+\bibitem{prost} \textit{The influence of the Christiansen effect on
+    I. R. spectra of powders}, R.~Prost, Clays and Clay MInerals,
+  1973, Vol. 21, pp. 363--368
+\bibitem{bala} \textit{New Christiansen filters}, K.~Balasubramanian,
+  M.~R.~Jacobson, and H.~A.~Macleod, Applied optics, Vol. 31, No 10., 1 April 1992
+
+\end{thebibliography}
+
+
+\graph{Зависимость ширины полосы пропускания от плотности упаковки
+  частиц}{graph_form_factor}{form_factor_peak_width.pdf}
+\graph{Зависимость ширины полосы пропускания от длины
+  фильтра}{graph_length_peak_width}{length_peak_width.pdf}
+\graph{Зависимость ширины полосы пропускания от стандартного
+  отклонения показателя преломления}{gr_n_deviation}{n_deviation_peak_width.pdf}
+\graph{Зависимость ширины полосы пропускания от размера частиц для
+  плотности упаковки $f=0.7$}{gr_part_size_dense}{particle_size_peak_width_dense.pdf}
+\graph{Зависимость ширины полосы пропускания от размера частиц для
+  плотности упаковки $f=0.3$}{gr_part_size_sparse}{particle_size_peak_width_sparse.pdf}
+\end{document}

utils/libconfig.py

         outfile = open(filename, 'w')
         for field, data in self.iteritems():
             try:
-                formatted_data = '%g' % float(data)
+                formatted_data = '%g' % (float(data))
             except ValueError:
                 formatted_data = '"%s"' % data
             outfile.write(field + " = " + formatted_data +";\n")