From 12e848205576cfd54d5c2ab0f3affca57c9f8361 Mon Sep 17 00:00:00 2001 From: Zsombor Barna Date: Wed, 28 Feb 2024 09:39:19 +0100 Subject: [PATCH] add old lab3 and basic szamszim --- modfizlab/old-lab3.org | 34 + szamszim/cpl/Makefile | 36 ++ szamszim/cpl/cpl.tar.bz2 | Bin 0 -> 13600 bytes szamszim/cpl/datafit.cpp | 178 ++++++ szamszim/cpl/datafit.hpp | 23 + szamszim/cpl/fft.cpp | 214 +++++++ szamszim/cpl/fft.hpp | 75 +++ szamszim/cpl/interpolate.cpp | 59 ++ szamszim/cpl/interpolate.hpp | 19 + szamszim/cpl/matrix.cpp | 895 ++++++++++++++++++++++++++ szamszim/cpl/matrix.hpp | 86 +++ szamszim/cpl/odeint.cpp | 136 ++++ szamszim/cpl/odeint.hpp | 44 ++ szamszim/cpl/optimize.cpp | 134 ++++ szamszim/cpl/optimize.hpp | 47 ++ szamszim/cpl/random.cpp | 73 +++ szamszim/cpl/random.hpp | 23 + szamszim/cpl/vector.cpp | 167 +++++ szamszim/cpl/vector.hpp | 124 ++++ szamszim/f1/Makefile | 30 + szamszim/f1/chez.c | 46 ++ szamszim/f1/euler-divergence.txt | 1001 ++++++++++++++++++++++++++++++ szamszim/f1/long.txt | 1001 ++++++++++++++++++++++++++++++ szamszim/f1/plot1.gnuplot | 19 + szamszim/f1/sho.cpp | 77 +++ szamszim/f1/sho.o | Bin 0 -> 23088 bytes szamszim/f1/short2.txt | 201 ++++++ szamszim/f2/Makefile | 30 + szamszim/f2/odeint.cpp | 10 + szamszim/f2/pendulum-gl.cpp | 246 ++++++++ szamszim/f2/pendulum.cpp | 76 +++ szamszim/f2/pendulum.hpp | 9 + szamszim/f2/plot.gnuplot | 19 + 33 files changed, 5132 insertions(+) create mode 100755 modfizlab/old-lab3.org create mode 100644 szamszim/cpl/Makefile create mode 100644 szamszim/cpl/cpl.tar.bz2 create mode 100644 szamszim/cpl/datafit.cpp create mode 100644 szamszim/cpl/datafit.hpp create mode 100644 szamszim/cpl/fft.cpp create mode 100644 szamszim/cpl/fft.hpp create mode 100644 szamszim/cpl/interpolate.cpp create mode 100644 szamszim/cpl/interpolate.hpp create mode 100644 szamszim/cpl/matrix.cpp create mode 100644 szamszim/cpl/matrix.hpp create mode 100644 szamszim/cpl/odeint.cpp create mode 100644 szamszim/cpl/odeint.hpp create mode 100644 szamszim/cpl/optimize.cpp create mode 100644 szamszim/cpl/optimize.hpp create mode 100644 szamszim/cpl/random.cpp create mode 100644 szamszim/cpl/random.hpp create mode 100644 szamszim/cpl/vector.cpp create mode 100644 szamszim/cpl/vector.hpp create mode 100644 szamszim/f1/Makefile create mode 100644 szamszim/f1/chez.c create mode 100644 szamszim/f1/euler-divergence.txt create mode 100644 szamszim/f1/long.txt create mode 100755 szamszim/f1/plot1.gnuplot create mode 100644 szamszim/f1/sho.cpp create mode 100644 szamszim/f1/sho.o create mode 100644 szamszim/f1/short2.txt create mode 100644 szamszim/f2/Makefile create mode 100644 szamszim/f2/odeint.cpp create mode 100644 szamszim/f2/pendulum-gl.cpp create mode 100644 szamszim/f2/pendulum.cpp create mode 100644 szamszim/f2/pendulum.hpp create mode 100755 szamszim/f2/plot.gnuplot diff --git a/modfizlab/old-lab3.org b/modfizlab/old-lab3.org new file mode 100755 index 0000000..558dad8 --- /dev/null +++ b/modfizlab/old-lab3.org @@ -0,0 +1,34 @@ +#+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 diff --git a/szamszim/cpl/Makefile b/szamszim/cpl/Makefile new file mode 100644 index 0000000..b2775a0 --- /dev/null +++ b/szamszim/cpl/Makefile @@ -0,0 +1,36 @@ +################################################## +# 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 $@ + diff --git a/szamszim/cpl/cpl.tar.bz2 b/szamszim/cpl/cpl.tar.bz2 new file mode 100644 index 0000000000000000000000000000000000000000..734df5cbb258364cca493bd88bcd1346c60521b5 GIT binary patch literal 13600 zcmV+*HQ&lYT4*^jL0KkKSz2l8HULH0|G)qL|GXO#lSbKp+BW3VIa%8Z=E$(wQ{V z(li4=GXjIT#R1JX2kK+&{-0iXa9 z5RxW}o|80&%6gxwdY`HuLnwNiP-yiuc}J?Dl)AI?Q{h?(j(xfU*Ig8GS=gd#=k^%epAF z9q!*~2e}EfA&D>&Dn(Ib6sL)@t1TJtQ0(BO4p+;^ctd<3) zTY%P9Lbeeo=tEk%L+fYNYNoUd5j^=pKK^mA@K_se*4iauQ;^dV*py?6uToK2{tl-k z7n31Ghag~##)}XNFWxF4q0Nd>#0`?GC?ZzGnPMiG34}qhMkZ7`t-(=E1u#&u)?%gz zX4O?i7Ezdn9>p{y1VKU4B^Z>|7{mK>)WES~0-F#(N;|i!%X#ohz8HL2mVHA9fpbL@-yCkJUVKUK7Y!bv0;9xQ(PAOo;Ktxp$ zf{aI0sBSxme1=;x`t5Wy8hN*1GSo2x8VF^J4<%rR?g3Ns{boW#Uh zh-ox279gTwqLxQaps3VPQDCuGcNrDSQPFNv*3DgHH5g2pO3ZDvNU&12mgN#vu|nHm z!x2GXDdRSlY;77gf;Cbiuo-NXicOpufXQmAwP~Nm=&Ux|3K0e3T{ymNId{!~GS?JU zJpf`rVg;56AdD6sZx&p2=d={C7@-6RECWePq8Ka#lWk9u*iLXf$-z8ofGQ+d5!EV1 zQX(p=BOZ(d5JGiHBo7E19;bZ)>4nqADCZchItbZCRLvHG>Ou|qfcxrHz04~%%M4y< zP||ceF<@`wL`nxk(r?y%b+0#kxwXW!Su`+ge~~c9kRW?`A`6qOTueBh&G#Z6B$CUsh0#U(BH+R=w(%rLd8lnE<$!MFM+mc185Yc?9KAULZD-! z-_f(6>C=K2KsaC+P*wAxKo09mEZg0T1kmG9W7=0Iq2uR3UUc z)q`e=+>3Nr_$dMXh7zjKl1QSp5meuMd9c7UiY2i}l;{K{{HC(#l}dt4yZSRAp>wM8 ztQ87kQ6UlpP>00W1aA(J36w#m5f}|(VY|C&iT1hY*~jc2@^I=mkKH_Hc*0mJ zh%0@LmZ!gu-0Y+DcuM~*NS|?3j0vb_NZwlhOG0n$xO7=u9*RQvh6kb_K7Jiw+UBNa zYaL2l*4rx(u@;D^Zj(YC2!sSlB#9A>kac8~T~eH?uf3sjS`-D*3nId(BzGQ6!fp0A zdO4L9e4hlM^!r2?lj?`s_tWQ%-&fJ>im&n0R|OySmY;UZWON;Qd{dN|p;Hsl@MF>k ztOWydUU+cz(r!^My{Dd@;TfnwD`z5c_)o9JHM+HGeEiVRFQjs;nVX@`bGk8m-u`F3 z?|hIa*1x6l;6Ny{4`G7=?(NYV0IF^dDFYL1g^ zNC?CMoah?W?P4-up6+}&7kg4OYKgC>#(nRO`C+w44;k(b-F;am zUAGpBnN}ho{~g7MTsA*h^%L&s+l;bZ%SlMZUU@p@ORkfTcv;hbpIJQ5SowGWiB5sFsXQl;bz3zT#PwUd^vcxw+RR32E{>tK7twg7S`}2)TR+|)CHd1SDDbd7{ z=@{zd{A{-USVZ7lGGdB_3uxP-4-SrUNDdM!(HnWL)j3t=UJrZoJG!HXo+;zi7i*om z(Qw-UV&Zy&q^l)Pu^B38-J6n$hYFM&v2k9x^6lef$LDRFm3eaH+)Jn9wT?fnv~nhf zsb{w(gQE%CV(P84qo9h}XP0xr;YaYiS5CF@UGrd&d`uLWz3`sK7{8Ev6Qo#IKh5rn zf#_1|!{}&W*Rs$$K}P;U-(vN1EPT@)UKeuX5fLPkNhFdn%SsB{3?^n|HR51%_~*M1hIR@xPerb^FL0UP zco~r*Dvj0=>{W3OVQ|Q*vcl-@-2);G8Wj~70$!MqEoxOySt@|*x?%Q=iI{O0VavWML?+m zM+$I2NQ`}FBiq{GW~On~1$JCpH%{*jYr6N1&LinWC`=CYm!4-317PdEFz zJOFGMBJCu3FnIYxWlL8sIL;p1h7EvB4kq+4CLmFeV2Haft&y{B2AClj8-*-)i64CH z>+;Awf<-5LhYGF1A0y@E+3}D)s(@k|LQ37Ul@g7_t&4hE5#5Rn2D%#JG*@#gi<5*w z(Rp){5CxcE&z#*kAgcpxsfhr>(rgrJlrD9sQI~~By5T9v<5qOfmyw>wVzGPQ#mp$? z<-v2_`q0n2x`Xh>4Z76=ObpP>%*dAtTnI^_rm`j+12#Og955EYhTPQN_i~R#<)lvt zXBU$V#Mp8T`5r>$wuyrssV35KLV#fy6kc@$TZ*#!w2f0`x=pOYXQ{l*(|jioLMVhX zz=g0L~9vV{y?h5Y?J9Ge~ijD?;U)X}X3vW!gw6C2UCHRc*& zvVgpYwZP08K(K%kFowq+)C2Y$Cyfmxtw*nA7$TW0j?*k0^DO*9pk%GB=}QpAh7w`{ z(3GRZu(F6iS*SE{sapp=2TrR=V`Vr(V^9_3VGNmS3;h+v(jp3tNon8J&sixXj(Y)x z6xuQ}Fu=gW3@}h;$Q^62EVorUU|k@@$mPO`&9fqe7d1CPVzxD#j_z;{optaVj^$1g z1~Sm#Zw00FBT`8$nodIuaws57Er)6jQAI+T4oey~U7I3zgTf-(4GX~9Yzit6Gm1%G zL49`UyL&dsX_HiXh>|Qr6g41x=R~Kd7c0|`Xvof^U~{fn_E_m1fJ368;#w$0MN-OY zEKw_7grv-?q7hPoayc|>oGL`B>(-5}ym1n9<@Yvaozq>ARut?+f*3aj%=+YIf^iZJ z$1xsuilQ1$_VJ5((U6vLSvAqo#-^v5rr_M z78cPh(~?($@Ko3{2~5+2mJSHi@nQCloQT;b85dE1>`Y--?8R42;t-9#B0N}jpon-KAbCi=TntG=yng0A)eF?Tl&$>! zw-W+*I*~|-2^a_?1^Zxs-@-}(KtdoSfr6g9Rj-{2p3L>T2c!a|VM-zECIXWND4{6E zD0)uQK?2)twSO1O3sMEEFbF^}i6RdONR+_lXd$8`vI{3+Zs(L*O{c%F7|(J&5I<8X3bt+!ETj{z}RZ49fx5a zDxJPW8c0oA;=_>KrzJxH9ls4&(O|WEXey=cGq0;*{KWM(!PL^+QkhG++GW@slcK07 zhE$RQh{`%f`=pnw&iOC!-P9pSKD-dYr`>Qkxdk9N#6_^^&kn7XOo?f34MPdK30l~p z=nV!PyCiL;yVT2zMz>A5k;ahuDG1t|Uq=$Cdb~T*6+ons?_~-G)&d&bu<@u4E zHj)nNiB>>zoK{YMez z8Z4@mr?%)3k$+&e-TM;`h0%5j^@K0<`0+|?<{XZ2q-RY|u3?ksQll70&^R@2qtLX& zox^XHJ_b42fXHj^jMdKt+_akw1AY!!KDSRQTH~`-+lJDHplj14j#sn@=VJ2MrqIkd zg!^C1PrvZTtSGY^H_HD#&?t=g;me+bqf{-UzIS_JT1{l;!#|qM@tNG=t<}Gb4 ze@1g!5Bet8A!qGF)X4YhLr-FSt(PYcdCR(PBqqKSD&3X`lXIYJlUwt(nNtDref~cF zt_rKcvy*y@+wy2CwmpZkY}NYg_$RwjJjv_F)%)8$X^eDyS$+O1@V$%g)r<{J=zU28 zk_V(jJ|I1&Vu~ubL?i)%B3ppTfV*~l9ZmgsH*+s%ztt$kB#!R@ccw_=1yvn<6jp*07c{T^mRe6Esaic@WSsPsC5ydkM-D7!)uPR6?u&U`@jqLN4|a8RH+yV-Ut`p{xR zA3&s}D9J$z#1&<=5VDD6wk0S9Eh`@;Vqs^x^|E`7O$w0S?=YdU!w6%4q4qg+lKuW` z@#mQ^#m077teUn65iRR?Ulv5mR@L07>6QDmgUw={)f8T@6oggDO|L& z!uU7AW+eiXlKsC;I!n^`c1~_u*Q1(#hvl5x>3#al4e)MCQLXw-Ory_jtm1b(XiUB)t-uVxH z!iq`C`&HLjL03u$4X4NkccGf%OdpdTM1J%MIu(bZ;UGjhM;7RAH+(GYXVe*%b_nUz zE_7|gA-xGAVF&K5wF+P5AmgVZW|9>SVIj&cH}Qo6n!cAeL{%YP;g^>~zFs3G9K#cH z$Z?XNdp02k{FwgYV(< z^xfZ_UEOPz&nvAE&~fwmocYstH`{YuwN-bpoEB72T{DE3tQIK2hHY{QF)Bjp>C{h; z-tj@|8N0*YpO4%Hdp%0}@JHDX zVDTCZ8(4b(j{Of0bcbH(x~~@Z`NGPoFSjYi%lTdYj9Lvj;&ZRjtd#mT#_g7MwpnV> zE$PB^pfawR?<*wxcZhsqh>$@%UL0PXo42^<_}lNy$gyjrmR#t_`yGPRqdVF&X$wver;evtqGATg(|4jaw&H4U|GiIX6Z1bRKf(AMBh#w*W z{b-?71Gd{$pft1uTtWi-h~j{>cElCqe*481d8udT7EI72f$Tq+?|LRlghhm8QbmX> zZ~>YKZ~H*Znz6qU&|jZ&h@*Y@apQdW&4tsQH5A8y zn%Y*#E=IkbnV3+Fw&iNWIC4mfnkdr-nDB6DG>PXf`O8{hj1~S~S97ShDb22;6&{Mt zObP7<2d`gC0nX*n_qy{=^-2`?OhN+LvChh$6V=JWP5dG%-DFV4=XAqEU`0zo<5S!DyulVYWlj#@RKRoN*r+oCjx~Z;5xMbRWoHsUB=>>U49`7^4ta&4b{5LpI zz4TqgY1?8@eDC~?W{lv)MIidz`K zj!>a^Z*eMFHN5_A1zOf(<@ac1A?HjtjR-hFhu0CURP>~Ia*XthawimXMgcaH(Yh7l z@O;Rm$xe-)BW>e_(4F{m6Ie7u-WCC_?%89lfq0n z6tTkwro-n*5x;~m5w3&h#LFy9nlKzR%TmQmt)1k~sJj7ZrhHZ{Ok3|c&hj>IbBbW2 z>r1O;n#t%dL;b57h4!YGdAsKkXKGF=X;)HLI!8mnfYr=bYnHiZSu^CEl@3QUUaSYYd*36& zADWpl?;oo?;@ZAtJn6oojvk>6+{VY5R8l($BG?3H#0C*_?vl4-6BZ;v5Vsh?8CJc4 z0tQz8HgA1>`z~?&S9xH*w|ULqPfb{&QUGHnK6{ml2Q!*{Z|Lip$w_a0)&2inYx}+@ z&lyYM-nmKQ;UI)_nM4#rQGo(ZqX^EW!?As84x0}q;(CCD!2e^6e8`a=VK%oqh@{A1 zPxR9V&m&0F8)W_`bKuS_+IZbI8MkG3b=u=eOUP>BgVfMETdhr+L!V?^C)6)O=BmGz z#+##>xhxc={4|UVUm6>OK%B{zmgyk^+Fb3`Whh)>A>;33a?e5!AkS-%IRO4`FG7f; z;v=-Ia+kL$@`-Nm64|5Mcc6VI`6t4L`*+83zX`;g)=7z=s|1Oqc12e45b1~#v(Xmp z=`Mq$>uX&x3nM&6;(EEl0++t@X-i&rM7r6fqR8`Jb=My2cU0w+sUQIz*%(D3z3^oL5(yqatseP;T&zn(yhSks9XN## zLW~7f0wkqYC0?Xdvi0LNUDkT_crS@F?|fZ8uz~rEn?25=Eb8ovqRBMKS-qN@iRgRS zVxzAQEr$#e+ZyKEY*NKBy4^~xTE?`a9YONyx8G<3o{|Bfq7g|XL6KGrgTtc7=79(Y z<{-Tw!5o>c7Ld6_X{5oWXwFPX+Bi}vDgi{4LC&FAx&OX2GLFC(5VJl8f+LQ3|UYW0K_b6%Qh(% z4HP;zC_&1iB)Nf-3m~x(O3E=&pAWWsU4HeMe?p}7wKXxumF=g*W!<1Hu-N7%y7#lB zvht5Q%NRQ-JyU3uAc)5S2XIB=@byZ28yc1CZ=26_(BC4H-9C?};_o$4(jpj)lq#)P z_3wm8`gxwX@>OkoYUARt=Q53U`@T1;s^N3GPgIYG&ChRzG*Z=9F*HSU4U$cJtYvL> zAxorl9lfNfdb)fduIXL^cFWIt=E22_gv8hINf`EX9_>c0ImQ`RBPS3i8(zIK3^QB$ zTRKBNYOAG`T@AM98Ai4Cm1`X9$B!%%T_T@qh7s6{CN*G&ZB}e>W{q~(_)4lA_dTEx z7{(}wFrADjQ)|C(nNQ1^2p=*11%i7FAw2{-hsgOy*1o`Xmu-JEF9v*uxY@=t)*0*E z7vgJbn9XJDFr>=x>a@%GgwpUu2jmSlC{BX=>MBNWI4y}*(O;Is$s zZwl_Be3ikdC_az&CJzfvDHvd@G7k*sBzM=69-qq%wb6dD&{Cd$y@LId@TpZe+4f1? zzR?HXUE;V^|3=6229vbjZy(A!4de+VTOWP3u-D)#QY7r}b`%zLhG$Zpl*k3*6_2W7t4XvUqw#jV2@-BT zfZV^JwCEt*9UwmgtnT0&(Ab=yoTvSJ-&5b`(WO;IR>0PKAS(-^VM8FIz!oTIDpFsZ zGdHVKb6}{dHT?rTd4<5}H%{ZgbcfxRZwJH-CDAPCfsN)?Yusg3XefKuJ;d|zB zQh-^RhvC2l<`w1Z1&9m43&UdEdnas+z&t{db%XbgT4c584O4?NPD==J`>NhAqJjP4=N+$jp|e9U=`ZU=H2s;@#8 zKGoW43>6mvTr9GJAc&^2;b=CSgvHh0aZo)o>&x!#N%eRBkJ$QKLqWQ#A}!=IGLUf` zf*l44No1M|+z<*D6-cl~0~nw^bY8UERtlvI4>@Rp?^Me4DD*mT_P$0p5rm>Ri-M*C z;-F%c7#|Pd?Pu3IA%2|E*g??*!Dgn)cfD7);kTQ&D?Z#dJy%ufZrzhJ(qzpwy=qp7 zz$j`mm8Nv|4NF*iQ_yT*53M%&a>-*?(5%_hCM47KscI|95sDepV)^9Z@}xeHJ)vg{ zDr*^8^or_eO(PGg724Q#H6jQz88vYbN|vMTm9nL>DO)#w!WMq({2S9MX)virmrlqK zS+T!%>#|2h>tm*ISIVPs+t|k~)#&>PsydsmBn1?ciAk@I*X7+_r zOGqea{i9R}Z|oEzC2NWIRm*^@@IRMes)*U=|qCqie0vvdj}jOtC9Dn(NBI7LI7 z!c@1$n3P_^GIaQeMifabxsTixDN~iP*cJoL*SF`LBR^ZzSMDg;Vii=jeao0JO5Tpk zg^I@NJhv&1#;%DPlzr>185yPqhE|%@#@(W^*(e~`npbRdJeq?jrmt2$vKE1$DquR{ ze4^*z%US`oW)rGv?Evm+d)rj4-+omx6$F1ZuhHzt7m9k@2|`h*>2oCv2Grzrr7V=RPgaP`E2sF6^y@E}@i z^Kse1_;qwqZ!Ns`9#+KeA6X0!QjqlCd|TeRqgCqfE%pbEi#E>2c=JrH zy%jD=Q59*8EMtoAqE7-UI%KaONKTegx=ZGU;m~}wepRX}ilVKx`%n(BJ2@P2-$_+r zb~@u5ZERx-hijt2q|8YzK+qJW1V=iA6I1{p7m>}^s}p&_uM?^1O<7Lz=!7-nqPd|) z_&MuurZ;h<>Ocz^LiAytmZ)pvgiPiMpI1JD`;JHpW|~%#;r_FR-nEiGfr~k> z&=x^+MWjwzLcq=zhYz#=-BzZ!$n>*JH&uc!KaDF*>LPt$N#SL?#d%SuF1r2J6 znSFFN7z(2ta?SFYzdk{!{Yk$~XU~lQ^)N?k%phev{oZ~oT=COiU_Pu$ z$aR9c5RWFG7qdazaH(aEOsb<0Fr@R@LG2YwBcMT`s}M?L>G-8!(luETmZI5!380d? z%T$>mn_V1gCJ>l`wa1?--!{C7Z96lFY}zOmWH!htlb2dMZtg$`B5M9`?O zzO}t9%`ZjaQ5=H>1^`5QWDCbd%s>IcI}q1XJ0b)~$vZW!Q0T1LlAv<-X@q1#JQsjs zd-tvfwg+Y2u^e2ikQ`8O$`rKPRXAozECU>LTb9yg!L+MRHaGmSDWw^pLZF_P!9qb| z!b|o0xeo%AD`)|P6H@UnG2bD)pf(yH+qO~B_Fe!iIzvc~EC%e&nEz1XAqk&gV=>Kc{>{Vj716?&_jsp+tR)-o;4-z z04jjad4!TwG4sSSKP!*e-dfIVILiK%%mTV^tF>D zi;F}LMAN2}Y^>#kQk|Wfl>hJd{H`GCb0J@lYS#*s zFyWb{;31G%S#L5&mF!YOZ?vonsYz0`aKnLl3F-XO>XkJ!S#0wJ!;+~9@!Pp9Gr0@a z**%olNpI=L>2Jp4+1~X}Llknp2<~_tLQ-_%_KxXyJG8B5r%rX+!-i{+wwFnpWH?g1 zA6t$jO1_)Ha6~w4YHG~`JsHpY?|`B)sNQmJLS;GQ%nI`mm6vQ8lS!hTHPyFW)qZal z(9o^NOgii-zm|E?5mSik5lJ)h=p&380K{o^1`QFW{Huz=*g<4`HQAII++-pP0bnjd zP|&T5$CXn)y4v4BBK^Q-IJ4}(TM)OlF|d*qtXZC06unjf@-Y$!D3B0DMS}p_!9F<` zP~6&x^VHwf15o6oAA*?%^aA`Z6|M9(wyBCPh>B>4p;1*2ImgNCb&2kAW&nkNMK|$! zJ}se;uTYVDNd{#_*W0#yxV;CRGR4;Fhr`UsrAmyYmD4I+W+RERYgBbJDj3z(`1jF~ zZxpPdP%_gF#?sY{g42wH^FYCePJ3WS+j-*%Fazcu(vjlq@a5@r8@gAQnD_Jsuug~) zEQBD8i`G1Pb!zTElI6WX!Du}b+-y2C6;(&tz@4cIcn~_!vU7e>2{t@1WcjMEZ<{i` zxqQ0UVaN{BpfF(XN^~7RlGH$Q`vIA8J}Mo6cAWm0!Ga$Odg#Oyj08mN?ZP4%iIPbW zz&=*7)Cg%(1raQ)f+Y6=50|hSA%O@&FBIbC0CDQhBxVw4{IjpZoc$)G!;SF)~ zhSYm_ZuL$>wvGg}kC1Q^amNAA!sc@bDWGqHKsH3h)#xV%`b0}Ps7-N@VG<({`vcr~ zarFe+9&MfE-6VQIMiZ0>8?pg|1W1Z7Cy?0TH#7&^oU^Ku2M?;7+3YB(tt=-JMk*&x zj*WG9^!f?i&Z7p}f?HW>WmM(IPROK~Vj+Ot5<(EHl~X&j0Wq8fR?rU)TTuK4^&5t= zwJkvrBT8MsE<0~ZCRw`@i8a;s7^!5&a$qbY+#85Wf=g6q*FwmG@R3HK#3F$P^sI2qUs_ZxK?@ATaRZIOod)x< z)bxW~>0^LL8(_i&ktU)_4bH?b4jzi;Zwh)*(QyofM~qz%LM(D%L&1S2;8!q>!MY%0 zMICqm-bWJ&Naz`a#0^wzk)&5K)W;=@MhgeX&Mg;E4HksPlxtrq5I9Ga@wg~C0D7F1 z9uwe;NU67rwhF3?HklYC6CgxJH1j!@Y~hq?W?vym3mM?^oRS_63=F8EWl4@W7hRZ| zz=}{Q(s`jU4k{o*72$!gGxOxM9BG^VItn|H4|E^_wk@b`7f24W3lA5!T*5u)UqeG5 zPkss!kdcL>vzQGiVNMOEC@Vrz0<;BGka4NjVTpm{u45Sj3{tZa+j_?2RQ8Z7yNqlT zn?Cpl8x^7?h_VZwkkd$mlb~?!_$f3A{=Rt*K-k{C0L;olt*KB$BXX&Ny1634mg$w1 zfmktvkT`N>VAyJiAf_{p=v0OauZwQkvrj>6L9(N&G~V+DnX)X(twjzjXgakMXr*O6 zbvcdOO)@OsS*>5GX1+?>v4oj3gdu{U4V~c4T)L8OQc75o%u(~Y(hb!0mjr>2tZ*5@YBLZ(M#OkYc3Ev!vM^U zvOozjSl<%XDfB2Nt9(3qkeaSxd1A~D5b;unQyC7BPowyY;FtHK6FSR9+CxyL9BU9 z_X6Gkz+Ob+wuO*nS78`L&dkGVXH?AEHlTyRdxSN5z}*SNhoyj|778kn1`3b9B@rWp z`z_|T9gqiTA~03Fv8G7CY$Ajp02GTBBuGFAanhspxq|BrAn}F@K!g&+*c!>j-}T`& zpSa{6G6y2>4WM$V4O`GbK0VvM+y3U%G(7YyG712q4H(j{Dv%I9%RoFuknP%IA5ORN zTWoinkQkf<${q6j$lzvViUv-89(t@YJE*bONM)^YeE@RULP{%XMXLh2LEIVmR13}O zh38;;aSI|SS~!(;gJag9SfGfLtX{L=C@rENOKiC$Vl+TY)r)PWDy$=9JNi(sT&H~N z54#FOAanE`sG!Ee)HH)MEeauO+Vbt~va$n9`q`z8>0@A#0tgE|N*axj3!}2yU1rgQ zp#82{=HTl$HWXnv7!`p*ED;$9vI_->l0Y~tHZX-5$rDGacV`DHza>F6$jm#e3dkUY z=%*tgJCa}_(wHikJ|)GK1ReP51;#^~<7++d$>SmKV%Ww>BZY(G(x94W2xy$LYyydZ zn2L2KaH&i==aWC1+Kp}HueU=TF2oSSbr8JHIO*U>DNf{wC}=||qK9kYC@tKlBH#_9 zGnbc_Ycp`|{^^3RQYH%`#9~E}P=Lj8YZEwa6*GwD*bV^= zkL3nP5`#Ik`~;sMF;r=$6haYFL2n)u$V8J=4~|Fz@FZ97zZ2xGPx6jLlTRL8UB# zh+L31!UYYB=eSik#tma0hcY(6_VYEoJr4dFF_!{xPA2(-=Fv>}a~3@4;%iI_dCgnq zB2F}388Xsph$CDQ&ok^DUHTVrKogu3500q;xH;%3V}42uvD02``dWo2G6n`wZ|G^M zolu%if$mw(RFF*?rRMHXXm)SevSfFUYS{0*526iAtf2rRSOT@6Ct8G1vIu4kCfv4a z%v?dT^bW7rcKo#}l1)bLUM57S`L9}I$eRaYN}?u>y9SvzZ#uXelsF$EYHvf&89sy# z+CWH$gn6O7Q06uP#Ze45hUXM$hMk2M>uv&|ISg9yo4MKAIYoijbk7oIXM6iQjtE5+Yq=St*)S#X{Q$TOh|s9dh}xG**f{~P;cjwmoeq{7>A0t%2p0h&?J zUJwjee5iW6x6eMBl!M`EBYOrh)TyFGWL1QU!U&Q`5&$e92*G4oA~GPsA_z|Mg{15O zi1PGVzgBlMRFSl4UIIb1foLn*=(lp{kR$^0LQ#kSt_!S3UCA!tgJR9NnGLi=0*P{? zIL@xp5bT55?%2+Q(treXfXX?OgmVIfP8vseer|G;VaKy(&582$qpw^Mj9Lg(P*sS* zV#ll@DLQOEkcTG#W(_WMSU`ru;951=&<~6NWf*6=$FH#eqzie79-_{GO{Y&b~Z#ZLo7fX_X4Xj8x77PjnvM=fX&%l^k9^glBCq6nLSjCj1 zHaU9$RY6wRRVan4ZB(kNt(Zby-v3tfI2a5df&wuG39!^pu5001wr#)TiLK%0-lg!E8$ByWYIyO3=HQfU_tC5)0m@*_6x}Oh!Z0yAP~@l5M<%k(oFe%4Kre8 z(w6oo5i<@#W$?zmhluoot!AM(AWYkr1!P1-sbN7^MU{>O6)iM2v`|eumX~47^o5N% z+dfXw4SuAw%S@3Z(;c2V z7GAa{qWDlMf|`k~fk(qJ3M2V8Wv^jfMB^o*udfe4rmbj_CNeM>of#SVKYFc8GhKS# zE$u?fV>h686>Qx9?I?GM$U;cMGAIy|2&1}4k{vb-p)MN*&<7L-j2nglWovjTh2 m!a%|a05D7~Lp3h +#include +#include +#include +#include +#include +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::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::epsilon(); + const double FPMIN=numeric_limits::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;i2) q=gammq(0.5*(ndata-2),0.5*chi2); + } +} + +} /* end namespace cpl */ diff --git a/szamszim/cpl/datafit.hpp b/szamszim/cpl/datafit.hpp new file mode 100644 index 0000000..75d2219 --- /dev/null +++ b/szamszim/cpl/datafit.hpp @@ -0,0 +1,23 @@ +#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 */ diff --git a/szamszim/cpl/fft.cpp b/szamszim/cpl/fft.cpp new file mode 100644 index 0000000..d4a76a7 --- /dev/null +++ b/szamszim/cpl/fft.cpp @@ -0,0 +1,214 @@ +#include +#include +#include +#include +#include +using namespace std; + +#include "fft.hpp" +#include "vector.hpp" + +namespace cpl { + +ComplexVector::ComplexVector(int dim) { + v = new std::complex [this->dim = dim]; + for (int i = 0; i < dim; i++) v[i] = 0.0; +} + +ComplexVector::ComplexVector(const ComplexVector& cv) { + v = new std::complex [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 [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 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 W(0, pi / n); + W = inverse ? std::exp(W) : std::exp(-W); + std::complex W_j(1, 0); + for (int j = 0; j < n; ++j) { + for (int i = j; i < N; i += 2 * n) { + std::complex 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 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* [rows]; + for (int i = 0; i < rows; i++) { + m[i] = new complex [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* [rows]; + for (int i = 0; i < rows; i++) { + m[i] = new complex [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* ComplexMatrix::operator[](const int row) const { + if (row < 0 || row >= rows) + error("bad row index", row); + return m[row]; +} + +complex* 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* [rows]; + for (int i = 0; i < rows; i++) + m[i] = new complex [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 */ diff --git a/szamszim/cpl/fft.hpp b/szamszim/cpl/fft.hpp new file mode 100644 index 0000000..605dc69 --- /dev/null +++ b/szamszim/cpl/fft.hpp @@ -0,0 +1,75 @@ +#ifndef CPL_FFT_HPP +#define CPL_FFT_HPP + +#include +#include + +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 operator[](const int i) const { return v[i]; } + std::complex& operator[](const int i) { return v[i]; } + ComplexVector& operator = (const ComplexVector& cv); + + private: + int dim; + std::complex *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 d=0); + ComplexMatrix(const ComplexMatrix& m); + ~ComplexMatrix(); + + int numRows() const { return rows; } + int numCols() const { return cols; } + const std::complex* operator[](const int row) const; + std::complex* operator[](const int row); + std::complex& operator()(int row, int col); + ComplexMatrix& operator=(const ComplexMatrix& m); + + private: + int rows; + int cols; + std::complex **m; +}; + +class FFT2 { + public: + FFT2() { } + void transform(ComplexMatrix& data); + void inverseTransform(ComplexMatrix& data); + + private: + +}; + +} /* end namespace cpl */ + +#endif /* CPL_FFT_HPP */ diff --git a/szamszim/cpl/interpolate.cpp b/szamszim/cpl/interpolate.cpp new file mode 100644 index 0000000..08d68f4 --- /dev/null +++ b/szamszim/cpl/interpolate.cpp @@ -0,0 +1,59 @@ +#include +#include +#include +#include +#include +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& xa, // vector of abscissa values + const vector& 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 c(n),d(n); + dif=abs(x-xa[0]); + for (i=0;i& xa, // vector of x values - input + const vector& 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 */ diff --git a/szamszim/cpl/matrix.cpp b/szamszim/cpl/matrix.cpp new file mode 100644 index 0000000..a9c6a9d --- /dev/null +++ b/szamszim/cpl/matrix.cpp @@ -0,0 +1,895 @@ +#include +#include +#include +#include +#include +#include +#include + +#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 big) big=temp; + if (big == 0.0) error("Singular matrix in routine ludcmp"); + vv[i]=1.0/big; + } + for (j=0;j= big) { + big=dum; + imax=i; + } + } + if (j != imax) { + for (k=0;k=0;i--) { + sum=B[i][k]; + for (j=i+1;j 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=0;i--) { + if (i < n-1) { + if (g != 0.0) { + for (j=l;j=0;i--) { + l=i+1; + g=w[i]; + for (j=l;j=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::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 stpmax) + for (i=0;i= 0.0) error("Roundoff problem in lineSearch"); + test=0.0; + for (i=0;i test) test=temp; + } + alamin=TOLX/test; + alam=1.0; + for (;;) { + for (i=0;i 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::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 test) test=temp; + } + if (test < TOLX) + return; + for (i=0;i test) test=temp; + } + if (test < gtol) + return; + for (i=0;i std::sqrt(EPS*sumdg*sumxi)) { + fac=1.0/fac; + fad=1.0/fae; + for (i=0;i +#include +#include +#include +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 */ diff --git a/szamszim/cpl/odeint.hpp b/szamszim/cpl/odeint.hpp new file mode 100644 index 0000000..4192ef8 --- /dev/null +++ b/szamszim/cpl/odeint.hpp @@ -0,0 +1,44 @@ +#ifndef CPL_ODEINT_HPP +#define CPL_ODEINT_HPP + +#include + +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 */ diff --git a/szamszim/cpl/optimize.cpp b/szamszim/cpl/optimize.cpp new file mode 100644 index 0000000..872d08a --- /dev/null +++ b/szamszim/cpl/optimize.cpp @@ -0,0 +1,134 @@ +#include +#include +#include +#include +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 */ diff --git a/szamszim/cpl/optimize.hpp b/szamszim/cpl/optimize.hpp new file mode 100644 index 0000000..1a33a41 --- /dev/null +++ b/szamszim/cpl/optimize.hpp @@ -0,0 +1,47 @@ +#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 */ diff --git a/szamszim/cpl/random.cpp b/szamszim/cpl/random.cpp new file mode 100644 index 0000000..42b03d1 --- /dev/null +++ b/szamszim/cpl/random.cpp @@ -0,0 +1,73 @@ +#include +#include +#include +#include +#include +#include +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 */ diff --git a/szamszim/cpl/random.hpp b/szamszim/cpl/random.hpp new file mode 100644 index 0000000..1ecc36c --- /dev/null +++ b/szamszim/cpl/random.hpp @@ -0,0 +1,23 @@ +#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 */ diff --git a/szamszim/cpl/vector.cpp b/szamszim/cpl/vector.cpp new file mode 100644 index 0000000..b3f105a --- /dev/null +++ b/szamszim/cpl/vector.cpp @@ -0,0 +1,167 @@ +#include +#include +#include +#include +#include +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 */ diff --git a/szamszim/cpl/vector.hpp b/szamszim/cpl/vector.hpp new file mode 100644 index 0000000..73a4dd2 --- /dev/null +++ b/szamszim/cpl/vector.hpp @@ -0,0 +1,124 @@ +#ifndef CPL_VECTOR_HPP +#define CPL_VECTOR_HPP + +#include + +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 */ diff --git a/szamszim/f1/Makefile b/szamszim/f1/Makefile new file mode 100644 index 0000000..d9006b1 --- /dev/null +++ b/szamszim/f1/Makefile @@ -0,0 +1,30 @@ +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 diff --git a/szamszim/f1/chez.c b/szamszim/f1/chez.c new file mode 100644 index 0000000..f8165a1 --- /dev/null +++ b/szamszim/f1/chez.c @@ -0,0 +1,46 @@ +#include + +#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); +} diff --git a/szamszim/f1/euler-divergence.txt b/szamszim/f1/euler-divergence.txt new file mode 100644 index 0000000..8cb014b --- /dev/null +++ b/szamszim/f1/euler-divergence.txt @@ -0,0 +1,1001 @@ +0 1 0.1 50.005 +0.00628319 1.00063 -0.528319 50.2024 +0.0125664 0.997309 -1.15703 50.4006 +0.0188496 0.990039 -1.78366 50.5996 +0.0251327 0.978832 -2.40572 50.7993 +0.0314159 0.963716 -3.02074 50.9999 +0.0376991 0.944736 -3.62626 51.2012 +0.0439823 0.921952 -4.21985 51.4034 +0.0502655 0.895438 -4.79913 51.6063 +0.0565487 0.865284 -5.36175 51.81 +0.0628319 0.831595 -5.90543 52.0146 +0.069115 0.79449 -6.42793 52.2199 +0.0753982 0.754102 -6.92713 52.4261 +0.0816814 0.710578 -7.40094 52.633 +0.0879646 0.664076 -7.84741 52.8408 +0.0942478 0.61477 -8.26466 53.0494 +0.100531 0.562841 -8.65094 53.2589 +0.106814 0.508486 -9.00458 53.4691 +0.113097 0.451908 -9.32407 53.6802 +0.119381 0.393324 -9.60801 53.8921 +0.125664 0.332955 -9.85515 54.1049 +0.131947 0.271033 -10.0643 54.3185 +0.13823 0.207797 -10.2346 54.5329 +0.144513 0.143491 -10.3652 54.7482 +0.150796 0.0783641 -10.4554 54.9643 +0.15708 0.0126711 -10.5046 55.1813 +0.163363 -0.0533312 -10.5126 55.3992 +0.169646 -0.119384 -10.4791 55.6179 +0.175929 -0.185225 -10.404 55.8375 +0.182212 -0.250596 -10.2877 56.0579 +0.188496 -0.315235 -10.1302 56.2792 +0.194779 -0.378885 -9.93214 56.5014 +0.201062 -0.441291 -9.69408 56.7244 +0.207345 -0.5022 -9.41681 56.9484 +0.213628 -0.561368 -9.10126 57.1732 +0.219911 -0.618553 -8.74855 57.3989 +0.226195 -0.673522 -8.3599 57.6255 +0.232478 -0.726048 -7.93671 57.853 +0.238761 -0.775916 -7.48052 58.0814 +0.245044 -0.822918 -6.993 58.3107 +0.251327 -0.866856 -6.47595 58.5409 +0.257611 -0.907546 -5.93128 58.772 +0.263894 -0.944813 -5.36106 59.004 +0.270177 -0.978497 -4.76741 59.237 +0.27646 -1.00845 -4.15261 59.4708 +0.282743 -1.03454 -3.51898 59.7056 +0.289027 -1.05665 -2.86895 59.9413 +0.29531 -1.07468 -2.20504 60.178 +0.301593 -1.08853 -1.5298 60.4155 +0.307876 -1.09815 -0.84585 60.654 +0.314159 -1.10346 -0.155864 60.8935 +0.320442 -1.10444 0.537461 61.1339 +0.326726 -1.10106 1.2314 61.3752 +0.333009 -1.09333 1.92322 61.6175 +0.339292 -1.08124 2.61018 61.8608 +0.345575 -1.06484 3.28954 62.105 +0.351858 -1.04417 3.9586 62.3502 +0.358142 -1.0193 4.61468 62.5963 +0.364425 -0.990306 5.25512 62.8435 +0.370708 -0.957287 5.87735 63.0916 +0.376991 -0.920359 6.47883 63.3406 +0.383274 -0.879651 7.05711 63.5907 +0.389557 -0.83531 7.60981 63.8417 +0.395841 -0.787496 8.13465 64.0938 +0.402124 -0.736385 8.62945 64.3468 +0.408407 -0.682164 9.09213 64.6008 +0.41469 -0.625037 9.52075 64.8559 +0.420973 -0.565216 9.91347 65.1119 +0.427257 -0.502928 10.2686 65.369 +0.43354 -0.438408 10.5846 65.627 +0.439823 -0.371903 10.8601 65.8861 +0.446106 -0.303667 11.0937 66.1462 +0.452389 -0.233963 11.2845 66.4074 +0.458673 -0.16306 11.4315 66.6695 +0.464956 -0.0912339 11.534 66.9327 +0.471239 -0.0187637 11.5913 67.197 +0.477522 0.0540668 11.6031 67.4623 +0.483805 0.126971 11.5691 67.7286 +0.490088 0.199662 11.4894 67.996 +0.496372 0.271852 11.3639 68.2644 +0.502655 0.343254 11.1931 68.5339 +0.508938 0.413582 10.9774 68.8045 +0.515221 0.482555 10.7176 69.0761 +0.521504 0.549896 10.4144 69.3488 +0.527788 0.615331 10.0689 69.6226 +0.534071 0.678596 9.68223 69.8974 +0.540354 0.739431 9.25586 70.1734 +0.546637 0.797587 8.79126 70.4504 +0.55292 0.852824 8.29012 70.7285 +0.559203 0.904913 7.75428 71.0078 +0.565487 0.953634 7.18571 71.2881 +0.57177 0.998783 6.58652 71.5695 +0.578053 1.04017 5.95897 71.8521 +0.584336 1.07761 5.30541 72.1357 +0.590619 1.11094 4.62833 72.4205 +0.596903 1.14002 3.9303 72.7064 +0.603186 1.16472 3.214 72.9934 +0.609469 1.18491 2.48219 73.2816 +0.615752 1.20051 1.73768 73.5709 +0.622035 1.21143 0.983382 73.8614 +0.628319 1.21761 0.222219 74.153 +0.634602 1.219 -0.542825 74.4457 +0.640885 1.21559 -1.30875 74.7396 +0.647168 1.20737 -2.07253 75.0347 +0.653451 1.19435 -2.83114 75.3309 +0.659734 1.17656 -3.58157 75.6283 +0.666018 1.15405 -4.32082 75.9268 +0.672301 1.12691 -5.04594 76.2266 +0.678584 1.0952 -5.75399 76.5275 +0.684867 1.05905 -6.44213 76.8296 +0.69115 1.01857 -7.10755 77.133 +0.697434 0.973913 -7.74753 77.4375 +0.703717 0.925234 -8.35946 77.7432 +0.71 0.87271 -8.9408 78.0501 +0.716283 0.816533 -9.48914 78.3582 +0.722566 0.756911 -10.0022 78.6676 +0.728849 0.694065 -10.4778 78.9781 +0.735133 0.628232 -10.9139 79.2899 +0.741416 0.559658 -11.3086 79.603 +0.747699 0.488604 -11.6602 79.9172 +0.753982 0.41534 -11.9672 80.2327 +0.760265 0.340148 -12.2282 80.5495 +0.766549 0.263316 -12.4419 80.8675 +0.772832 0.185141 -12.6074 81.1867 +0.779115 0.105927 -12.7237 81.5072 +0.785398 0.0259814 -12.7902 81.829 +0.791681 -0.0543821 -12.8066 82.152 +0.797965 -0.134848 -12.7724 82.4764 +0.804248 -0.2151 -12.6877 82.802 +0.810531 -0.294819 -12.5525 83.1289 +0.816814 -0.373688 -12.3673 83.457 +0.823097 -0.451394 -12.1325 83.7865 +0.82938 -0.527625 -11.8489 84.1173 +0.835664 -0.602074 -11.5174 84.4494 +0.841947 -0.674439 -11.1391 84.7828 +0.84823 -0.744428 -10.7153 85.1175 +0.854513 -0.811754 -10.2476 85.4535 +0.860796 -0.876142 -9.73752 85.7909 +0.86708 -0.937324 -9.18702 86.1296 +0.873363 -0.995048 -8.59809 86.4696 +0.879646 -1.04907 -7.97288 86.8109 +0.885929 -1.09917 -7.31373 87.1537 +0.892212 -1.14512 -6.6231 87.4977 +0.898495 -1.18673 -5.9036 87.8432 +0.904779 -1.22383 -5.15795 88.1899 +0.911062 -1.25624 -4.389 88.5381 +0.917345 -1.28381 -3.59968 88.8876 +0.923628 -1.30643 -2.79304 89.2386 +0.929911 -1.32398 -1.97219 89.5909 +0.936195 -1.33637 -1.1403 89.9445 +0.942478 -1.34354 -0.300638 90.2996 +0.948761 -1.34542 0.543531 90.6561 +0.955044 -1.34201 1.38889 91.014 +0.961327 -1.33328 2.2321 91.3733 +0.967611 -1.31926 3.06982 91.7341 +0.973894 -1.29997 3.89874 92.0962 +0.980177 -1.27547 4.71553 92.4598 +0.98646 -1.24585 5.51694 92.8248 +0.992743 -1.21118 6.29972 93.1913 +0.999026 -1.1716 7.06073 93.5592 +1.00531 -1.12724 7.79687 93.9285 +1.01159 -1.07825 8.50513 94.2993 +1.01788 -1.02481 9.18261 94.6716 +1.02416 -0.967111 9.82652 95.0454 +1.03044 -0.905369 10.4342 95.4206 +1.03673 -0.839809 11.003 95.7973 +1.04301 -0.770675 11.5307 96.1755 +1.04929 -0.698225 12.0149 96.5552 +1.05558 -0.622733 12.4536 96.9364 +1.06186 -0.544485 12.8449 97.319 +1.06814 -0.463778 13.187 97.7032 +1.07442 -0.380921 13.4784 98.089 +1.08071 -0.296234 13.7178 98.4762 +1.08699 -0.210043 13.9039 98.865 +1.09327 -0.122682 14.0359 99.2553 +1.09956 -0.034492 14.1129 99.6471 +1.10584 0.0541823 14.1346 100.041 +1.11212 0.142993 14.1006 100.435 +1.11841 0.231589 14.0107 100.832 +1.12469 0.319621 13.8652 101.23 +1.13097 0.406739 13.6644 101.63 +1.13726 0.492595 13.4088 102.031 +1.14354 0.576845 13.0993 102.434 +1.14982 0.659151 12.7369 102.838 +1.15611 0.739179 12.3227 103.244 +1.16239 0.816605 11.8583 103.652 +1.16867 0.891113 11.3452 104.061 +1.17496 0.962397 10.7853 104.472 +1.18124 1.03016 10.1806 104.884 +1.18752 1.09413 9.53334 105.298 +1.19381 1.15403 8.84587 105.714 +1.20009 1.20961 8.12078 106.131 +1.20637 1.26063 7.36076 106.55 +1.21265 1.30688 6.56868 106.971 +1.21894 1.34815 5.74754 107.393 +1.22522 1.38427 4.90047 107.817 +1.2315 1.41506 4.03071 108.243 +1.23779 1.44038 3.1416 108.67 +1.24407 1.46012 2.23658 109.099 +1.25035 1.47418 1.31916 109.53 +1.25664 1.48246 0.392905 109.962 +1.26292 1.48493 -0.538555 110.396 +1.2692 1.48155 -1.47157 110.832 +1.27549 1.4723 -2.40245 111.27 +1.28177 1.45721 -3.32753 111.709 +1.28805 1.4363 -4.24312 112.15 +1.29434 1.40964 -5.14557 112.593 +1.30062 1.37731 -6.03127 113.037 +1.3069 1.33941 -6.89666 113.484 +1.31319 1.29608 -7.73824 113.932 +1.31947 1.24746 -8.55259 114.381 +1.32575 1.19372 -9.3364 114.833 +1.33204 1.13506 -10.0864 115.286 +1.33832 1.07169 -10.7996 115.741 +1.3446 1.00383 -11.473 116.198 +1.35088 0.931743 -12.1037 116.657 +1.35717 0.855693 -12.6891 117.118 +1.36345 0.775965 -13.2268 117.58 +1.36973 0.692858 -13.7143 118.044 +1.37602 0.606689 -14.1497 118.51 +1.3823 0.517784 -14.5309 118.978 +1.38858 0.426484 -14.8562 119.448 +1.39487 0.33314 -15.1242 119.919 +1.40115 0.238112 -15.3335 120.393 +1.40743 0.141769 -15.4831 120.868 +1.41372 0.0444854 -15.5722 121.345 +1.42 -0.0533574 -15.6001 121.824 +1.42628 -0.151376 -15.5666 122.305 +1.43257 -0.249184 -15.4715 122.788 +1.43885 -0.346394 -15.3149 123.273 +1.44513 -0.44262 -15.0973 123.759 +1.45142 -0.537479 -14.8192 124.248 +1.4577 -0.630591 -14.4815 124.738 +1.46398 -0.72158 -14.0852 125.231 +1.47027 -0.81008 -13.6319 125.725 +1.47655 -0.895732 -13.1229 126.222 +1.48283 -0.978185 -12.5601 126.72 +1.48911 -1.0571 -11.9455 127.22 +1.4954 -1.13216 -11.2813 127.722 +1.50168 -1.20304 -10.5699 128.227 +1.50796 -1.26945 -9.81401 128.733 +1.51425 -1.33112 -9.01639 129.241 +1.52053 -1.38777 -8.18002 129.751 +1.52681 -1.43916 -7.30806 130.264 +1.5331 -1.48508 -6.40381 130.778 +1.53938 -1.52532 -5.4707 131.294 +1.54566 -1.55969 -4.51232 131.812 +1.55195 -1.58804 -3.53233 132.333 +1.55823 -1.61024 -2.53454 132.855 +1.56451 -1.62616 -1.52279 133.38 +1.5708 -1.63573 -0.501045 133.906 +1.57708 -1.63888 0.526715 134.435 +1.58336 -1.63557 1.55645 134.966 +1.58965 -1.62579 2.58411 135.499 +1.59593 -1.60955 3.60563 136.033 +1.60221 -1.5869 4.61694 136.57 +1.6085 -1.55789 5.61402 137.11 +1.61478 -1.52262 6.59287 137.651 +1.62106 -1.48119 7.54956 138.194 +1.62734 -1.43376 8.48022 138.74 +1.63363 -1.38047 9.38107 139.288 +1.63991 -1.32153 10.2484 139.838 +1.64619 -1.25714 11.0788 140.39 +1.65248 -1.18753 11.8687 140.944 +1.65876 -1.11295 12.6148 141.5 +1.66504 -1.03369 13.3141 142.059 +1.67133 -0.950038 13.9636 142.62 +1.67761 -0.862302 14.5605 143.183 +1.68389 -0.770816 15.1023 143.748 +1.69018 -0.675925 15.5866 144.315 +1.69646 -0.577992 16.0113 144.885 +1.70274 -0.477389 16.3745 145.457 +1.70903 -0.374505 16.6745 146.031 +1.71531 -0.269737 16.9098 146.608 +1.72159 -0.163489 17.0792 147.187 +1.72788 -0.0561774 17.182 147.768 +1.73416 0.0517801 17.2173 148.351 +1.74044 0.159959 17.1847 148.937 +1.74673 0.267934 17.0842 149.525 +1.75301 0.375278 16.9159 150.115 +1.75929 0.481563 16.6801 150.708 +1.76558 0.586367 16.3775 151.303 +1.77186 0.68927 16.0091 151.9 +1.77814 0.789858 15.576 152.5 +1.78442 0.887725 15.0797 153.102 +1.79071 0.982474 14.5219 153.706 +1.79699 1.07372 13.9046 154.313 +1.80327 1.16108 13.23 154.922 +1.80956 1.24421 12.5005 155.534 +1.81584 1.32275 11.7187 156.148 +1.82212 1.39638 10.8876 156.764 +1.82841 1.46479 10.0102 157.383 +1.83469 1.52769 9.08987 158.004 +1.84097 1.5848 8.13 158.628 +1.84726 1.63588 7.13424 159.254 +1.85354 1.68071 6.10638 159.883 +1.85982 1.71908 5.05036 160.514 +1.86611 1.75081 3.97023 161.148 +1.87239 1.77576 2.87017 161.784 +1.87867 1.79379 1.75443 162.423 +1.88496 1.80481 0.627354 163.064 +1.89124 1.80875 -0.506643 163.708 +1.89752 1.80557 -1.64312 164.354 +1.90381 1.79525 -2.77759 165.003 +1.91009 1.77779 -3.90558 165.654 +1.91637 1.75326 -5.0226 166.308 +1.92265 1.7217 -6.1242 166.965 +1.92894 1.68322 -7.20598 167.624 +1.93522 1.63794 -8.26357 168.286 +1.9415 1.58602 -9.29272 168.95 +1.94779 1.52763 -10.2892 169.617 +1.95407 1.46298 -11.2491 170.287 +1.96035 1.3923 -12.1683 170.959 +1.96664 1.31585 -13.0431 171.634 +1.97292 1.23389 -13.8699 172.312 +1.9792 1.14675 -14.6452 172.992 +1.98549 1.05473 -15.3657 173.675 +1.99177 0.958184 -16.0284 174.361 +1.99805 0.857474 -16.6304 175.049 +2.00434 0.752982 -17.1692 175.74 +2.01062 0.645105 -17.6423 176.434 +2.0169 0.534255 -18.0476 177.13 +2.02319 0.420858 -18.3833 177.83 +2.02947 0.305352 -18.6478 178.532 +2.03575 0.188185 -18.8396 179.236 +2.04204 0.0698122 -18.9579 179.944 +2.04832 -0.0493036 -19.0017 180.654 +2.0546 -0.168695 -18.9707 181.368 +2.06088 -0.287892 -18.8648 182.084 +2.06717 -0.406422 -18.6839 182.802 +2.07345 -0.523817 -18.4285 183.524 +2.07973 -0.639606 -18.0994 184.249 +2.08602 -0.753328 -17.6975 184.976 +2.0923 -0.864525 -17.2242 185.706 +2.09858 -0.972747 -16.681 186.439 +2.10487 -1.07756 -16.0698 187.175 +2.11115 -1.17853 -15.3927 187.914 +2.11743 -1.27524 -14.6522 188.656 +2.12372 -1.3673 -13.851 189.401 +2.13 -1.45433 -12.9919 190.149 +2.13628 -1.53596 -12.0781 190.899 +2.14257 -1.61185 -11.113 191.653 +2.14885 -1.68168 -10.1003 192.41 +2.15513 -1.74514 -9.04364 193.169 +2.16142 -1.80196 -7.94713 193.932 +2.1677 -1.8519 -6.81493 194.697 +2.17398 -1.89472 -5.65135 195.466 +2.18027 -1.93022 -4.46086 196.238 +2.18655 -1.95825 -3.24807 197.013 +2.19283 -1.97866 -2.01766 197.79 +2.19911 -1.99134 -0.774433 198.571 +2.2054 -1.9962 0.476761 199.355 +2.21168 -1.99321 1.73101 200.142 +2.21796 -1.98233 2.98338 200.932 +2.22425 -1.96359 4.22892 201.725 +2.23053 -1.93702 5.46268 202.522 +2.23681 -1.90269 6.67974 203.321 +2.2431 -1.86072 7.87524 204.124 +2.24938 -1.81124 9.04436 204.93 +2.25566 -1.75441 10.1824 205.739 +2.26195 -1.69044 11.2847 206.551 +2.26823 -1.61953 12.3469 207.367 +2.27451 -1.54195 13.3644 208.185 +2.2808 -1.45798 14.3333 209.007 +2.28708 -1.36792 15.2494 209.832 +2.29336 -1.27211 16.1088 210.661 +2.29965 -1.17089 16.9081 211.492 +2.30593 -1.06466 17.6438 212.327 +2.31221 -0.953798 18.3128 213.165 +2.3185 -0.838735 18.9121 214.007 +2.32478 -0.719907 19.4391 214.852 +2.33106 -0.597768 19.8914 215.7 +2.33734 -0.472787 20.267 216.552 +2.34363 -0.345446 20.564 217.407 +2.34991 -0.216238 20.7811 218.265 +2.35619 -0.0856666 20.917 219.126 +2.36248 0.0457585 20.9708 219.992 +2.36876 0.177522 20.942 220.86 +2.37504 0.309104 20.8305 221.732 +2.38133 0.439986 20.6363 222.607 +2.38761 0.569648 20.3598 223.486 +2.39389 0.697572 20.0019 224.368 +2.40018 0.823248 19.5636 225.254 +2.40646 0.94617 19.0463 226.143 +2.41274 1.06584 18.4518 227.036 +2.41903 1.18178 17.7822 227.933 +2.42531 1.29351 17.0396 228.832 +2.43159 1.40057 16.2269 229.736 +2.43788 1.50253 15.3469 230.643 +2.44416 1.59895 14.4028 231.553 +2.45044 1.68945 13.3982 232.467 +2.45673 1.77363 12.3367 233.385 +2.46301 1.85115 11.2223 234.307 +2.46929 1.92166 10.0591 235.232 +2.47558 1.98486 8.85173 236.16 +2.48186 2.04048 7.60461 237.093 +2.48814 2.08826 6.32254 238.029 +2.49442 2.12798 5.01044 238.968 +2.50071 2.15947 3.67339 239.912 +2.50699 2.18255 2.31656 240.859 +2.51327 2.1971 0.945225 241.81 +2.51956 2.20304 -0.435255 242.764 +2.52584 2.20031 -1.81947 243.723 +2.53212 2.18887 -3.20196 244.685 +2.53841 2.16876 -4.57727 245.651 +2.54469 2.14 -5.93994 246.621 +2.55097 2.10267 -7.28454 247.594 +2.55726 2.0569 -8.60569 248.572 +2.56354 2.00283 -9.89808 249.553 +2.56982 1.94064 -11.1565 250.538 +2.57611 1.87054 -12.3758 251.527 +2.58239 1.79278 -13.5511 252.52 +2.58867 1.70764 -14.6776 253.517 +2.59496 1.61542 -15.7505 254.518 +2.60124 1.51645 -16.7655 255.523 +2.60752 1.41111 -17.7183 256.532 +2.61381 1.29979 -18.605 257.544 +2.62009 1.18289 -19.4216 258.561 +2.62637 1.06086 -20.1649 259.582 +2.63265 0.934158 -20.8314 260.607 +2.63894 0.80327 -21.4184 261.635 +2.64522 0.668695 -21.9231 262.668 +2.6515 0.530948 -22.3432 263.705 +2.65779 0.390561 -22.6768 264.746 +2.66407 0.248078 -22.9222 265.792 +2.67035 0.104054 -23.0781 266.841 +2.67664 -0.0409503 -23.1435 267.894 +2.68292 -0.186365 -23.1178 268.952 +2.6892 -0.331618 -23.0007 270.014 +2.69549 -0.476136 -22.7923 271.08 +2.70177 -0.619344 -22.4931 272.15 +2.70805 -0.760672 -22.104 273.224 +2.71434 -0.899556 -21.626 274.303 +2.72062 -1.03544 -21.0608 275.386 +2.7269 -1.16777 -20.4103 276.473 +2.73319 -1.29601 -19.6765 277.564 +2.73947 -1.41964 -18.8622 278.66 +2.74575 -1.53815 -17.9702 279.76 +2.75204 -1.65106 -17.0038 280.865 +2.75832 -1.7579 -15.9664 281.974 +2.7646 -1.85822 -14.8619 283.087 +2.77088 -1.9516 -13.6943 284.204 +2.77717 -2.03764 -12.4681 285.326 +2.78345 -2.11598 -11.1878 286.453 +2.78973 -2.18628 -9.85828 287.584 +2.79602 -2.24822 -8.4846 288.719 +2.8023 -2.30153 -7.07201 289.859 +2.80858 -2.34597 -5.62591 291.003 +2.81487 -2.38131 -4.1519 292.152 +2.82115 -2.4074 -2.65567 293.305 +2.82743 -2.42409 -1.14306 294.463 +2.83372 -2.43127 0.38004 295.626 +2.84 -2.42888 1.90765 296.793 +2.84628 -2.4169 3.43376 297.965 +2.85257 -2.39532 4.95234 299.141 +2.85885 -2.3642 6.45737 300.322 +2.86513 -2.32363 7.94284 301.507 +2.87142 -2.27372 9.40282 302.698 +2.8777 -2.21465 10.8314 303.893 +2.88398 -2.14659 12.2229 305.092 +2.89027 -2.06979 13.5717 306.297 +2.89655 -1.98452 14.8722 307.506 +2.90283 -1.89107 16.1191 308.72 +2.90911 -1.78979 17.3073 309.939 +2.9154 -1.68105 18.4318 311.162 +2.92168 -1.56524 19.4881 312.391 +2.92796 -1.44279 20.4715 313.624 +2.93425 -1.31416 21.3781 314.862 +2.94053 -1.17984 22.2038 316.105 +2.94681 -1.04033 22.9451 317.353 +2.9531 -0.896162 23.5988 318.606 +2.95938 -0.747887 24.1618 319.864 +2.96566 -0.596074 24.6317 321.127 +2.97195 -0.441308 25.0063 322.394 +2.97823 -0.284189 25.2836 323.667 +2.98451 -0.125328 25.4621 324.945 +2.9908 0.0346556 25.5409 326.228 +2.99708 0.195134 25.5191 327.516 +3.00336 0.355475 25.3965 328.809 +3.00965 0.515045 25.1731 330.107 +3.01593 0.673213 24.8495 331.41 +3.02221 0.829347 24.4265 332.718 +3.0285 0.982823 23.9054 334.032 +3.03478 1.13303 23.2879 335.351 +3.04106 1.27935 22.576 336.675 +3.04734 1.4212 21.7722 338.004 +3.05363 1.558 20.8792 339.338 +3.05991 1.68918 19.9003 340.678 +3.06619 1.81422 18.8389 342.023 +3.07248 1.93259 17.699 343.373 +3.07876 2.0438 16.4847 344.728 +3.08504 2.14737 15.2006 346.089 +3.09133 2.24288 13.8514 347.456 +3.09761 2.32991 12.4421 348.827 +3.10389 2.40809 10.9782 350.205 +3.11018 2.47707 9.46514 351.587 +3.11646 2.53654 7.90876 352.975 +3.12274 2.58623 6.31501 354.369 +3.12903 2.62591 4.69003 355.768 +3.13531 2.65538 3.04012 357.172 +3.14159 2.67448 1.3717 358.582 +3.14788 2.6831 -0.308721 359.998 +3.15416 2.68116 -1.99456 361.419 +3.16044 2.66862 -3.67918 362.846 +3.16673 2.64551 -5.35593 364.278 +3.17301 2.61185 -7.01815 365.716 +3.17929 2.56776 -8.65922 367.16 +3.18557 2.51335 -10.2726 368.61 +3.19186 2.44881 -11.8518 370.065 +3.19814 2.37434 -13.3904 371.526 +3.20442 2.2902 -14.8822 372.993 +3.21071 2.1967 -16.3212 374.465 +3.21699 2.09415 -17.7015 375.943 +3.22327 1.98293 -19.0172 377.428 +3.22956 1.86344 -20.2632 378.918 +3.23584 1.73612 -21.434 380.413 +3.24212 1.60145 -22.5248 381.915 +3.24841 1.45992 -23.531 383.423 +3.25469 1.31207 -24.4483 384.937 +3.26097 1.15846 -25.2727 386.456 +3.26726 0.999662 -26.0006 387.982 +3.27354 0.836295 -26.6287 389.514 +3.27982 0.668982 -27.1542 391.052 +3.28611 0.498368 -27.5745 392.595 +3.29239 0.325112 -27.8876 394.145 +3.29867 0.149889 -28.0919 395.701 +3.30496 -0.0266182 -28.1861 397.263 +3.31124 -0.203717 -28.1694 398.832 +3.31752 -0.38071 -28.0414 400.406 +3.32381 -0.556899 -27.8022 401.987 +3.33009 -0.731585 -27.4523 403.574 +3.33637 -0.904073 -26.9926 405.167 +3.34265 -1.07367 -26.4245 406.767 +3.34894 -1.2397 -25.7499 408.373 +3.35522 -1.40149 -24.971 409.985 +3.3615 -1.55839 -24.0904 411.603 +3.36779 -1.70976 -23.1113 413.228 +3.37407 -1.85497 -22.037 414.86 +3.38035 -1.99343 -20.8715 416.497 +3.38664 -2.12457 -19.619 418.142 +3.39292 -2.24784 -18.2841 419.792 +3.3992 -2.36272 -16.8717 421.45 +3.40549 -2.46873 -15.3872 423.114 +3.41177 -2.56541 -13.836 424.784 +3.41805 -2.65234 -12.2241 426.461 +3.42434 -2.72915 -10.5576 428.145 +3.43062 -2.79549 -8.84282 429.835 +3.4369 -2.85105 -7.08636 431.532 +3.44319 -2.89557 -5.29499 433.235 +3.44947 -2.92884 -3.47565 434.946 +3.45575 -2.95068 -1.63541 436.663 +3.46204 -2.96096 0.218559 438.387 +3.46832 -2.95958 2.07898 440.117 +3.4746 -2.94652 3.93854 441.855 +3.48088 -2.92177 5.7899 443.599 +3.48717 -2.88539 7.6257 445.35 +3.49345 -2.83748 9.43865 447.109 +3.49973 -2.77818 11.2215 448.874 +3.50602 -2.70767 12.9671 450.646 +3.5123 -2.62619 14.6683 452.425 +3.51858 -2.53403 16.3184 454.211 +3.52487 -2.4315 17.9106 456.004 +3.53115 -2.31896 19.4384 457.804 +3.53743 -2.19683 20.8954 459.612 +3.54372 -2.06554 22.2757 461.426 +3.55 -1.92558 23.5735 463.248 +3.55628 -1.77746 24.7834 465.077 +3.56257 -1.62174 25.9002 466.913 +3.56885 -1.459 26.9192 468.756 +3.57513 -1.28987 27.8359 470.607 +3.58142 -1.11497 28.6464 472.464 +3.5877 -0.934977 29.3469 474.33 +3.59398 -0.750585 29.9344 476.202 +3.60027 -0.562502 30.406 478.082 +3.60655 -0.371456 30.7594 479.97 +3.61283 -0.178189 30.9928 481.864 +3.61911 0.0165449 31.1048 483.767 +3.6254 0.211982 31.0944 485.677 +3.63168 0.407354 30.9612 487.594 +3.63796 0.601888 30.7052 489.519 +3.64425 0.794815 30.327 491.452 +3.65053 0.985365 29.8277 493.392 +3.65681 1.17278 29.2085 495.34 +3.6631 1.3563 28.4717 497.295 +3.66938 1.53519 27.6195 499.258 +3.67566 1.70873 26.6549 501.229 +3.68195 1.87621 25.5812 503.208 +3.68823 2.03694 24.4024 505.195 +3.69451 2.19027 23.1225 507.189 +3.7008 2.33555 21.7464 509.191 +3.70708 2.47219 20.2789 511.202 +3.71336 2.5996 18.7256 513.22 +3.71965 2.71726 17.0922 515.246 +3.72593 2.82465 15.3849 517.28 +3.73221 2.92132 13.6101 519.322 +3.7385 3.00683 11.7746 521.372 +3.74478 3.08081 9.88534 523.431 +3.75106 3.14292 7.94961 525.497 +3.75734 3.19287 5.97485 527.572 +3.76363 3.23041 3.96871 529.654 +3.76991 3.25535 1.93898 531.745 +3.77619 3.26753 -0.106419 533.845 +3.78248 3.26687 -2.15947 535.952 +3.78876 3.2533 -4.2121 538.068 +3.79504 3.22683 -6.25621 540.192 +3.80133 3.18752 -8.28369 542.325 +3.80761 3.13547 -10.2865 544.466 +3.81389 3.07084 -12.2565 546.615 +3.82018 2.99383 -14.186 548.773 +3.82646 2.9047 -16.0671 550.94 +3.83274 2.80375 -17.8922 553.115 +3.83903 2.69133 -19.6538 555.298 +3.84531 2.56784 -21.3448 557.491 +3.85159 2.43373 -22.9582 559.691 +3.85788 2.28947 -24.4874 561.901 +3.86416 2.13562 -25.9259 564.119 +3.87044 1.97272 -27.2678 566.346 +3.87673 1.80139 -28.5073 568.582 +3.88301 1.62227 -29.6391 570.827 +3.88929 1.43605 -30.6584 573.08 +3.89557 1.24341 -31.5607 575.343 +3.90186 1.04511 -32.342 577.614 +3.90814 0.8419 -32.9986 579.895 +3.91442 0.634564 -33.5276 582.184 +3.92071 0.423904 -33.9263 584.482 +3.92699 0.210738 -34.1927 586.79 +3.93327 -0.00410059 -34.3251 589.106 +3.93956 -0.219771 -34.3225 591.432 +3.94584 -0.435426 -34.1844 593.767 +3.95212 -0.650213 -33.9108 596.111 +3.95841 -0.863281 -33.5023 598.464 +3.96469 -1.07378 -32.9599 600.827 +3.97097 -1.28088 -32.2852 603.199 +3.97726 -1.48373 -31.4804 605.58 +3.98354 -1.68153 -30.5481 607.971 +3.98982 -1.87347 -29.4916 610.371 +3.99611 -2.05877 -28.3145 612.781 +4.00239 -2.23667 -27.0209 615.2 +4.00867 -2.40645 -25.6156 617.629 +4.01496 -2.5674 -24.1036 620.067 +4.02124 -2.71884 -22.4904 622.515 +4.02752 -2.86016 -20.7821 624.972 +4.0338 -2.99073 -18.985 627.44 +4.04009 -3.11002 -17.1059 629.917 +4.04637 -3.2175 -15.1518 632.404 +4.05265 -3.3127 -13.1302 634.9 +4.05894 -3.3952 -11.0488 637.407 +4.06522 -3.46462 -8.91549 639.923 +4.0715 -3.52064 -6.73861 642.449 +4.07779 -3.56298 -4.52652 644.986 +4.08407 -3.59142 -2.28784 647.532 +4.09035 -3.60579 -0.0312827 650.088 +4.09664 -3.60599 2.2343 652.655 +4.10292 -3.59195 4.50002 655.231 +4.1092 -3.56368 6.75691 657.818 +4.11549 -3.52122 8.99603 660.415 +4.12177 -3.4647 11.2085 663.022 +4.12805 -3.39427 13.3854 665.64 +4.13434 -3.31017 15.5181 668.268 +4.14062 -3.21267 17.5979 670.906 +4.1469 -3.1021 19.6165 673.554 +4.15319 -2.97884 21.5656 676.214 +4.15947 -2.84334 23.4373 678.883 +4.16575 -2.69608 25.2238 681.563 +4.17204 -2.5376 26.9178 684.254 +4.17832 -2.36847 28.5122 686.955 +4.1846 -2.18932 30.0004 689.667 +4.19088 -2.00082 31.376 692.39 +4.19717 -1.80368 32.6331 695.123 +4.20345 -1.59864 33.7664 697.868 +4.20973 -1.38648 34.7709 700.623 +4.21602 -1.16801 35.642 703.389 +4.2223 -0.944062 36.3759 706.166 +4.22858 -0.715505 36.9691 708.953 +4.23487 -0.483222 37.4186 711.752 +4.24115 -0.248113 37.7223 714.562 +4.24743 -0.0110975 37.8781 717.383 +4.25372 0.226898 37.8851 720.215 +4.26 0.464937 37.7426 723.058 +4.26628 0.702081 37.4504 725.913 +4.27257 0.937389 37.0093 728.779 +4.27885 1.16992 36.4203 731.656 +4.28513 1.39876 35.6852 734.544 +4.29142 1.62298 34.8064 737.444 +4.2977 1.84167 33.7866 740.356 +4.30398 2.05396 32.6295 743.278 +4.31027 2.25898 31.3389 746.213 +4.31655 2.45588 29.9196 749.159 +4.32283 2.64388 28.3765 752.116 +4.32911 2.82217 26.7153 755.085 +4.3354 2.99003 24.9421 758.066 +4.34168 3.14674 23.0634 761.059 +4.34796 3.29165 21.0862 764.064 +4.35425 3.42414 19.018 767.08 +4.36053 3.54364 16.8666 770.108 +4.36681 3.64961 14.64 773.149 +4.3731 3.7416 12.3469 776.201 +4.37938 3.81918 9.99599 779.265 +4.38566 3.88198 7.59634 782.342 +4.39195 3.92971 5.15721 785.43 +4.39823 3.96212 2.6881 788.531 +4.40451 3.97901 0.198633 791.644 +4.4108 3.98025 -2.30145 794.769 +4.41708 3.96579 -4.80232 797.907 +4.42336 3.93562 -7.2941 801.057 +4.42965 3.88979 -9.76692 804.219 +4.43593 3.82842 -12.2109 807.394 +4.44221 3.7517 -14.6164 810.582 +4.4485 3.65986 -16.9737 813.782 +4.45478 3.55321 -19.2732 816.994 +4.46106 3.43211 -21.5058 820.22 +4.46734 3.29699 -23.6622 823.458 +4.47363 3.14832 -25.7338 826.709 +4.47991 2.98662 -27.7119 829.972 +4.48619 2.81251 -29.5885 833.249 +4.49248 2.6266 -31.3557 836.539 +4.49876 2.42958 -33.006 839.841 +4.50504 2.2222 -34.5325 843.157 +4.51133 2.00523 -35.9288 846.485 +4.51761 1.77948 -37.1887 849.827 +4.52389 1.54581 -38.3068 853.182 +4.53018 1.30513 -39.2781 856.55 +4.53646 1.05833 -40.0981 859.932 +4.54274 0.806391 -40.7631 863.327 +4.54903 0.550269 -41.2697 866.735 +4.55531 0.290963 -41.6155 870.157 +4.56159 0.0294857 -41.7983 873.592 +4.56788 -0.233141 -41.8168 877.041 +4.57416 -0.495883 -41.6703 880.503 +4.58044 -0.757706 -41.3588 883.979 +4.58673 -1.01757 -40.8827 887.469 +4.59301 -1.27444 -40.2433 890.973 +4.59929 -1.5273 -39.4426 894.49 +4.60557 -1.77513 -38.4829 898.021 +4.61186 -2.01692 -37.3676 901.567 +4.61814 -2.25171 -36.1003 905.126 +4.62442 -2.47853 -34.6855 908.699 +4.63071 -2.69647 -33.1282 912.287 +4.63699 -2.90462 -31.434 915.888 +4.64327 -3.10212 -29.609 919.504 +4.64956 -3.28816 -27.6598 923.134 +4.65584 -3.46196 -25.5938 926.778 +4.66212 -3.62277 -23.4186 930.437 +4.66841 -3.76991 -21.1424 934.11 +4.67469 -3.90275 -18.7737 937.798 +4.68097 -4.02071 -16.3215 941.5 +4.68726 -4.12326 -13.7952 945.217 +4.69354 -4.20994 -11.2045 948.949 +4.69982 -4.28034 -8.55929 952.695 +4.70611 -4.33412 -5.86988 956.456 +4.71239 -4.371 -3.14667 960.232 +4.71867 -4.39077 -0.40029 964.023 +4.72496 -4.39328 2.35851 967.829 +4.73124 -4.37847 5.11889 971.65 +4.73752 -4.3463 7.86997 975.486 +4.7438 -4.29685 10.6008 979.337 +4.75009 -4.23025 13.3006 983.203 +4.75637 -4.14668 15.9586 987.084 +4.76265 -4.04641 18.564 990.981 +4.76894 -3.92977 21.1064 994.894 +4.77522 -3.79715 23.5756 998.821 +4.7815 -3.64902 25.9614 1002.76 +4.78779 -3.4859 28.2541 1006.72 +4.79407 -3.30837 30.4444 1010.7 +4.80035 -3.11709 32.5231 1014.69 +4.80664 -2.91274 34.4816 1018.69 +4.81292 -2.69608 36.3118 1022.72 +4.8192 -2.46793 38.0058 1026.75 +4.82549 -2.22913 39.5564 1030.81 +4.83177 -1.98059 40.957 1034.88 +4.83805 -1.72325 42.2015 1038.96 +4.84434 -1.45809 43.2842 1043.06 +4.85062 -1.18613 44.2003 1047.18 +4.8569 -0.90841 44.9456 1051.31 +4.86319 -0.626009 45.5164 1055.47 +4.86947 -0.340021 45.9097 1059.63 +4.87575 -0.0515615 46.1234 1063.82 +4.88203 0.23824 46.1558 1068.01 +4.88832 0.528245 46.0061 1072.23 +4.8946 0.81731 45.6742 1076.46 +4.90088 1.10429 45.1606 1080.71 +4.90717 1.38804 44.4668 1084.98 +4.91345 1.66743 43.5947 1089.26 +4.91973 1.94135 42.547 1093.56 +4.92602 2.20868 41.3272 1097.88 +4.9323 2.46834 39.9394 1102.22 +4.93858 2.71929 38.3885 1106.57 +4.94487 2.96049 36.6799 1110.94 +4.95115 3.19096 34.8198 1115.32 +4.95743 3.40974 32.8149 1119.72 +4.96372 3.61592 30.6725 1124.14 +4.97 3.80864 28.4005 1128.58 +4.97628 3.98709 26.0075 1133.04 +4.98257 4.1505 23.5023 1137.51 +4.98885 4.29817 20.8945 1142 +4.99513 4.42945 18.1939 1146.51 +5.00142 4.54377 15.4108 1151.04 +5.0077 4.6406 12.5558 1155.58 +5.01398 4.71949 9.64005 1160.14 +5.02027 4.78006 6.67471 1164.72 +5.02655 4.822 3.67131 1169.32 +5.03283 4.84506 0.641562 1173.94 +5.03911 4.84909 -2.40268 1178.57 +5.0454 4.834 -5.44946 1183.22 +5.05168 4.79976 -8.48675 1187.9 +5.05796 4.74643 -11.5025 1192.59 +5.06425 4.67416 -14.4848 1197.29 +5.07053 4.58315 -17.4217 1202.02 +5.07681 4.47369 -20.3013 1206.77 +5.0831 4.34613 -23.1122 1211.53 +5.08938 4.20091 -25.843 1216.31 +5.09566 4.03853 -28.4825 1221.11 +5.10195 3.85957 -31.02 1225.94 +5.10823 3.66467 -33.445 1230.78 +5.11451 3.45453 -35.7476 1235.63 +5.1208 3.22992 -37.9182 1240.51 +5.12708 2.99167 -39.9476 1245.41 +5.13336 2.74067 -41.8273 1250.33 +5.13965 2.47787 -43.5493 1255.26 +5.14593 2.20424 -45.1062 1260.22 +5.15221 1.92083 -46.4912 1265.19 +5.1585 1.62871 -47.6981 1270.19 +5.16478 1.32902 -48.7214 1275.2 +5.17106 1.02289 -49.5565 1280.24 +5.17734 0.711521 -50.1992 1285.29 +5.18363 0.39611 -50.6462 1290.36 +5.18991 0.0778908 -50.8951 1295.46 +5.19619 -0.241893 -50.944 1300.57 +5.20248 -0.561983 -50.7921 1305.71 +5.20876 -0.881119 -50.439 1310.86 +5.21504 -1.19804 -49.8853 1316.04 +5.22133 -1.51148 -49.1326 1321.23 +5.22761 -1.82018 -48.1829 1326.45 +5.23389 -2.12293 -47.0392 1331.69 +5.24018 -2.41848 -45.7054 1336.94 +5.24646 -2.70566 -44.1858 1342.22 +5.25274 -2.98329 -42.4858 1347.52 +5.25903 -3.25023 -40.6113 1352.84 +5.26531 -3.5054 -38.5691 1358.18 +5.27159 -3.74774 -36.3666 1363.54 +5.27788 -3.97624 -34.0119 1368.93 +5.28416 -4.18994 -31.5135 1374.33 +5.29044 -4.38794 -28.8809 1379.76 +5.29673 -4.56941 -26.1239 1385.2 +5.30301 -4.73355 -23.2528 1390.67 +5.30929 -4.87965 -20.2787 1396.16 +5.31557 -5.00706 -17.2127 1401.67 +5.32186 -5.11522 -14.0666 1407.21 +5.32814 -5.2036 -10.8527 1412.76 +5.33442 -5.27179 -7.58314 1418.34 +5.34071 -5.31943 -4.27078 1423.94 +5.34699 -5.34627 -0.928483 1429.56 +5.35327 -5.3521 2.43068 1435.2 +5.35956 -5.33683 5.7935 1440.87 +5.36584 -5.30043 9.14673 1446.56 +5.37212 -5.24296 12.4771 1452.27 +5.37841 -5.16456 15.7713 1458 +5.38469 -5.06547 19.0163 1463.76 +5.39097 -4.94598 22.1991 1469.54 +5.39726 -4.8065 25.3067 1475.34 +5.40354 -4.6475 28.3267 1481.16 +5.40982 -4.46951 31.2468 1487.01 +5.41611 -4.27319 34.0551 1492.88 +5.42239 -4.05921 36.74 1498.77 +5.42867 -3.82837 39.2905 1504.69 +5.43496 -3.5815 41.6959 1510.63 +5.44124 -3.31951 43.9463 1516.6 +5.44752 -3.04339 46.032 1522.58 +5.4538 -2.75416 47.9442 1528.59 +5.46009 -2.45292 49.6747 1534.63 +5.46637 -2.14081 51.2159 1540.69 +5.47265 -1.81901 52.561 1546.77 +5.47894 -1.48876 53.7039 1552.88 +5.48522 -1.15132 54.6393 1559.01 +5.4915 -0.808016 55.3627 1565.16 +5.49779 -0.460161 55.8704 1571.34 +5.50407 -0.109117 56.1596 1577.54 +5.51035 0.243744 56.2281 1583.77 +5.51664 0.597035 56.075 1590.02 +5.52292 0.949365 55.6998 1596.3 +5.5292 1.29934 55.1033 1602.6 +5.53549 1.64556 54.2869 1608.93 +5.54177 1.98666 53.253 1615.28 +5.54805 2.32126 52.0047 1621.66 +5.55434 2.64801 50.5463 1628.06 +5.56062 2.9656 48.8825 1634.49 +5.5669 3.27274 47.0191 1640.94 +5.57319 3.56817 44.9628 1647.42 +5.57947 3.85068 42.7209 1653.92 +5.58575 4.1191 40.3014 1660.45 +5.59203 4.37232 37.7133 1667.01 +5.59832 4.60928 34.9661 1673.59 +5.6046 4.82898 32.07 1680.19 +5.61088 5.03048 29.0358 1686.83 +5.61717 5.21292 25.8751 1693.49 +5.62345 5.3755 22.5997 1700.17 +5.62973 5.5175 19.2222 1706.88 +5.63602 5.63827 15.7555 1713.62 +5.6423 5.73727 12.2128 1720.39 +5.64858 5.814 8.60799 1727.18 +5.65487 5.86809 4.95494 1734 +5.66115 5.89922 1.26791 1740.84 +5.66743 5.90719 -2.43868 1747.72 +5.67372 5.89187 -6.15027 1754.62 +5.68 5.85322 -9.85224 1761.54 +5.68628 5.79132 -13.5299 1768.5 +5.69257 5.70631 -17.1687 1775.48 +5.69885 5.59843 -20.7541 1782.49 +5.70513 5.46803 -24.2717 1789.53 +5.71142 5.31553 -27.7074 1796.59 +5.7177 5.14144 -31.0472 1803.68 +5.72398 4.94636 -34.2777 1810.8 +5.73027 4.73099 -37.3856 1817.95 +5.73655 4.49609 -40.3581 1825.13 +5.74283 4.24251 -43.1831 1832.34 +5.74911 3.97118 -45.8488 1839.57 +5.7554 3.68311 -48.3439 1846.83 +5.76168 3.37935 -50.6581 1854.12 +5.76796 3.06106 -52.7814 1861.44 +5.77425 2.72942 -54.7047 1868.79 +5.78053 2.3857 -56.4197 1876.17 +5.78681 2.03121 -57.9186 1883.58 +5.7931 1.6673 -59.1949 1891.01 +5.79938 1.29536 -60.2425 1898.48 +5.80566 0.916848 -61.0564 1905.97 +5.81195 0.53322 -61.6325 1913.5 +5.81823 0.145972 -61.9675 1921.05 +5.82451 -0.243382 -62.0592 1928.63 +5.8308 -0.633311 -61.9063 1936.25 +5.83708 -1.02228 -61.5084 1943.89 +5.84336 -1.40875 -60.866 1951.57 +5.84965 -1.79118 -59.9809 1959.27 +5.85593 -2.16805 -58.8555 1967.01 +5.86221 -2.53785 -57.4932 1974.77 +5.8685 -2.89909 -55.8987 1982.57 +5.87478 -3.25031 -54.0771 1990.39 +5.88106 -3.59009 -52.0349 1998.25 +5.88734 -3.91704 -49.7792 2006.14 +5.89363 -4.22981 -47.318 2014.06 +5.89991 -4.52712 -44.6603 2022.01 +5.90619 -4.80772 -41.8159 2029.99 +5.91248 -5.07046 -38.7951 2038.01 +5.91876 -5.31422 -35.6092 2046.05 +5.92504 -5.53796 -32.2702 2054.13 +5.93133 -5.74072 -28.7906 2062.24 +5.93761 -5.92161 -25.1836 2070.38 +5.94389 -6.07985 -21.4629 2078.56 +5.95018 -6.2147 -17.6429 2086.76 +5.95646 -6.32556 -13.7381 2095 +5.96274 -6.41187 -9.76359 2103.27 +5.96903 -6.47322 -5.73489 2111.57 +5.97531 -6.50925 -1.66765 2119.91 +5.98159 -6.51973 2.42224 2128.28 +5.98788 -6.50451 6.51871 2136.68 +5.99416 -6.46356 10.6056 2145.12 +6.00044 -6.39692 14.6668 2153.59 +6.00673 -6.30476 18.6861 2162.09 +6.01301 -6.18736 22.6475 2170.62 +6.01929 -6.04506 26.5351 2179.19 +6.02557 -5.87833 30.3333 2187.8 +6.03186 -5.68774 34.0268 2196.43 +6.03814 -5.47395 37.6005 2205.1 +6.04442 -5.23769 41.0399 2213.81 +6.05071 -4.97983 44.3308 2222.55 +6.05699 -4.70129 47.4598 2231.32 +6.06327 -4.4031 50.4137 2240.13 +6.06956 -4.08634 53.1802 2248.98 +6.07584 -3.7522 55.7477 2257.85 +6.08212 -3.40192 58.1053 2266.77 +6.08841 -3.03684 60.2428 2275.72 +6.09469 -2.65832 62.1509 2284.7 +6.10097 -2.26781 63.8212 2293.72 +6.10726 -1.86681 65.2461 2302.78 +6.11354 -1.45686 66.419 2311.87 +6.11982 -1.03954 67.3344 2320.99 +6.12611 -0.616463 67.9876 2330.16 +6.13239 -0.189284 68.3749 2339.36 +6.13867 0.240328 68.4938 2348.59 +6.14496 0.670688 68.3428 2357.86 +6.15124 1.1001 67.9214 2367.17 +6.15752 1.52686 67.2302 2376.52 +6.1638 1.94928 66.2709 2385.9 +6.17009 2.36567 65.0461 2395.32 +6.17637 2.77437 63.5597 2404.77 +6.18265 3.17373 61.8165 2414.27 +6.18894 3.56213 59.8224 2423.8 +6.19522 3.93801 57.5842 2433.37 +6.2015 4.29982 55.1099 2442.97 +6.20779 4.64609 52.4083 2452.62 +6.21407 4.97538 49.489 2462.3 +6.22035 5.28633 46.3629 2472.02 +6.22664 5.57763 43.0414 2481.78 +6.23292 5.84807 39.5369 2491.58 +6.2392 6.09649 35.8624 2501.42 +6.24549 6.32182 32.0319 2511.29 +6.25177 6.52308 28.0598 2521.2 +6.25805 6.69938 23.9612 2531.16 +6.26434 6.84994 19.7519 2541.15 +6.27062 6.97404 15.4479 2551.18 +6.2769 7.0711 11.066 2561.25 +6.28319 7.14063 6.62311 2571.37 diff --git a/szamszim/f1/long.txt b/szamszim/f1/long.txt new file mode 100644 index 0000000..297bffc --- /dev/null +++ b/szamszim/f1/long.txt @@ -0,0 +1,1001 @@ +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 +1.26292 0.996504 -0.548973 49.8017 +1.2692 0.989121 -1.1751 49.6084 +1.27549 0.977833 -1.79658 49.4217 +1.28177 0.962684 -2.41097 49.2444 +1.28805 0.943735 -3.01584 49.0794 +1.29434 0.92106 -3.60881 48.9293 +1.30062 0.894749 -4.18753 48.7965 +1.3069 0.864906 -4.74971 48.683 +1.31319 0.831648 -5.29315 48.5906 +1.31947 0.795107 -5.81569 48.5209 +1.32575 0.755427 -6.31527 48.4748 +1.33204 0.712765 -6.78992 48.4532 +1.33832 0.667288 -7.23776 48.4563 +1.3446 0.619178 -7.65703 48.4841 +1.35088 0.568623 -8.04607 48.5362 +1.35717 0.515823 -8.40335 48.6118 +1.36345 0.460987 -8.72745 48.7096 +1.36973 0.404331 -9.0171 48.8282 +1.37602 0.346078 -9.27115 48.9656 +1.3823 0.28646 -9.48859 49.1197 +1.38858 0.22571 -9.66858 49.288 +1.39487 0.16407 -9.8104 49.4679 +1.40115 0.101782 -9.91349 49.6566 +1.40743 0.0390915 -9.97744 49.8511 +1.41372 -0.0237529 -10.002 50.0482 +1.42 -0.0865036 -9.98708 50.245 +1.42628 -0.148913 -9.93272 50.4383 +1.43257 -0.210734 -9.83916 50.625 +1.43885 -0.271723 -9.70675 50.8022 +1.44513 -0.33164 -9.53602 50.9671 +1.45142 -0.390247 -9.32765 51.1172 +1.4577 -0.447314 -9.08245 51.2499 +1.46398 -0.502615 -8.80139 51.3633 +1.47027 -0.555931 -8.48559 51.4556 +1.47655 -0.607053 -8.13629 51.5253 +1.48283 -0.655778 -7.75487 51.5712 +1.48911 -0.701915 -7.34283 51.5928 +1.4954 -0.74528 -6.9018 51.5895 +1.50168 -0.785703 -6.43353 51.5616 +1.50796 -0.823024 -5.93986 51.5094 +1.51425 -0.857096 -5.42274 51.4337 +1.52053 -0.887785 -4.88421 51.3358 +1.52681 -0.914968 -4.32639 51.2172 +1.5331 -0.93854 -3.7515 51.0797 +1.53938 -0.958406 -3.1618 50.9256 +1.54566 -0.974488 -2.55962 50.7572 +1.55195 -0.986724 -1.94733 50.5772 +1.55823 -0.995064 -1.32735 50.3885 +1.56451 -0.999475 -0.702134 50.1941 +1.5708 -0.999941 -0.074145 49.9969 +1.57708 -0.99646 0.554137 49.8001 +1.58336 -0.989044 1.18023 49.6069 +1.58965 -0.977724 1.80167 49.4202 +1.59593 -0.962544 2.41599 49.243 +1.60221 -0.943564 3.02077 49.0781 +1.6085 -0.920858 3.61363 48.9282 +1.61478 -0.894518 4.19222 48.7955 +1.62106 -0.864646 4.75426 48.6822 +1.62734 -0.831361 5.29754 48.59 +1.63363 -0.794793 5.8199 48.5204 +1.63991 -0.755088 6.31928 48.4745 +1.64619 -0.712402 6.79372 48.4531 +1.65248 -0.666903 7.24133 48.4564 +1.65876 -0.618772 7.66036 48.4845 +1.66504 -0.568197 8.04914 48.5368 +1.67133 -0.51538 8.40615 48.6125 +1.67761 -0.460528 8.72998 48.7105 +1.68389 -0.403858 9.01933 48.8293 +1.69018 -0.345593 9.27309 48.9668 +1.69646 -0.285964 9.49023 49.121 +1.70274 -0.225207 9.6699 49.2894 +1.70903 -0.16356 9.81141 49.4694 +1.71531 -0.101267 9.91417 49.6582 +1.72159 -0.0385747 9.9778 49.8527 +1.72788 0.0242699 10.002 50.0498 +1.73416 0.0870188 9.98679 50.2466 +1.74044 0.149424 9.93211 50.4398 +1.74673 0.211239 9.83823 50.6265 +1.75301 0.272221 9.7055 50.8036 +1.75929 0.332128 9.53446 50.9684 +1.76558 0.390723 9.32578 51.1183 +1.77186 0.447776 9.08028 51.2509 +1.77814 0.503062 8.79893 51.3642 +1.78442 0.556361 8.48285 51.4563 +1.79071 0.607464 8.13328 51.5257 +1.79699 0.656169 7.7516 51.5715 +1.80327 0.702283 7.33932 51.5928 +1.80956 0.745625 6.89806 51.5894 +1.81584 0.786023 6.42957 51.5613 +1.82212 0.823318 5.9357 51.5089 +1.82841 0.857363 5.41839 51.433 +1.83469 0.888023 4.87969 51.3349 +1.84097 0.915177 4.32173 51.2161 +1.84726 0.938718 3.74671 51.0785 +1.85354 0.958554 3.15689 50.9242 +1.85982 0.974605 2.55462 50.7558 +1.86611 0.986808 1.94226 50.5757 +1.87239 0.995116 1.32223 50.3869 +1.87867 0.999495 0.696975 50.1924 +1.88496 0.999929 0.0689739 49.9953 +1.89124 0.996415 -0.5593 49.7985 +1.89752 0.988967 -1.18537 49.6053 +1.90381 0.977615 -1.80675 49.4187 +1.91009 0.962403 -2.421 49.2416 +1.91637 0.943392 -3.0257 49.0768 +1.92265 0.920656 -3.61845 48.927 +1.92894 0.894286 -4.19692 48.7945 +1.93522 0.864386 -4.75881 48.6813 +1.9415 0.831073 -5.30192 48.5893 +1.94779 0.794479 -5.8241 48.5199 +1.95407 0.754749 -6.32329 48.4743 +1.96035 0.712039 -6.79751 48.453 +1.96664 0.666518 -7.2449 48.4566 +1.97292 0.618365 -7.66368 48.4848 +1.9792 0.567772 -8.05221 48.5373 +1.98549 0.514937 -8.40895 48.6132 +1.99177 0.460069 -8.7325 48.7114 +1.99805 0.403385 -9.02157 48.8303 +2.00434 0.345108 -9.27502 48.968 +2.01062 0.285469 -9.49186 49.1223 +2.0169 0.224703 -9.67123 49.2909 +2.02319 0.16305 -9.81241 49.471 +2.02947 0.100753 -9.91486 49.6598 +2.03575 0.038058 -9.97816 49.8543 +2.04204 -0.0247869 -10.0021 50.0515 +2.04832 -0.0875339 -9.9865 50.2482 +2.0546 -0.149935 -9.9315 50.4414 +2.06088 -0.211745 -9.83729 50.628 +2.06717 -0.272719 -9.70425 50.805 +2.07345 -0.332615 -9.5329 50.9697 +2.07973 -0.391199 -9.32391 51.1195 +2.08602 -0.448239 -9.07811 51.2519 +2.0923 -0.503509 -8.79647 51.365 +2.09858 -0.556791 -8.48011 51.4569 +2.10487 -0.607875 -8.13027 51.5262 +2.11115 -0.656559 -7.74833 51.5718 +2.11743 -0.702651 -7.3358 51.5929 +2.12372 -0.745969 -6.89431 51.5893 +2.13 -0.786343 -6.42561 51.5609 +2.13628 -0.823612 -5.93153 51.5083 +2.14257 -0.857629 -5.41404 51.4323 +2.14885 -0.888261 -4.87518 51.334 +2.15513 -0.915386 -4.31707 51.2151 +2.16142 -0.938897 -3.74191 51.0773 +2.1677 -0.958701 -3.15199 50.9229 +2.17398 -0.974721 -2.54962 50.7543 +2.18027 -0.986893 -1.93718 50.5742 +2.18655 -0.995168 -1.3171 50.3854 +2.19283 -0.999515 -0.691817 50.1908 +2.19911 -0.999916 -0.0638029 49.9936 +2.2054 -0.996369 0.564463 49.7969 +2.21168 -0.988889 1.1905 49.6037 +2.21796 -0.977505 1.81184 49.4172 +2.22425 -0.962262 2.42602 49.2402 +2.23053 -0.94322 3.03063 49.0755 +2.23681 -0.920454 3.62327 48.9258 +2.2431 -0.894055 4.20161 48.7935 +2.24938 -0.864126 4.76336 48.6805 +2.25566 -0.830785 5.30631 48.5886 +2.26195 -0.794165 5.8283 48.5195 +2.26823 -0.754409 6.32729 48.474 +2.27451 -0.711675 6.8013 48.453 +2.2808 -0.666132 7.24846 48.4567 +2.28708 -0.617959 7.667 48.4851 +2.29336 -0.567346 8.05528 48.5378 +2.29965 -0.514493 8.41175 48.614 +2.30593 -0.45961 8.73502 48.7123 +2.31221 -0.402911 9.0238 48.8314 +2.3185 -0.344623 9.27696 48.9692 +2.32478 -0.284973 9.49349 49.1237 +2.33106 -0.224199 9.67254 49.2923 +2.33734 -0.162539 9.81341 49.4725 +2.34363 -0.100238 9.91554 49.6613 +2.34991 -0.0375413 9.97852 49.8559 +2.35619 0.0253038 10.0021 50.0531 +2.36248 0.088049 9.98621 50.2498 +2.36876 0.150447 9.93089 50.443 +2.37504 0.21225 9.83636 50.6295 +2.38133 0.273216 9.703 50.8064 +2.38761 0.333103 9.53133 50.971 +2.39389 0.391675 9.32204 51.1206 +2.40018 0.448701 9.07594 51.253 +2.40646 0.503955 8.79401 51.3659 +2.41274 0.55722 8.47737 51.4576 +2.41903 0.608285 8.12726 51.5267 +2.42531 0.656949 7.74506 51.5721 +2.43159 0.703019 7.33228 51.593 +2.43788 0.746314 6.89056 51.5892 +2.44416 0.786662 6.42164 51.5606 +2.45044 0.823905 5.92737 51.5078 +2.45673 0.857895 5.40969 51.4316 +2.46301 0.888498 4.87066 51.3331 +2.46929 0.915594 4.3124 51.214 +2.47558 0.939075 3.73712 51.0761 +2.48186 0.958849 3.14708 50.9216 +2.48814 0.974837 2.54462 50.7529 +2.49442 0.986977 1.93211 50.5727 +2.50071 0.99522 1.31197 50.3838 +2.50699 0.999534 0.686658 50.1892 +2.51327 0.999903 0.0586318 49.992 +2.51956 0.996324 -0.569626 49.7953 +2.52584 0.988811 -1.19563 49.6022 +2.53212 0.977395 -1.81692 49.4157 +2.53841 0.962121 -2.43104 49.2388 +2.54469 0.943048 -3.03556 49.0742 +2.55097 0.920252 -3.62809 48.9247 +2.55726 0.893823 -4.2063 48.7924 +2.56354 0.863865 -4.76791 48.6796 +2.56982 0.830497 -5.31069 48.588 +2.57611 0.79385 -5.83251 48.519 +2.58239 0.75407 -6.3313 48.4737 +2.58867 0.711312 -6.80509 48.4529 +2.59496 0.665746 -7.25202 48.4568 +2.60124 0.617552 -7.67032 48.4855 +2.60752 0.56692 -8.05834 48.5384 +2.61381 0.51405 -8.41455 48.6147 +2.62009 0.45915 -8.73754 48.7132 +2.62637 0.402438 -9.02603 48.8324 +2.63265 0.344137 -9.27889 48.9704 +2.63894 0.284477 -9.49512 49.125 +2.64522 0.223695 -9.67386 49.2937 +2.6515 0.162029 -9.81441 49.474 +2.65779 0.0997236 -9.91622 49.6629 +2.66407 0.0370245 -9.97888 49.8575 +2.67035 -0.0258208 -10.0021 50.0547 +2.67664 -0.0885641 -9.98591 50.2514 +2.68292 -0.150958 -9.93027 50.4445 +2.6892 -0.212756 -9.83542 50.631 +2.69549 -0.273713 -9.70174 50.8078 +2.70177 -0.333591 -9.52976 50.9723 +2.70805 -0.392151 -9.32016 51.1218 +2.71434 -0.449163 -9.07376 51.254 +2.72062 -0.504402 -8.79155 51.3667 +2.7269 -0.55765 -8.47462 51.4583 +2.73319 -0.608696 -8.12424 51.5272 +2.73947 -0.657339 -7.74178 51.5723 +2.74575 -0.703387 -7.32877 51.5931 +2.75204 -0.746658 -6.88682 51.589 +2.75832 -0.786981 -6.41768 51.5603 +2.7646 -0.824198 -5.9232 51.5073 +2.77088 -0.858161 -5.40534 51.4309 +2.77717 -0.888736 -4.86614 51.3322 +2.78345 -0.915802 -4.30774 51.213 +2.78973 -0.939253 -3.73232 51.0749 +2.79602 -0.958996 -3.14217 50.9202 +2.8023 -0.974953 -2.53962 50.7514 +2.80858 -0.98706 -1.92703 50.5711 +2.81487 -0.995272 -1.30685 50.3822 +2.82115 -0.999554 -0.681499 50.1876 +2.82743 -0.999889 -0.0534606 49.9904 +2.83372 -0.996278 0.574788 49.7937 +2.84 -0.988733 1.20077 49.6006 +2.84628 -0.977285 1.82201 49.4142 +2.85257 -0.961979 2.43605 49.2374 +2.85885 -0.942875 3.04048 49.073 +2.86513 -0.920049 3.63291 48.9235 +2.87142 -0.893591 4.21099 48.7914 +2.8777 -0.863604 4.77245 48.6788 +2.88398 -0.830209 5.31507 48.5873 +2.89027 -0.793536 5.83671 48.5185 +2.89655 -0.75373 6.3353 48.4734 +2.90283 -0.710948 6.80888 48.4528 +2.90911 -0.66536 7.25558 48.457 +2.9154 -0.617145 7.67364 48.4858 +2.92168 -0.566494 8.06141 48.5389 +2.92796 -0.513606 8.41734 48.6154 +2.93425 -0.458691 8.74005 48.7141 +2.94053 -0.401965 9.02826 48.8335 +2.94681 -0.343651 9.28082 48.9716 +2.9531 -0.283982 9.49674 49.1263 +2.95938 -0.223191 9.67517 49.2952 +2.96566 -0.161519 9.81541 49.4755 +2.97195 -0.0992091 9.91689 49.6645 +2.97823 -0.0365077 9.97923 49.8591 +2.98451 0.0263377 10.0022 50.0563 +2.9908 0.0890792 9.98562 50.253 +2.99708 0.151469 9.92965 50.4461 +3.00336 0.213261 9.83448 50.6325 +3.00965 0.274211 9.70048 50.8092 +3.01593 0.334078 9.52819 50.9736 +3.02221 0.392627 9.31828 51.123 +3.0285 0.449625 9.07159 51.255 +3.03478 0.504849 8.78908 51.3676 +3.04106 0.558079 8.47187 51.4589 +3.04734 0.609106 8.12122 51.5276 +3.05363 0.657729 7.73851 51.5726 +3.05991 0.703754 7.32525 51.5931 +3.06619 0.747002 6.88306 51.5889 +3.07248 0.7873 6.41371 51.5599 +3.07876 0.824491 5.91903 51.5067 +3.08504 0.858426 5.40099 51.4301 +3.09133 0.888973 4.86163 51.3313 +3.09761 0.91601 4.30307 51.2119 +3.10389 0.93943 3.72752 51.0737 +3.11018 0.959142 3.13726 50.9189 +3.11646 0.975068 2.53461 50.75 +3.12274 0.987144 1.92196 50.5696 +3.12903 0.995323 1.30172 50.3806 +3.13531 0.999572 0.676339 50.186 +3.14159 0.999876 0.0482895 49.9888 +3.14788 0.996232 -0.579951 49.7921 +3.15416 0.988655 -1.2059 49.599 +3.16044 0.977175 -1.82709 49.4127 +3.16673 0.961837 -2.44107 49.236 +3.17301 0.942702 -3.04541 49.0717 +3.17929 0.919846 -3.63773 48.9224 +3.18557 0.893358 -4.21568 48.7904 +3.19186 0.863343 -4.777 48.6779 +3.19814 0.82992 -5.31945 48.5867 +3.20442 0.793221 -5.84091 48.518 +3.21071 0.75339 -6.3393 48.4732 +3.21699 0.710584 -6.81267 48.4527 +3.22327 0.664974 -7.25914 48.4571 +3.22956 0.616738 -7.67696 48.4861 +3.23584 0.566068 -8.06447 48.5394 +3.24212 0.513162 -8.42014 48.6161 +3.24841 0.458231 -8.74257 48.715 +3.25469 0.401491 -9.03048 48.8346 +3.26097 0.343166 -9.28275 48.9728 +3.26726 0.283486 -9.49836 49.1277 +3.27354 0.222687 -9.67648 49.2966 +3.27982 0.161008 -9.8164 49.477 +3.28611 0.0986945 -9.91757 49.6661 +3.29239 0.0359909 -9.97958 49.8607 +3.29867 -0.0268547 -10.0022 50.058 +3.30496 -0.0895943 -9.98532 50.2546 +3.31124 -0.15198 -9.92902 50.4477 +3.31752 -0.213766 -9.83353 50.634 +3.32381 -0.274708 -9.69922 50.8106 +3.33009 -0.334566 -9.52661 50.9749 +3.33637 -0.393102 -9.3164 51.1241 +3.34265 -0.450087 -9.06941 51.256 +3.34894 -0.505295 -8.78661 51.3684 +3.35522 -0.558508 -8.46912 51.4596 +3.3615 -0.609516 -8.1182 51.5281 +3.36779 -0.658118 -7.73523 51.5729 +3.37407 -0.704122 -7.32172 51.5932 +3.38035 -0.747346 -6.87931 51.5887 +3.38664 -0.787619 -6.40974 51.5596 +3.39292 -0.824784 -5.91486 51.5062 +3.3992 -0.858692 -5.39664 51.4294 +3.40549 -0.88921 -4.85711 51.3304 +3.41177 -0.916217 -4.2984 51.2108 +3.41805 -0.939608 -3.72272 51.0725 +3.42434 -0.959289 -3.13235 50.9176 +3.43062 -0.975183 -2.52961 50.7486 +3.4369 -0.987227 -1.91688 50.5681 +3.44319 -0.995374 -1.29659 50.379 +3.44947 -0.999591 -0.67118 50.1844 +3.45575 -0.999862 -0.0431184 49.9871 +3.46204 -0.996186 0.585113 49.7905 +3.46832 -0.988576 1.21104 49.5975 +3.4746 -0.977065 1.83218 49.4112 +3.48088 -0.961695 2.44608 49.2346 +3.48717 -0.94253 3.05033 49.0704 +3.49345 -0.919643 3.64254 48.9212 +3.49973 -0.893125 4.22037 48.7894 +3.50602 -0.863082 4.78154 48.6771 +3.5123 -0.829631 5.32383 48.586 +3.51858 -0.792906 5.8451 48.5176 +3.52487 -0.753049 6.3433 48.4729 +3.53115 -0.71022 6.81645 48.4527 +3.53743 -0.664588 7.2627 48.4572 +3.54372 -0.616331 7.68027 48.4865 +3.55 -0.565641 8.06752 48.54 +3.55628 -0.512718 8.42293 48.6169 +3.56257 -0.457771 8.74508 48.7159 +3.56885 -0.401017 9.0327 48.8356 +3.57513 -0.34268 9.28467 48.974 +3.58142 -0.28299 9.49998 49.129 +3.5877 -0.222182 9.67779 49.2981 +3.59398 -0.160498 9.81739 49.4786 +3.60027 -0.0981799 9.91824 49.6677 +3.60655 -0.0354742 9.97992 49.8624 +3.61283 0.0273716 10.0022 50.0596 +3.61911 0.0901093 9.98502 50.2562 +3.6254 0.152491 9.9284 50.4492 +3.63168 0.214271 9.83258 50.6355 +3.63796 0.275205 9.69795 50.8121 +3.64425 0.335053 9.52504 50.9762 +3.65053 0.393578 9.31452 51.1253 +3.65681 0.450549 9.06723 51.257 +3.6631 0.505741 8.78414 51.3692 +3.66938 0.558937 8.46637 51.4602 +3.67566 0.609926 8.11518 51.5286 +3.68195 0.658507 7.73195 51.5731 +3.68823 0.704489 7.3182 51.5933 +3.69451 0.747689 6.87556 51.5886 +3.7008 0.787938 6.40577 51.5593 +3.70708 0.825076 5.91069 51.5057 +3.71336 0.858957 5.39228 51.4287 +3.71965 0.889446 4.85258 51.3295 +3.72593 0.916425 4.29373 51.2098 +3.73221 0.939785 3.71792 51.0713 +3.7385 0.959435 3.12744 50.9162 +3.74478 0.975298 2.52461 50.7471 +3.75106 0.98731 1.91181 50.5666 +3.75734 0.995425 1.29146 50.3775 +3.76363 0.999609 0.66602 50.1827 +3.76991 0.999848 0.0379472 49.9855 +3.77619 0.996139 -0.590276 49.7889 +3.78248 0.988498 -1.21617 49.5959 +3.78876 0.976954 -1.83726 49.4097 +3.79504 0.961553 -2.4511 49.2332 +3.80133 0.942356 -3.05526 49.0691 +3.80761 0.919439 -3.64736 48.92 +3.81389 0.892892 -4.22506 48.7884 +3.82018 0.862821 -4.78608 48.6763 +3.82646 0.829342 -5.32821 48.5853 +3.83274 0.79259 -5.8493 48.5171 +3.83903 0.752709 -6.3473 48.4726 +3.84531 0.709856 -6.82024 48.4526 +3.85159 0.664201 -7.26625 48.4574 +3.85788 0.615924 -7.68358 48.4868 +3.86416 0.565215 -8.07058 48.5405 +3.87044 0.512274 -8.42572 48.6176 +3.87673 0.457312 -8.74759 48.7168 +3.88301 0.400543 -9.03492 48.8367 +3.88929 0.342194 -9.28659 48.9752 +3.89557 0.282494 -9.5016 49.1303 +3.90186 0.221678 -9.6791 49.2995 +3.90814 0.159988 -9.81838 49.4801 +3.91442 0.0976652 -9.9189 49.6692 +3.92071 0.0349574 -9.98027 49.864 +3.92699 -0.0278885 -10.0022 50.0612 +3.93327 -0.0906243 -9.98471 50.2579 +3.93956 -0.153002 -9.92777 50.4508 +3.94584 -0.214776 -9.83163 50.637 +3.95212 -0.275702 -9.69669 50.8135 +3.95841 -0.33554 -9.52346 50.9775 +3.96469 -0.394053 -9.31263 51.1264 +3.97097 -0.45101 -9.06504 51.258 +3.97726 -0.506187 -8.78166 51.3701 +3.98354 -0.559366 -8.46362 51.4609 +3.98982 -0.610336 -8.11216 51.529 +3.99611 -0.658896 -7.72867 51.5734 +4.00239 -0.704856 -7.31467 51.5933 +4.00867 -0.748033 -6.8718 51.5885 +4.01496 -0.788256 -6.4018 51.5589 +4.02124 -0.825368 -5.90652 51.5051 +4.02752 -0.859222 -5.38793 51.428 +4.0338 -0.889683 -4.84806 51.3286 +4.04009 -0.916632 -4.28906 51.2087 +4.04637 -0.939962 -3.71312 51.0701 +4.05265 -0.959581 -3.12253 50.9149 +4.05894 -0.975413 -2.5196 50.7457 +4.06522 -0.987393 -1.90673 50.5651 +4.0715 -0.995475 -1.28634 50.3759 +4.07779 -0.999627 -0.660861 50.1811 +4.08407 -0.999833 -0.0327761 49.9839 +4.09035 -0.996092 0.595438 49.7873 +4.09664 -0.988419 1.2213 49.5943 +4.10292 -0.976843 1.84234 49.4082 +4.1092 -0.961411 2.45611 49.2318 +4.11549 -0.942183 3.06018 49.0678 +4.12177 -0.919236 3.65217 48.9189 +4.12805 -0.892659 4.22975 48.7874 +4.13434 -0.862559 4.79062 48.6754 +4.14062 -0.829053 5.33258 48.5847 +4.1469 -0.792275 5.85349 48.5166 +4.15319 -0.752368 6.35129 48.4724 +4.15947 -0.709492 6.82402 48.4526 +4.16575 -0.663814 7.26981 48.4575 +4.17204 -0.615516 7.68689 48.4872 +4.17832 -0.564788 8.07363 48.5411 +4.1846 -0.51183 8.4285 48.6183 +4.19088 -0.456852 8.75009 48.7177 +4.19717 -0.40007 9.03714 48.8377 +4.20345 -0.341708 9.28851 48.9765 +4.20973 -0.281998 9.50321 49.1317 +4.21602 -0.221174 9.6804 49.301 +4.2223 -0.159477 9.81937 49.4816 +4.22858 -0.0971506 9.91957 49.6708 +4.23487 -0.0344406 9.98061 49.8656 +4.24115 0.0284054 10.0022 50.0628 +4.24743 0.0911393 9.9844 50.2595 +4.25372 0.153513 9.92714 50.4523 +4.26 0.215281 9.83068 50.6385 +4.26628 0.276199 9.69542 50.8149 +4.27257 0.336027 9.52188 50.9788 +4.27885 0.394528 9.31074 51.1276 +4.28513 0.451472 9.06285 51.259 +4.29142 0.506633 8.77919 51.3709 +4.2977 0.559794 8.46086 51.4615 +4.30398 0.610745 8.10913 51.5295 +4.31027 0.659285 7.72539 51.5737 +4.31655 0.705223 7.31115 51.5934 +4.32283 0.748376 6.86804 51.5883 +4.32911 0.788575 6.39782 51.5586 +4.3354 0.82566 5.90235 51.5046 +4.34168 0.859486 5.38357 51.4272 +4.34796 0.889919 4.84354 51.3277 +4.35425 0.916839 4.28439 51.2076 +4.36053 0.940139 3.70832 51.0688 +4.36681 0.959727 3.11761 50.9136 +4.3731 0.975527 2.5146 50.7442 +4.37938 0.987475 1.90166 50.5635 +4.38566 0.995525 1.28121 50.3743 +4.39195 0.999645 0.655701 50.1795 +4.39823 0.999819 0.0276049 49.9823 +4.40451 0.996045 -0.6006 49.7856 +4.4108 0.988339 -1.22643 49.5928 +4.41708 0.976731 -1.84743 49.4067 +4.42336 0.961268 -2.46112 49.2303 +4.42965 0.942009 -3.06511 49.0665 +4.43593 0.919032 -3.65699 48.9177 +4.44221 0.892426 -4.23443 48.7864 +4.4485 0.862297 -4.79516 48.6746 +4.45478 0.828764 -5.33696 48.584 +4.46106 0.791959 -5.85768 48.5162 +4.46734 0.752028 -6.35529 48.4721 +4.47363 0.709127 -6.8278 48.4525 +4.47991 0.663427 -7.27336 48.4577 +4.48619 0.615108 -7.6902 48.4875 +4.49248 0.564361 -8.07669 48.5416 +4.49876 0.511386 -8.43128 48.619 +4.50504 0.456392 -8.7526 48.7186 +4.51133 0.399596 -9.03936 48.8388 +4.51761 0.341222 -9.29043 48.9777 +4.52389 0.281502 -9.50483 49.133 +4.53018 0.22067 -9.6817 49.3024 +4.53646 0.158967 -9.82035 49.4831 +4.54274 0.0966359 -9.92023 49.6724 +4.54903 0.0339237 -9.98095 49.8672 +4.55531 -0.0289223 -10.0023 50.0645 +4.56159 -0.0916542 -9.98409 50.2611 +4.56788 -0.154024 -9.9265 50.4539 +4.57416 -0.215786 -9.82973 50.64 +4.58044 -0.276696 -9.69414 50.8163 +4.58673 -0.336514 -9.52029 50.9801 +4.59301 -0.395003 -9.30885 51.1288 +4.59929 -0.451933 -9.06067 51.26 +4.60557 -0.507079 -8.77671 51.3717 +4.61186 -0.560223 -8.4581 51.4622 +4.61814 -0.611155 -8.1061 51.53 +4.62442 -0.659674 -7.7221 51.5739 +4.63071 -0.705589 -7.30762 51.5934 +4.63699 -0.748719 -6.86428 51.5882 +4.64327 -0.788893 -6.39385 51.5582 +4.64956 -0.825952 -5.89817 51.504 +4.65584 -0.859751 -5.37921 51.4265 +4.66212 -0.890155 -4.83901 51.3268 +4.66841 -0.917045 -4.27971 51.2066 +4.67469 -0.940315 -3.70352 51.0676 +4.68097 -0.959873 -3.1127 50.9122 +4.68726 -0.975641 -2.50959 50.7428 +4.69354 -0.987558 -1.89658 50.562 +4.69982 -0.995575 -1.27608 50.3727 +4.70611 -0.999663 -0.65054 50.1779 +4.71239 -0.999804 -0.0224337 49.9806 +4.71867 -0.995998 0.605762 49.784 +4.72496 -0.98826 1.23157 49.5912 +4.73124 -0.97662 1.85251 49.4052 +4.73752 -0.961125 2.46614 49.2289 +4.7438 -0.941835 3.07003 49.0652 +4.75009 -0.918827 3.6618 48.9166 +4.75637 -0.892192 4.23912 48.7854 +4.76265 -0.862035 4.7997 48.6737 +4.76894 -0.828474 5.34133 48.5834 +4.77522 -0.791643 5.86188 48.5157 +4.7815 -0.751686 6.35928 48.4718 +4.78779 -0.708762 6.83158 48.4524 +4.79407 -0.66304 7.27691 48.4578 +4.80035 -0.614701 7.69351 48.4879 +4.80664 -0.563934 8.07973 48.5421 +4.81292 -0.510941 8.43407 48.6198 +4.8192 -0.455931 8.7551 48.7195 +4.82549 -0.399121 9.04157 48.8399 +4.83177 -0.340736 9.29234 48.9789 +4.83805 -0.281005 9.50644 49.1344 +4.84434 -0.220165 9.683 49.3038 +4.85062 -0.158456 9.82133 49.4847 +4.8569 -0.0961212 9.92089 49.674 +4.86319 -0.0334069 9.98129 49.8688 +4.86947 0.0294392 10.0023 50.0661 +4.87575 0.0921692 9.98378 50.2627 +4.88203 0.154535 9.92587 50.4555 +4.88832 0.216291 9.82877 50.6414 +4.8946 0.277193 9.69287 50.8177 +4.90088 0.337001 9.5187 50.9813 +4.90717 0.395478 9.30696 51.1299 +4.91345 0.452395 9.05847 51.261 +4.91973 0.507525 8.77423 51.3726 +4.92602 0.560651 8.45534 51.4629 +4.9323 0.611564 8.10307 51.5304 +4.93858 0.660063 7.71881 51.5742 +4.94487 0.705956 7.30408 51.5935 +4.95115 0.749062 6.86052 51.588 +4.95743 0.78921 6.38987 51.5579 +4.96372 0.826244 5.89399 51.5035 +4.97 0.860015 5.37485 51.4258 +4.97628 0.890391 4.83449 51.3259 +4.98257 0.917252 4.27504 51.2055 +4.98885 0.940491 3.69871 51.0664 +4.99513 0.960018 3.10778 50.9109 +5.00142 0.975755 2.50459 50.7413 +5.0077 0.987639 1.8915 50.5605 +5.01398 0.995625 1.27095 50.3711 +5.02027 0.99968 0.64538 50.1763 +5.02655 0.999789 0.0172625 49.979 +5.03283 0.99595 -0.610923 49.7824 +5.03911 0.98818 -1.2367 49.5897 +5.0454 0.976508 -1.85759 49.4037 +5.05168 0.960981 -2.47115 49.2275 +5.05796 0.941661 -3.07495 49.0639 +5.06425 0.918623 -3.66661 48.9154 +5.07053 0.891958 -4.2438 48.7844 +5.07681 0.861772 -4.80423 48.6729 +5.0831 0.828184 -5.3457 48.5827 +5.08938 0.791327 -5.86606 48.5153 +5.09566 0.751345 -6.36327 48.4716 +5.10195 0.708397 -6.83535 48.4524 +5.10823 0.662653 -7.28045 48.4579 +5.11451 0.614292 -7.69681 48.4882 +5.1208 0.563507 -8.08278 48.5427 +5.12708 0.510497 -8.43684 48.6205 +5.13336 0.455471 -8.7576 48.7205 +5.13965 0.398647 -9.04378 48.841 +5.14593 0.34025 -9.29426 48.9801 +5.15221 0.280509 -9.50804 49.1357 +5.1585 0.219661 -9.68429 49.3053 +5.16478 0.157945 -9.82231 49.4862 +5.17106 0.0956064 -9.92155 49.6756 +5.17734 0.0328901 -9.98162 49.8704 +5.18363 -0.0299561 -10.0023 50.0677 +5.18991 -0.0926841 -9.98346 50.2643 +5.19619 -0.155046 -9.92523 50.457 +5.20248 -0.216796 -9.82781 50.6429 +5.20876 -0.27769 -9.69159 50.8191 +5.21504 -0.337488 -9.51711 50.9826 +5.22133 -0.395953 -9.30506 51.1311 +5.22761 -0.452856 -9.05628 51.262 +5.23389 -0.50797 -8.77174 51.3734 +5.24018 -0.561079 -8.45257 51.4635 +5.24646 -0.611973 -8.10004 51.5309 +5.25274 -0.660451 -7.71552 51.5745 +5.25903 -0.706322 -7.30055 51.5936 +5.26531 -0.749404 -6.85676 51.5879 +5.27159 -0.789528 -6.38589 51.5575 +5.27788 -0.826535 -5.88982 51.503 +5.28416 -0.860279 -5.37049 51.425 +5.29044 -0.890626 -4.82996 51.325 +5.29673 -0.917458 -4.27036 51.2044 +5.30301 -0.940667 -3.69391 51.0652 +5.30929 -0.960163 -3.10287 50.9095 +5.31557 -0.975868 -2.49958 50.7399 +5.32186 -0.987721 -1.88642 50.5589 +5.32814 -0.995674 -1.26582 50.3695 +5.33442 -0.999697 -0.64022 50.1747 +5.34071 -0.999773 -0.0120913 49.9774 +5.34699 -0.995902 0.616085 49.7808 +5.35327 -0.988099 1.24183 49.5881 +5.35956 -0.976396 1.86267 49.4022 +5.36584 -0.960838 2.47616 49.2261 +5.37212 -0.941486 3.07987 49.0626 +5.37841 -0.918418 3.67142 48.9143 +5.38469 -0.891724 4.24848 48.7834 +5.39097 -0.86151 4.80877 48.6721 +5.39726 -0.827894 5.35007 48.5821 +5.40354 -0.79101 5.87025 48.5148 +5.40982 -0.751004 6.36726 48.4713 +5.41611 -0.708032 6.83913 48.4523 +5.42239 -0.662266 7.284 48.4581 +5.42867 -0.613884 7.70011 48.4886 +5.43496 -0.56308 8.08583 48.5432 +5.44124 -0.510052 8.43962 48.6212 +5.44752 -0.455011 8.7601 48.7214 +5.4538 -0.398173 9.04599 48.842 +5.46009 -0.339763 9.29617 48.9813 +5.46637 -0.280013 9.50965 49.137 +5.47265 -0.219156 9.68558 49.3067 +5.47894 -0.157435 9.82328 49.4877 +5.48522 -0.0950917 9.9222 49.6772 +5.4915 -0.0323732 9.98195 49.8721 +5.49779 0.030473 10.0023 50.0693 +5.50407 0.093199 9.98314 50.2659 +5.51035 0.155557 9.92459 50.4586 +5.51664 0.217301 9.82685 50.6444 +5.52292 0.278187 9.69031 50.8205 +5.5292 0.337975 9.51552 50.9839 +5.53549 0.396428 9.30317 51.1322 +5.54177 0.453317 9.05408 51.263 +5.54805 0.508416 8.76926 51.3742 +5.55434 0.561507 8.44981 51.4642 +5.56062 0.612382 8.097 51.5313 +5.5669 0.66084 7.71223 51.5747 +5.57319 0.706688 7.29701 51.5936 +5.57947 0.749747 6.85299 51.5877 +5.58575 0.789845 6.38191 51.5572 +5.59203 0.826826 5.88563 51.5024 +5.59832 0.860542 5.36612 51.4243 +5.6046 0.890861 4.82543 51.3241 +5.61088 0.917663 4.26569 51.2033 +5.61717 0.940843 3.6891 51.064 +5.62345 0.960308 3.09795 50.9082 +5.62973 0.975982 2.49457 50.7385 +5.63602 0.987802 1.88135 50.5574 +5.6423 0.995724 1.26069 50.3679 +5.64858 0.999714 0.635059 50.173 +5.65487 0.999757 0.00692014 49.9758 +5.66115 0.995854 -0.621246 49.7792 +5.66743 0.988019 -1.24696 49.5865 +5.67372 0.976284 -1.86775 49.4007 +5.68 0.960694 -2.48117 49.2247 +5.68628 0.941312 -3.08479 49.0613 +5.69257 0.918213 -3.67623 48.9131 +5.69885 0.89149 -4.25316 48.7824 +5.70513 0.861247 -4.8133 48.6713 +5.71142 0.827604 -5.35444 48.5814 +5.7177 0.790694 -5.87444 48.5144 +5.72398 0.750662 -6.37125 48.4711 +5.73027 0.707667 -6.8429 48.4523 +5.73655 0.661878 -7.28754 48.4582 +5.74283 0.613476 -7.70341 48.4889 +5.74911 0.562652 -8.08887 48.5438 +5.7554 0.509607 -8.44239 48.622 +5.76168 0.45455 -8.76259 48.7223 +5.76796 0.397698 -9.04819 48.8431 +5.77425 0.339277 -9.29807 48.9825 +5.78053 0.279516 -9.51125 49.1384 +5.78681 0.218652 -9.68687 49.3082 +5.7931 0.156924 -9.82426 49.4893 +5.79938 0.0945769 -9.92285 49.6788 +5.80566 0.0318564 -9.98228 49.8737 +5.81195 -0.0309899 -10.0023 50.071 +5.81823 -0.0937138 -9.98282 50.2675 +5.82451 -0.156068 -9.92394 50.4602 +5.8308 -0.217806 -9.82588 50.6459 +5.83708 -0.278684 -9.68903 50.8219 +5.84336 -0.338461 -9.51393 50.9852 +5.84965 -0.396903 -9.30127 51.1334 +5.85593 -0.453778 -9.05188 51.264 +5.86221 -0.508861 -8.76677 51.3751 +5.8685 -0.561935 -8.44704 51.4648 +5.87478 -0.612791 -8.09397 51.5318 +5.88106 -0.661228 -7.70894 51.575 +5.88734 -0.707054 -7.29348 51.5937 +5.89363 -0.750089 -6.84922 51.5876 +5.89991 -0.790163 -6.37793 51.5568 +5.90619 -0.827117 -5.88145 51.5019 +5.91248 -0.860806 -5.36176 51.4236 +5.91876 -0.891096 -4.8209 51.3232 +5.92504 -0.917869 -4.26101 51.2023 +5.93133 -0.941018 -3.68429 51.0628 +5.93761 -0.960452 -3.09303 50.9069 +5.94389 -0.976095 -2.48956 50.737 +5.95018 -0.987884 -1.87627 50.5559 +5.95646 -0.995773 -1.25556 50.3664 +5.96274 -0.99973 -0.629898 50.1714 +5.96903 -0.999741 -0.00174894 49.9741 +5.97531 -0.995805 0.626407 49.7776 +5.98159 -0.987938 1.25209 49.585 +5.98788 -0.976171 1.87283 49.3992 +5.99416 -0.96055 2.48618 49.2233 +6.00044 -0.941137 3.08971 49.0601 +6.00673 -0.918008 3.68104 48.912 +6.01301 -0.891255 4.25784 48.7814 +6.01929 -0.860984 4.81783 48.6704 +6.02557 -0.827313 5.35881 48.5808 +6.03186 -0.790377 5.87862 48.5139 +6.03814 -0.75032 6.37523 48.4708 +6.04442 -0.707301 6.84667 48.4522 +6.05071 -0.66149 7.29108 48.4584 +6.05699 -0.613067 7.70671 48.4893 +6.06327 -0.562224 8.09191 48.5443 +6.06956 -0.509162 8.44517 48.6227 +6.07584 -0.454089 8.76508 48.7232 +6.08212 -0.397224 9.0504 48.8442 +6.08841 -0.33879 9.29998 48.9837 +6.09469 -0.279019 9.51285 49.1397 +6.10097 -0.218147 9.68816 49.3096 +6.10726 -0.156413 9.82523 49.4908 +6.11354 -0.0940621 9.9235 49.6803 +6.11982 -0.0313395 9.9826 49.8753 +6.12611 0.0315068 10.0023 50.0726 +6.13239 0.0942287 9.9825 50.2691 +6.13867 0.156579 9.92329 50.4617 +6.14496 0.21831 9.82491 50.6474 +6.15124 0.27918 9.68774 50.8233 +6.15752 0.338948 9.51233 50.9865 +6.1638 0.397378 9.29936 51.1345 +6.17009 0.454238 9.04968 51.265 +6.17637 0.509306 8.76428 51.3759 +6.18265 0.562363 8.44427 51.4654 +6.18894 0.6132 8.09093 51.5322 +6.19522 0.661616 7.70564 51.5752 +6.2015 0.70742 7.28994 51.5937 +6.20779 0.750431 6.84545 51.5874 +6.21407 0.79048 6.37394 51.5565 +6.22035 0.827408 5.87727 51.5013 +6.22664 0.861069 5.35739 51.4228 +6.23292 0.891331 4.81637 51.3223 +6.2392 0.918074 4.25633 51.2012 +6.24549 0.941193 3.67949 51.0616 +6.25177 0.960597 3.08812 50.9055 +6.25805 0.976207 2.48456 50.7356 +6.26434 0.987964 1.87119 50.5544 +6.27062 0.995821 1.25043 50.3648 +6.2769 0.999746 0.624737 50.1698 +6.28319 0.999725 -0.00342225 49.9725 diff --git a/szamszim/f1/plot1.gnuplot b/szamszim/f1/plot1.gnuplot new file mode 100755 index 0000000..c90304e --- /dev/null +++ b/szamszim/f1/plot1.gnuplot @@ -0,0 +1,19 @@ +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 + + + diff --git a/szamszim/f1/sho.cpp b/szamszim/f1/sho.cpp new file mode 100644 index 0000000..b9664b7 --- /dev/null +++ b/szamszim/f1/sho.cpp @@ -0,0 +1,77 @@ +#include +#include +#include + +#include + +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 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 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; +} diff --git a/szamszim/f1/sho.o b/szamszim/f1/sho.o new file mode 100644 index 0000000000000000000000000000000000000000..9f61eb905ac57e19b07c4a118121c57fcae15c89 GIT binary patch literal 23088 zcmd5^eQ;dWb${|=Fa(fHfS8XYYvF)Ruqycr8%&Uq*IpTa2un7$vHiSSt*kB5?y|cU z5?^*>)1V4xL;SXj)C~cA%Gfgrk1?u2T3`r;B7Lt%aiki@=6Vuowv{jQ#)bqRV zocGSY4{2rcN6+}_-uJub>z@1Z-aB_;Z({ALX*D&LP>uB|%excQvI@6)^=4CTw(hq6 z)bP%ud&qo}%Nw|yPh|tw|Afm%Dx0{zkjh0|U(DqaDsSZaQZAQqc@vj6Q@Nb$w@~>h zuEWJ$g+6mh0=NT+j6lT()s}2bCMS-cDr) z*EdnQnd^6QxrNHDT;E3JU0mNziV}KFh-4uvF_Tr)UtJ*EB;>X7Bm}GB0eG0$c|$l$OmI- zLJN~S>iSS;x-Ncmw=Y9D6m}w_{xbEwB1}%t+LQa=O(usPJmQJXOuRuiK3h9YSWgUd zhfcrxUhVZzmmqsi+bi9*r;nP!N)9JZb|-6dN0SE)dvfUbNqG6v`SeF_l0z>iS({-U z`tjNF95aGs@S&rQ)uFc?RTn68#}n^h(f;yH`g1%nN`I|izCeGH!<$|qe~&=&W$_LD zqB*lAK|Q3&_y4g!`xB!zPWYB3Mr$8kJU%|YpSmc}JGGy^oS>1BX`fV_a`e_&^H+BKZua@jPI*PM+ppoo zn6nR>UV3-?->nlQIow>AeE4YjV-$d~wYV8xH-`GRVx&Gf^aGlYE5=WZyn1YWyv97S ze3G@FKgu?=R-W;WSnXZ(Hu0sJ>3?5iov5Kb6aHamo<*u?Bbu2Q)$8cM(eiY&!x?aR z%Pd+TqwcD}zKS)6B{=>XZvd5k9XD}F9&et33i?9zI@I}DL#CbnJL66d($hz4XVO#m zN9aDW%X)Q0XR(?yf|G~KiO;<^n?1zzzA)zQ1#Pp8QrQ~X>L@A^L$v?FF{}0o2`Hu} zx!?;gPgZdzXq?r%5vJ*l0IL`tqg!v^8u}N}#o?6knTzp(^w6fah7)hmEIDugSWWrL zQ%fu4rDJucmQ>)hwN+S+ow|H@(_1vZc(^gO%s`vHP;cctHP~_G>39DDW9rO~fO_RU zjd0vJU(rt*A~m`nI84s0lewfc^YBKtU}Su;-9Hv97wOO5`}U7bD<|pC-hvnJOhsJ1 zc(1F#0|7ZTzj6T7(lH@Cb&-}I@ePB!N`gJBoR}4N{#oY@-Y}MdHwlL=^E%_MVY4D6 zWYx~`Fdsu#SF~G06U?b$D>)Ra*ki=OikM)3n26d|^>-yg7y|I>mX@31b5~|Mv#DIX zX<@^HhQ|5J=5fukE@d4x_ieOl2I^|Ay70UiBlOO5HQ~5Gm*CHei|hV;TFXV}t|JDz zMcZ=VL>J?aw6}Wg+kq2BwWF=(zP7c-@_Apx?NwqwPB(La5aBsxp4DPwk7wg37>LY< zD(&rF`>1KPtOF>eecE2bVETu5T%xG`q&AYaEXT~TA!$n8;cqHPlS0=MHE-}SH8^KXRP@y+hGQy|x8Lt; zpLUO@<8|g~@$n223;Oszw<90n+UWM=FU4yLQJ4>qU&QTKf=R}FlPBLXL3^A;%G~~p zmd7|g>FGE$fe!d^gxjOW>5WO+$G9E&5w{QILy2QF%S=Qk7vr|8aazjlVRIV(ZQ=GP z{|2}nxfWMH*?)lS#|b*b?O`#7NEBgupJV;!cvySP-BBuAW7f5m0N=%|c@el-T@`Kk z7OJs4oR^ug2wd)^baGbLyWcp1Lut%p$|}^Go0x^nS&(AZ{2)k3%w*SsO)(4EwjjkU zWZ!}mvlawFLShy&bU}()i-RB`F$>wdAjPbuL6DG`g(w9nW}XOP_{~9(keFp=d7v$3 zsidW9jae(yomY)nE#8f)#jLO=G3+9P7^Ilh8UzW6S;#B~DQ2w;f`r5@WE_JOvyg!d zQp_?>jzB}qLUuApF>6y0BqU}bYZ;`NwIv7=60^1iaI@}1zQA`G&T|9&g9gEyM>pUFN&3lHK~L;WtkN(_;25s{OaYy>%z}hiSpqt95sik+Q7$ zgwLMV$#<3|9xoz1%x-Ihp0vTeb!y&kNbAFd&$Q+S^ba#G_N2Ur_IdD&ZNU#Rzv$)d z)_Q_*?BO86zexBcYX3aoG@Ex4;BqhJH<^D1l_2x}MCG>`-^sXnZzB9BbT2-FU+5PR z_X)t1@qIHMqAw%M)vB1YR<@w{MF-PZ-?WpJjdh+~Ber z_}!oUzmg_Wz2)uC=KYWGyurQw*}QWR{xySp`}6&Tfu>Wgr(W&PEe>qu4K8=iSreMR zWpHnQ7X346i>X)p^OdZpnef@xcD}Q$IM>z~esA9v{(B7W?c2{V|34et+qZ@Pw+8q2 zZ5fxV&VxT~u3f=93@)?Sa-U?MGq|^JOFT~++}pQhp0A;WU$6FU(SN|;-oCwxJveRf zZLUA!|7En*hOLWR2%l~3bor^9hE^LKCsq!LY+Hn$LBpRCc8iA3A&ioG;p`;%*0AlHHxl7c+FTzP;SE$X>N4 zmqrTI?o5?3D;FpFotCBDjJM)5I^i|=MOvanS{Ao%Y;Uyt6NwESHYv^B+t&u=?d1hc z*?h^^vwO3h#Z+<7&gQaZw!F+hGG|)SZwZ5Taw9Tin|BV&4s0GUwhdXW5V0Zi_-Z$7{2xG`QB0jV(sPLQY-bY zDczGQ+T~&@TQ0Sx6N!z4Eu=xqrtLh9dM4G^n(nX@?ThU+mF837k_tJ#vvSn7KU%WR}>}u`8MB9j2 zy`-4w%lBt&8u?5ylS^k>X(A>Tl?Mx%gyV&OXj|M61X_4L1ftfmpu+R8<6@PtME|_?<(9&)A}f9qWXybWoahwE~a{&O*EA& zw|1`WY)VATmvBooch=dvGiAGwhq|IK9A-XKcR5&gv9*7rnG$Pghovwr%j8m>y%{^( z-MXbzq$!Xs+5P!!S1T>+#r-%H!|LXFyrHAfURl`C(S*Gf!XevcYa50G)1FLr=V{OA zu(4=5oS8V0pfgoP0AuKTR76BzqQP7=!C+R*8*k(1B}0L9JSRVzY&+G}RivjvzPt5y z$4^=`8)>eX0DB9^X`5i>ERwiw^(Q$i8@oPQbSGd;mL#0R0|DP*?7D zwVwY`xQ^TTl)nJ^ou^n6>W@TH3XZiS;NdY)7`=CiY_j*mahIQ+r4 zx$ys;kDqU?mlXeM#s58pCl&rfr3cv=(LbvATNVDU!jZKRemPHqpDhZ%7yywUze3?_ zs9p3JX8MC=gE&ZDN%983u@4Bvxgz|M&usMJ*YSO`!Z8np+PETmOz}x5& z9CZ-zJLnR;6Ce@rjdTg#2apKs6j?LD9{@-M9KX*NT;`EqzvN4SZ&LhMP#pw3s$ZL_ zP4q7WGZFYPwt{Q@TNK_zbrAgAN|*5829OAPkgXFOF#!SBc6$x#;CE47__hBH3cptI z-=J_=gV3+zujk2j#Xm>s*`aW4H!5z(DvF;vZnna;UG0Aqj%OJN{_8j=3>u6xVh94i zj&sJK!8rFv;EyPL7uBU-f2r`4!nOZ64HAAGHyt0nj&wZreC||w^n6PzT-)tZxVEd~ zT+KYucC~)3hfEp5{)5=oIk)aP8t`t$P)NUk?^d`}f%niArbqTp;YA7l1qz6B9qSHb zFVxHx;FPj2JE=p1y22pFyd zxFOQ-`M~9zPQSMUmvg!vrRQQOV0uL1;Foi;#}r;G%v2s$xQq#!Usw1{#s9v-H zm%D(qLaX)oC?Q-p46>q*l!IrW#CYo{m4LisKZk zp|e!74B8;aF%8AM(*#~ZBT_AN3?OniSP08CY8I0$t2|gje!XFLI&EhL(wRcpCVzV~ zjzNFo%p6B(bAASW4&!vs^WJ}g5o*_}l%gK_#LQ8N>q4+Q=D46zNkB*C`G6=>3AZZG zD@3yDgxloj3!+jvGlf$If4}|16OkDXCO#t$90K~!_UvIwiq7sBB4DDbbgfQ_h6iE7 zjOf`JrD2-zaH2vh+5h98Q&d=p*lEI#EUJ^#XBK{mDkte#9bZ)1ir613Pcpm}IN^DQ z$97hy8WoPwke>AcM}?xQH=>X;#)ty|--aCAa3l5<4A<0;J;I_omBWv-i-unG5n-6p zstdV1W-$sa%Wf$(uiIJY_O(98_TJB87V-pa#D01T@`Vujt|`dNIkt}9)+xwKo=wjm zhN|BY#@0{L>-!kq7n5#|Z=LxpZkd%<|WW05>x1vLcrTmSz`m zDs`M(+Qo`{2%bZi#B&z&-`aVL6%SjCeB4-PLas#94eaG9hY+@V?J@KjqM+CIh{(K?|v=I z`|ZzQM|JKPqJBHM2>Hvalz%Kjeh$mu1fz=l_Mf9R#AKSUVJFLDY&H4g_cLk?i+?A} z%RCbuJ~P$@BpfDxAIs}pF?9Ilr<1Umf4+tpcv&YB`#n@|q>I(M_&m1n6StmVMxASi zKKQDsndSZC|18Vjp_!@fw{vrZ|8e#|$?`o^;w+HWy8iEqu>V_Tl=CiD6R?l}zKZ?@ z*3TS1Tdde=*BT0be!!dxZS+EI--$zlk~o zo&NFrewFqgjFA5c%VU4j7Tz%Y=5%x z-$3$V@!!RU_4p@3?0+-Deu3@lGeE?FrdtI)7h(Tdwy*c^)gGk_y4mI_Sf)>_hj)O zim<*e%l-W?MA(0k?N7G<-yLE9l@R;LzxeI{Q-u9rvwi>jlfVD}7-4^m?d$PJ{>5*9 z8NDZm#s7T%vBqTM|4@YdM_7K99z3f1{U0Rxu<@%~;0PvJzK6;!bou>%Kf?Y_wm;eW zIUQkt3)|Q6+ZJN~6XX#36!WeB0k$8fHcjHcU~vM%`d zM+N!g-_G)W|DR&{x(aE}zQF4v{6E6}>-C5Hjo<%lg#U-xK5Z`* +#include +#include +#include +#include +using namespace std; + +#ifdef __APPLE__ // Mac OS X uses different header +# include +#else // Unix and Windows +# include +# include +# include +#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(); +} + diff --git a/szamszim/f2/pendulum.cpp b/szamszim/f2/pendulum.cpp new file mode 100644 index 0000000..dd94e47 --- /dev/null +++ b/szamszim/f2/pendulum.cpp @@ -0,0 +1,76 @@ +#include +#include +#include + +#include + +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; +} + diff --git a/szamszim/f2/pendulum.hpp b/szamszim/f2/pendulum.hpp new file mode 100644 index 0000000..0dd0d48 --- /dev/null +++ b/szamszim/f2/pendulum.hpp @@ -0,0 +1,9 @@ +#ifndef PENDULUM_H +#define PENDULUM_H +struct State { + double t; + double theta; + double omega; +}; + +#endif diff --git a/szamszim/f2/plot.gnuplot b/szamszim/f2/plot.gnuplot new file mode 100755 index 0000000..c90304e --- /dev/null +++ b/szamszim/f2/plot.gnuplot @@ -0,0 +1,19 @@ +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 + + +