From 0515717e6d5b57feb7d6aafc176e60dbe0b1d6e2 Mon Sep 17 00:00:00 2001 From: Jean-Marc Petit Date: Wed, 25 Feb 2026 14:46:37 +0100 Subject: [PATCH 1/2] Adding licence files and .gitignore. 2026-02-25T13:47 UT. --- .gitignore | 8 ++++++++ eupl1.1.-en_0_0.pdf | Bin 0 -> 33994 bytes eupl1.1.-licence-en_0.pdf | Bin 0 -> 34271 bytes 3 files changed, 8 insertions(+) create mode 100644 .gitignore create mode 100644 eupl1.1.-en_0_0.pdf create mode 100644 eupl1.1.-licence-en_0.pdf diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..e911088 --- /dev/null +++ b/.gitignore @@ -0,0 +1,8 @@ +.idea +__pycache__ +build +ModelUsed.dat +OSSOS-model.out +*.log +*.ipynb +*.so \ No newline at end of file diff --git a/eupl1.1.-en_0_0.pdf b/eupl1.1.-en_0_0.pdf new file mode 100644 index 0000000000000000000000000000000000000000..d69ea0d0460755bd04c0f8a6f023f1f48ca7faf7 GIT binary patch literal 33994 zcma&Nb9^Ar_AMIQ_Qaalwrxyo+tx%A+qP}nwr$(V&3w-}?|1L#ornIbyQ{it?X}mg zuI}1JA}1t5O;5uDOR~E+H3tjD$cj&gZ=-Jx%f&@2WoB)pXK&_iWQb2sD}~R<$j;70 zD~r#_PAiPhz)H_VM=SGPp!8k9!p1}^hOdFo#KMA4&&r_1!vkw%ZTPnt{Quhl3&q6p zKYa)p8Q2&a$?4f!IehmhXk+PQW$l2^L@Q@+WN2pKXk(90&%nb&D`IBpXkjBu6pV~%RW<1G8Sv@932<<1x5=0g)e|2Hv9>Ulpf~`g04mgA%>9!h{R~jNJr0L30RAPfgw)x z$yoUchM)_uqYViN3-!(8znl6uqYjRG_KvRhM#ivEbg)n)B*L;Huu%UO%dk+)f3f}B zgM#?%4F8ooe0qAOZ@8Z~nK7KVeZo`KP~I7}Sy>HaDGmaL<_lhLUW)4D}rK@EQMcnXI6k zy^XDry`!1YU%|-=iaKa8{nfm}x6f$h^h}JT^z@A^9lpK&PyK)0akFu9q?NTc`jfBg5~I`it!U)drsLgI4hyqk*H6y^+yBs`>}Uzr-Z?jo8x0 zEBp|T^&Ug9lyixyFm1B!CycM z-vuJyX#Q4a(E0A2=^sE!O5dmd#L!<(?QIMcjT|*-ziFXWGIDj)`g`ktTp%s@UHKo= z|KUWz=5Ht8Zu}QH|APA0t^WksUw8ZuDNKJY{x>QAWrPqhvv+XBXJh`4##tE|nEwu0 z-(%L_fR)uZS9G-c-xXh?6ch9{`eJIZ!*qp?kAOd>LWSR1UZ4eTdLJCR5Ca}# z{5u5_>Fdoai*RKQOFX6{tbyyOTD=`*M1MoJG;K{`%@a!FG>z7J%o>`Xwkr^;vHZ(m zF-x~&o|+6wnc3}b!D#^nRX;?dl5&2tk@wT?X=5O5&twaSRvXVIQ(I>X*M?`SfK2%BQpC^K4+_u}@XDtK<3Y#<yVK({s%w*eelM-fM|+;g%?jpYjbns{dHw8X~~D%e{Zb`IA=mVYp99$uCHz1 zSvJ$_;oY^V0D_^74Be)SUAy!3!1^c2c_&9?kdBb5i>!`Ld+WyN<;lasMv{*9m&KPb z%>5r_;iI)cz0|B(@+>4pIDs9cLlNHr1 z?e>id=`F8^I2OL+c1{vU$s}RWQwiW`=qK_^@6|T13M+@;>H%ug%=##LbFT|I7Q_Yz zjWyPV=}}gfYwtVrGOQm@z7FDxF_fD*jR990ri<+N)W0^+I#Sink{`v7rUSOjA%fMK zmK{&sA37)-bJ^`_z%TXlLmwxq|5VtN@mqXkr9P-Of?RrGHb4qlgUn zzNPJt4{Z>SYxxU{s$;He3^c zp>b=ufJTxB5KjxbJ`kQm#W(089y95e zg+~fO6>sh#5eq7u1U|m()sM~MBb1SvIaq(GI=F1;VE_%n|lTo5tAqa z%oW%Tx*T{{^WPW1)w0fbJ$}T>&@{LK|6-sNmtL=r(8&1JP!qP$C+(l1dg}kMQn=JZ z;8-Yqlnhq1oZ4)c@D5dwwXYz46U}bHy7&{Ty7LT8MSgSFdQrtO`ePkgzT7SlEyyyD zCfvIf@r{zE$J2tUjbuHV&UPCnnUz3gAf-6C(_ft15dO5{k-087K#*8GYdP$JnPWEM z+rJS<=h7K=V)KU9PV{i3FJ;V$bS<9iXX)^R7h#K`Dq5qRfAUt;#SezJ6^_(ZNA5XH zR|hQ>uU%~CUU%}(z?vvH`djeXNzl&1aAU@^_avmFLKD8>#C9M&ez@3_M-xP z66Pm$Z*)1~eC8|T@g<(h6&%ck#ovo1z={}#rb~MQ-JM&HjT4z(ETM9ltc7qK8tqfb zdNOKxQOxAI>u7eec2?Z>GSSSSj0BSuu`}dKit~hz`i8|7Sj9_)MEqTpjfDVFf_6v6 z*BPn}rcs|r?8(Yzuv7Wt!Qx7dQCN+h6pYau3?Qlc)G}6cG7?4@9R`5{v}4-6Y|Rc! zKDuL+Sr6;#4h;#h!qzO#1p~o6Yiml0g&!6c<-FE>OmSyqUOt$8zxp?9Sv+}te;UwL zJ}~lN<_nW@vg<`RH{iN~nABv*Y~F_oE?{OyhioVsNwM*TJ& z&Wl@WR$yIkVrA7hqlJJHIv;!~=dqcz45!=#mF6@rJ{7Q++b6zdW_^q`anq)jcG!<8 z$wZavTwL?(BrO8_bOrZxH`87eoFz1=DYy%|Ra&s~P!M^@+-`eXsBJ>3Z*aUkYZWP6 zl=UWh(In&2_H4$v+Ucfee1h{Lm#M0AGq`fWaX)6Viqp&tnWqM!y|i=1yzT5wk+-7M zW6u5IkT!y4mTK9dI;nl>ie*iV6T-S4()6XyX4~cw_N6{Ri?z;x-VEI1_~~I5o#QF+ z6=UkOU?)U1%Kbrl*m84oAG>YKiMRFH)!-8({DuO0uv0$k1v*cAe(>T6u1pYQP56&A zu%lKZ4Hz%rsn$Bcv68zPX-vM2_?FC|WHtxqmW$7i)@FXz(RmEzHjQXvd6mdPGre#k z+|y2Nv5+tn8EK}q9OQ_#wdsyp(;5T8P$0Gv8p>nbG@v?huea_z7u|{L<6Nt8BTgXH zcaE+e(7*ypX#kWwB-B%H#1lf~t7`uYcn6ewDc7#+_u)6XhhDmSa8yG0bVPT0onq)- znJ`icQIBUxdPjX6gu@6Ueo8Psx)2y?6<5I= zmo^Ityw7kSO&wEj5g#A-Z(HLYosnCiPyFCCuTUEhq+(C}JpvhgI|6iRt8=3#R=5VA z%_l<(=ynlJ9!K4`D75QYZW}G1&E$@l`sq)MA{$dc1qAR9bEmhCyuy%mpxqvb5yS@e zZ{VI^-0ayRjtCR0xKZrqvM$(d5>b}VZu>QOPP=gx8d-EGz2u~9x`8qEk&tzz?YiOq zv7EZ$3a`JFMM}JTubqaMHR(L4G%54u1EIC?#-{--!e`NFsM8Gtz15^y%yIuwrl5#zaB?|3?9to^?&_S-SkV=jlTz+mHgxLqooNdoBA31)C z|1in{S&>TvqkMNzYW^gcsXilV^@$Hv@yptq%fw$2)B_~IgQj*!D!qC)w*>featqvsEOjE19EFdXv~2?~a=*raY`k>U-|v?_a*zjdB!_OO{CFuNM&9u= z**KnHz;)k73R&1F1W(-3u6!qX)t#EqAU}Bb1pN0K>?a= zU}{LGtAOsNt2Lzc)l7Toz8Jl#MF_BO-6$m{`Mms!>$u0SWnZKtZfg997_Lfvzr7^3 zD>K9)pJKz@!@kcIo)7H`cg!I0%c)x$~`mWB4p*x!5K!L0lFxMEeK=(v})R<#D4JPwaj;AmhJ z95^yYb|QZwuW%b)_ZT1|SK+|fSn3;mP2=Udx*IOFO+G`j?5C)z^znf)6W?{YudBGn z?*tA_WQ%HazI6tB6a`ZFkQTin?H&Kw&E5b@ydMNMK6dlRtxcd7a;iLQiE@ACd0_Fm z(LCW)+RQ*_@Ggntg^`Z)TOWfSDOOR{%-_L%C}g|}0;;_dL}>=*>K z%t?cvZZgfmR(1f?fMHB#$R%ggoY&ukhu1C8D=19ykq1rjjsRvzMZpM-xZ!a>eK(Fb zn}lPvV%*$IqahnhOeI55XhqDeHZ@=&i8pu8c2;Kw!i?%Zru9pIb<7ktsT@x4@wZqs zPv603DhWyfSF~zKrd0s7m{%d1F{g4sMU(a~Th6x<;S8lve!9$jM5LAKK+w>K?3k2J z-yw29wuip8{aWFlW1jQ~J?hzWueN^1yKZGd!EF_b5rX$CSOAMa$(uQ^d2kZ9WpU+# zrB~yUOy@8u4YfUS2AjIhIR|b&eFj$G8`q#LO|`rBI$2;t65JJS*Tqj+OD;yhQ`D$BeXcI)DMVY()ZUOu$qDcd& z2(G7Q)p*_{c?nVSoGf78ViG?$5d2a2%K0=MM)3;3?@ZwUG{Sgq_@pw@VY$iZHbK`w ze6*q=AyT0dax_)>xXmI}9haO5&cjQRAYorYAcGU45jj~W-he(La^Fw!ek_tfCJ-$T zv5t9)f>vy69M`kjiVn0+r@NQ8o)7U*h^OB68kfEzbDImQ974wk)AMd{esk6z(`K)j z#nBs7^=T`eJuJA)`l(l@c!M*>O3O<0K2akh>n6PYGHH#?$!HkC%6C262K@t?P4=~4 z3y7J)g7%&4*A83$Edhx4j0ZdW8u-6!2Th*{;)DI&N`T-4_frk1L!QJWg}UoQynXOC zv_QVDKRj|+Jcu-Lwy6kYGqK0>$I&oq{hvh% zzK5xRFkCrFNL|ozA@2TO-LC;eWtK`wB3WK>Ih1DK9t;V(Dt3BiM`TN$cF*sBW@zBN zP|037#3r6}YWQ{I&Gk?)66fU_%o$DpF6JEio4T`3D^#S zb0)h<@MbQ|gY+2;)IbAJ!w5i@UAq&{o{?*gfMU$nI%kiUX4eUU(hDlAOc$#r5?G+y zry%rdgzxF>&e5u6f6hiW9v{LaSXD7TNrLM3(_?01wZv@odM+UK_+zhW2NcvKwp=8e z<_xk6GY2Z2hEg|)1T@sFV^RD33_tkA1H3cp7jEY5Mqps>ae4QOu1g|qD6y%Lk6eT` zlgCk4GM-qW_@LN>=MZJyKpcmxbK@cT8q-?#LnRJ9grp&d`S&XI;i+^a2yrf+T@Rnf zDH>(G7g1Dj*jzyY!ZVBzmz5f{2hvr&G~xCN)i{6t7PN?Lka~Q<(eu9L zU^`?fz27tYF7|i|xTORCiBQsgoTb7L%SvoV$_OVD#~m%S!!5)8 z+bj={`4`oH8UBXhXoBv==ch&Jf(5aq zUZ|3s2ya#$$s*Z82J?3V71hgw*Z}KoRGDK}@Q=^WGo6XCCDpQ>eOYP1c;e%%--3~j zCQfxrzFbqwoqM)w)MM&X<{+1(n3f`jj2C8!#^4~GS*Gi)DNlF!%?(e*HiT%SWv5!J zF+M4*6uo4?&b#lQjn?}#GU;%R>Zj@f8@DcH6MLa+)3PA^2`I`fC|A=i3v*N zDtVtSoc|mP0$f#Cotd-Oe^!Bcg4eAq|CZG~%Sy`_K>`Uo_zzC;QXSuqmoI<^0 zWE&Ww&YEZGJIqY6xtdf>*3_SmZ%AvaXXgYHtP0!~IUNN>oU1#>|53?cYXv=+Jt3_I zJcCABzCuKek(Yf0kC9j3sGA73M+%e5P&Tc*QHn;R$rC;RX%OFOOOP@?HQfgFw0SKC z!=cZy3#&?y6(F|iXyu2|hcsPctac^ZDM$><5epolp$=+6XOE}!Z0O$WaRN1?F=cj8 zZ2O~a2LAAmk)KRW+7|A3h|nZaj6?f!%y0|d$ac(d1fEfknu0a-;5cJ*Hc(!?8#&&O zT-Fl{IuQz^9*lx|P=okZI8bNb2h$o{o0fVf5EUg|lsiTZbi9d)3F}7)hdKf(UWaXE z4E5{$!L#6kZJHJ%q6sNW)OT-^M4AldU$MC~x)f{oNHRe(K*^P$?924c3#$<=w~i2% z@$VSh*gbD6%zR8R40uew;NeJ%m5D$rrw!UUI;=b+? zRD5yfA=Gg^lZLUri?+qwl=q5+q^BbFA?Q3kmf}m(ek~*f;-N^)gdK2^NYrDS5rzw< zK#ED$z4>+(+tBtt!oa0wLUW>Bw+S0ayARC~BRk&7^j(QQhqABPFNxzwm~JgC67|A&26s7dRh>FvHGv#6kdgBtz1mc0a`d67U9r6_}P!29~e=R`k&GI=*K+!h!(&(7DV z+St%%*lwPOfqFzhXIu?Y# zEz3~3qS)F0x!qov{xW`ioDn0&`Z=~`f_DNQir{UCR}wwRe}CiTHx7pF22zvSf+m8;&lC#W#%NU#bjLqv z!fg>fy(x*Vp^k`%F3$8MR@+PQ0u=3qGEYX{!iaz~u1NxH7-jx~fk-P>?nW$LhA??8 zaAtrFuDzaGN`xE6=ztum8x~MoEeQr=W~Ju~Xn=A~ZqkJ@z*A9wY8>+`9W))h9s6fe z2qfWHEp`9&Z)n!6VqWA1+s4L0HmeJU!=7?kio){f?MUS`&7WG8)t`) z=kP!X;g}bKRXgqqImxO}I#66_#*kS@i5+(>`CjNeIFgc__v~5gGm_%hoh4IwQ71qJX6W>BneVpi(z5$QD?RDrsZ0<<7=$x>u^ zI`*r9#^~E3fy-(QSCmj|6HjK%GsVF=$a}qVgqpogkoiCihDrqoFkhC zyUEDa1W$Vja@OW8rS2m2J(_x|&)-Pz+3N7yKohuEsglry18BXLRDIuR0w_$bj#$Y& zVGTAh!qUL&#Di{AEq7Ru#>%hJfFHY0ELifywz?CY-ft@%UC|jPSC`~*-26Ud> z6`KFt|G``p6ZY9p99t*5Xpo9A8oDr#p*h(FHh8??x#1` z`%P3EEvd-0PPE4?bF2r89PQvS8lz zr?>%J?)fP{7|(HkY@BCmwtpPDzMLj;@)%_xZ(8XLdv0=jdEK0?>K5h(P=R9n#0&sV=v2T97k6pH%rLm${ zRh$<2g_}A&onP)g985)3X?O5gTR**DJ|A|G5M3n!)v_kpTvggQ`=nSr3zWuIT-P^U z$$GAKWYJXN>1=3cdh>R)w{&>ijtxY_AAuu>Of9o~zUVhb>vnA2rupFLly@{mV!a-0 z&^CF5Ke~LSq}NpO1jX5N2NS$jKlJ$&(aIJnx0lO}UNTjBsjjpp)SZ};!WZdcS^fCL zcAyoVc=3R2Iq*}r`*RY(-m|!6Kk%DzxufL5cyR0L%|yjf0-SU8D(n5kxo+x$RNd!m zf2|6S`!x#?VC@4T=)WwU|Fio0|GfJ7@9YHzdIrY-S>qH7^;BK6UK>X8d{P}Ja>AEy zj|1`t=-j7j7poo=dD{npPQ{X3m&2hpF%>>KlYP)6pZ${ofF}Hhn|onFnR^WUai{Ng z?BW9#dU+^{)urkDbPoMpJ@d!!IB?dpo4R3cjZJS@H4;S7uMegWj5uYFhIiPerd#tC zZE&ja59OZ4DOJB8zbQ0SE&a;a49|$dYWrn#qYun3_XZC~-|m0ZR1$n5H1epW>Tc1#^}{|ZHeGK~;!y{(I;&P~ukxla?l%Go>cPJGd;zRkp% zG~Qb3O7rqgJ!kZWPT$A1odv|nDA;gAAm$+A6^jd>g(8Z;;!yz{`!gREhi`nwJ*#oW z_A`O940o&AMjplIl+$k8h*OgLQB@->k^5d&TB~NtX+R@`8#Gdd!*&NehJSCSjBQhP znN>Czu<`z3E&M`aB$Les(`Jy@iwb}+lL!B!K%N7m=QaI@);eAO6%3I4%ld$o% zDX6F)!cPIACeA)CUJjflB!Z(fRlbh}f<}RGs6pPDtuPOR%N5I4Z;dKs&d2?YKc6^G zEfL)voXZzvij)^U3PC>nGohF@_A|e~EGW(yZbpsEmoSYywArwqJj^}Ma-{osXt2kbsxL?%iPqk6c;b)BJoADG0J2q##4hQ2tfA zM8Z}Ax5Qo-tw$aKX9fb!?B%I-^fQe+xns_}#iJ2pj%s}Fxwg-LcnWAqPiLTxBN3O} zi7nfaw1GX&tHIzneC>_bq{W3Zto`f-BBXzy2{--aC1r?Y%{iC%xOghpz+Tnasm%ym z{-KXWgM>L_#^k8FAIPZof$B9iL6G6Y3tKS)zfuA z!42t0=(Jvd?BK~5uRJEaWfX*UEa91Yo0e3@-v?P|k3;V`EDn0Ap;qPf7zr9hie?lX z05_G1Nt$BUR!v}<;W<*3hj_L6ty2=A7^ZJI$ z{H*K?qN&->m}mGl=!BTRVn9|@j-+ES6Z6%ve^+Bv;V$BflpO)#XL!|YP7h@fwO*n^#JqJ_p45rJZGZL4se(M=z>tw>>U z!ayW{H7XIqlVl>Z_MS^ef081WcMsP-XPr0PxF93hC?bvHII^uEpn++I^%6C&gQ_Qk z<)mIuH%T=yb-`&8e^L{8nQUq_cU;)xG?xweR;vg98{FFQ)RCpNhjDVQFQ4SAyou9@ zO>$%8>mtV9!M+-%10iE#=o*y)k5j7Pxa|+!x-D+#7l;{Z4MsyQLf3WNE-tUS^hZCY zQ%n%LcdV|xllQI{$yPNEV|z_O_9`D~wQRitZBbk9aLwI0ITqiLWPYu2F}cr4uGGDR z06)af)E<%XNvmUG%<)wHsA5KhXKh)sQWb8VV7Ucb_jpH4LKA6x8U4dBxD8Rde<^r!F_rz_L`8f@3~id<<-3WTXDD$2!+-J z0k8Y7E7osn=JT4_%V4*2w$n60tNf+r;zeBC^s#>sZqXSK`-5-+;X*IHD5uhc79NVt;&a9wavnKZ`5E@(DH+n0)Pj~;YFqUKh*DPtKP*^xMk}Y#AuT7d zE?nP8zHj%CIs?&AN{%O}LcwwQRio4Iwo5SERRF)JM?OXyDP$HCju6DL>a>t$lO2bT zenLzUkf^)i&~u2MJrvFadENXO4Pxm~O!8=j0IE^hy5G+*p)--GKmk+e%^iqg78wpI zVr7IpJS_`i6<8^v*a#OgAZIq!ka#yEQF2p|szbnIR1(A1!=;r6c966*q9y-_rOli#C-J87^ zhi#*0w)$pqhqTqp5Vc?{Am^39*Gm@wFhHk^)z^QUkudyEJ{JGqgiY)3R4!T}BWE)M zBLz``|CWjPpIHhs>)+I@|6`_-p6I5e=z=!nV>6N&olxCp0IT6!4HnI8sTYF_A`c-B z0f;0ZwhYQI1d1Xo^t&{HVbRJ_L{>P1rqZASZqe&>7Wjeg^q8$4Z7DJXolP`D#-woW z>c`a+`;*Jn)y7uGQ^p-YHd#Sw!vY;NK*VD9ltlOd8Y+r@J|qAz6h06l0E1fv4K*hL z_;&1|Wc9Cg$oH^JZSZKb4nnEXZo2EQQ0k=dB*;UpAiW$|DEbb_Q+@h7Os`M>hw%&-BIyYAB23Rm@AV*-S8dU^zWe#~$N zzE-7m$rWx+&c%(lWb(g5+gR})sd7X4x}sH0yf!t?QUionfpp&UDNLJUsB27pxB{O|!mc!FBOnwNDK>uBYAV+Kn=#)QtM3oRU^mma?4I z{K}B3?ElS+yQ_oyeCOtBLNONaUVs)g+*c%3p}cdkNI4JY)PSaF@ui6DW3(*G#Jkko zY)}PQwG3WOa|>|W1PyR73!qRByXCiePn{48&T!4VwX^xYJ_Y&uq`1Q~m#0#Z8Wv}G zplhGkY3)6AuMPy14Me&P_)`}?-H#Ii$Oi$am!ByFd;}kN-xsJFWtg9R8i2bPoH`rV z(GRa0&l0FlmlFfHt;gFIAv_?e2l86f7S9v1(oaJdNh%;%7mN$E3Lk=yA0`A^ogXgx z=RW9*pguX$pJ0U?s$u|2Kc#3i1_A16sABNb;Gi7bY3t*sY58LiN4QoHFHmWI`0T6f zBs!ouJv>YBww}&vj7UGlZ5&s~a(|0$(q?3BaQL3QYo2F1Paqz6odB|(;A<)%1gu_B z97F;_5n?e8g^+jx_k7e75wI9}Ltc9r_ApI;(iqWU==czOc^K5(9aE&#u+bT;dejTy z;2EuBy?Z8i&v1Re+5ZfTz2s&w331XxdAwc@y z$wE~4DhO$z9sL#k!2NsulKpM{>PVr6!ncIN(X-+Vv3~Qys{*T{ui=gd|3XhntEJ^mVHZEl;owJg45Qhg%#PkTNSKH>5-)InwK5naOQuMV^B)o% z3Lg@o1s)8>RvVV9EJB=-pMhBAoCx+-XcRs}z4Y=YPfrW@dyW}0m3^5}ntz7YS|X}c+Sr}X~u zh5yCzRS=jV93D8UA2Hw^HikAuWkr=ooshVZ2u@X{-dW#UzFR?A!J-zSu~8#Z!Ck>v zaaPf!o>k***jESH6tm{IhEv~S2EFFKYP`a-=B?||O1u9t;>fh4c9*#eiN7mF0KFpfQdWZX8ZTBqt6Tjt)Cyx_t&mIvm@E zT_HzU##lOdx-7?3Gh1_hOMUaXQ_G2j9co)*o3%6d+3UFUcH&jzsMY#KXSV`ADih=* zlyw{9dfvLly68rgI|6q(L9U4$$*~M-NgBx`OQL;;ZIs(Y5T^fGC1c?~e=_5u>Ul zzWo>;`PQKNyy-ma%;e0U11&j2}>ESOd&A|@zrQYL36SMNf$|zSi#sNk}R2=BEuq>BG0(f zlu}vlM6e=~A}gu3j%lYH=N(rLv7|y9sjfJJ;o^~JK~j1~A_w9-TJ_xW#PYh8_4D0x z?iE=FeTP?j>x0%KEC<;``nx~3&PUv1QDawAS(^K8rYpI3po^fV{zQTKz3sgjfpuYf z=r+H&v{01hl^wl8b;bK)9V9D~RH#4llhvty@+b+AkI=3t-Ud)FXi9nd8hntN-w%EL=8gR<)4j&CHfu%b))yA@IQUZ zN2@rp2`nU-L>q9d3LYIByqp9dNHUTtC0W&xRq9HNMY3|Y)S(w>O1Sm6cX$=i!!HLD z21{CpTVq;}9^4*+v6f+1gPACQ@aBTi? z8p;f4(OB-bQE~j4J5*&-eH5(ic7<*sY!STp*<5M4A*5rw)%ei=u==pXDrB|1Fm-8j z+SNvHwYLLkneDOtCOwcg2@3{iy(P6_xz5EN|0)ZR52*2(+?+z35}op%YMQ!`CXiN@ zj*}jm{*+;uF`LPkS($~K6_)jyZI->9Bbw8gOPZUQ2bSlYcbKn|KUBb0P*R9h7*_aE zWL>mftWexv!dg;Nic=a@22|!$c2cfcK3ySH(OgMgnNx*U6;=&U?NohQqg%6BD_PrH z$5vNaPg0-WfZXt>5wOv%@w&;hX}ejyd8S38rMH!{wV{ost*D)#J*@+^BdQas^LOW0 zmwVSkw{7=%k4ev7uTJk;pIYC1zkL7XfaJjNpvYj~5dTouF!yl#2**gvDBEb$7|U4W zIP-YJ1oK40B+F#u6zf#;H2ZYx4A)HOEbna3oY36hy!iasg3Q9qqVnSMlIGI(vf=Xa zisj19s`KjGn(sR3`kxJijl@ly&D<@Dt*UL7?T#Ixo$+18-L*Y~y|aCX{r7|4hp>l< zM?a2Ajv0>IPee|pPc=>t&TP-#&Vw!xFVZi`FYB-PuEwv`t`BbPZ@zBB?=bI*?wRg; zALJf3AFUqWp8h;zJ{P~Ry$rpoz8<}~yo0|de^7k1eoB6>e_4Hf0l@$~{l0sH`fsz% zzdJPkKmUKx{~vQobKPwK2s!}3Fo4ki*5cpo3;#b_q-Xu!JoC@yAbNZzMg}I9|NRG> zrZS{sq$-y8X4|_DaWVCn?tYXG8Zi=QPCg-`eprS)k|a{ifWm;OS%v7isGxiKUb3kZ zlVh<};f|7t5^6)*qP=-miG1E5&druvE~CThOO`{N$JKWl1kaVn)|)1TFT!|$h(gw9 zD0MG5d`eGBpXExPiPTd`F{c$^3A}675!#`2`eI%^9)c2}E7UFN_-F9{+uc5}T8yUPnj97i2tD?v$tL#xIo~!r@K} zjadHssILdK{nB6h3Ad&e<+FzzsV(;7QU{gF%XHbQrSKHjEl_P#(KOT@^7iDQyXkwU z4UHxRF27e%0p4OPA0B?u!MOxnkgjd;!B^tK`($@(Ej5^ZwS6(+V+|~ai@1_c*r6O8 z`)%!_-Vgk8?4$4uIZY-Z)7=(TJ?EKKEnJwOema;3{J7h0Si1yWcV0m0LGM@ADBc|3W^oDkX$oS{ zGPkBCFz6vynR=qs7KjrKFKF^d&)&QQkNZf(LQsrI0i^0A0=lKr@;Np{dh@=s1FJULw~ zR^?73>l84s$4gN?!$*+qBbHTq?VPqSLnKQ2#D0ic+plD$0jCm{ibeCv9dcDnsLt}+ z?%P#d$SzV+0K-qRnopQ>qlYmKK~~FP17Ee!BbOl_niGyiUDJ zGA^7msam?E8MwGR0L1~BXclU8YeI1CqlUTG3ur9L(v+ETy`Ul~(u3RRQK(dFFX}=& zGhE8r;vC(;ehwo03dP*9xn~WRYLNd7L3-hf9HwWdVci!eQ4%DSX^_uQE|pMqGbgA~ z$xuO`*V597T6B2e(*Mq{qQ%BANVnoQTqe{wW#5UoGjBeOpo*?`(CC>AvO0=HM zI3`i==hLr-A+k4&-$!;Z%50~;GxDQKa?#f1Hdv#%8sXVBU>8eIvar{<93vVaR7qXH zsnA!>HJ=hwRAyZeeIUA>H_VljR%Ub+C^#i@;?tULaH^M?W_SNXGt+i2?=)rTp{-ET_z1g`bzsbJoppA|q z-SCRMw;7i8iNsxGm|9_ATx!vdj@UX)ljt{|EDAe~d=K@JFp&_vN~UU9f>$)uX2qe! zy%1*<{VIfYDwbC6Py-|5T#9{K(`q1KeOkBpA zuZcT%_u%;iPP`^MZ_%7kd&dcel!cPHG+aN$|)8jcpGsZJPfBr$*ruMsghDk zYV`t>IWKQ718!hSA{r{6Gpow21pJ^kW!A|kknjvuCoG9%an#iXKzajKu3R%;&r~|1 zZsrG3p0~*hhf}F!=`n&0SpNcfU6mN*F39f{z&D*;3XRlb(0<0+BYa8zT7N#})s%{N# zDcDM3Au97@v*+m3F|C5o#$-3@seHZAHcw61JLig^nmLUECuf_{C6_=-RM7Y|d1o#w zJtBN^6O2}5qTi-@ljKsMXHF%-SVF59q)2BP-vM{{g|_5HN-&r+r$y87ekt*-RXCma z{<&?Ny3YJoT5miJFnGS><-Pq7+{ndRgVkKC@)+7g98wxdKdaR&8x0F-YfP73JZ)_c z?b%YSkp$Bi)WD)bL!+Si5Us6V?vI=_UxRdrGf-)2W=P4;P}$1X=NoP$A%UKNMH$odUj zY$4hR&Og%bN~0onn|EzYT(02s9XuLqrdUW_s4UF0ATf4*ae=`f*};H&ests7AHM{I zM6m=4gvK_8mWA-ag-9$E=1UNp2m7Lk=F2Ye8a3jzl9J}9??P^`WCg1ylvDnY9+y#Q$x2* z#Zq^9a~e0QK{dE9^aZwYmc>KM0|F(T&6tg&EJr0gY=nyRv!RZtp+l_G4*IXMVap=rx=3ZE;D@d|>HIYC{&{pBC&CJlnQE56@~ zzou=jT{W6u7~@m%)4SXbZuc|&G^2(XNJG>kV?DF^FK?ZR^_P|$d|lUUbF1xIELkUP zma`BWTh<`_y1nWFJ8{M8seJGH^msP=QpjO*^`aby>MjsqRifk_`HA(yh$RQLM#dwk zFzrZQcEpkGJpVkm85LfXO?tUwE;A}dv9N8N4W#x4TEWXT92i=SF%x)@9~a#G`?frwkHOZ zwp>|{Uk_R#kv5uh>8Qj$-$vrS+C8+G)A^h%C!}GoOvsI1`+S1X31L*Vt$n#K5`taq zG5m*qC+>BjC{Kt;p>zap1#gs`_|4NgOF$4a=7O;ndSNs$m%LFBllo_9>nhf(=SpWx z`D2p(DvuaUVmiV47jfVUZ>Rq1>uz8V5hDZrkJP&B%anq>J0#j$C$pS(qRj!x)t{T# zL|S7ZF0v&|8)M^BD5dvq_fPof`IJsf(klg)1-kvu?GM&x?t3Czk z;OT@ii^!WOkGyrKxdMI!ypO`ok0)$biJFjm^YVNzexz236+hpm`xY`*;`@l#8<=+Q0FE+TSk^q-wxhMa1osA z_ou-sv?2c&MZE`qRH(7SOcz))7!yDT9|Q)-ZuhTk4L(*!t8JsqzCic**A5LcK3>V7 z&d@^ui|({l;nJl?1q@HJrkYLthM?1)DZlbb4d9KsNiwS7O`MV@$i2r|i2bn{ZvgTe zyaIDMUb#P?l;d}3YcQ(?ol_o8JiuK)oXf#{uG&r5Tb=5=scem0JTTiuI+x{5I$akQh3msyfu=@$ zXo00R=#}0tvibN$lszc^ulBw=EXrp6Uqn4O z32Et)4(Sf18|en=?(W~Jujh^9d*b?i|2!95yU)(tHTT?e&$F}7%xBL?7-{dzBNWEhOLfqFIVKrWtf(biRA`st!tOp=h;E{%}$Bqk14&+UOYVGc>_= z;m#|WQ@XcLV%_jfADz1oN_FFP_BMY4rF$N+)m~p?c}})n!4R_<_i{7+odD*#g#WYa zc<;y}*!P0{2vQs7S82zqsj&sC_RcU<;|;1{QY3 zVt{yB9=7<)k$>!Jx`0g4wmd`OI=ER;+9?v4dtIxgpyydYyTE|2wN^7*%t47}!sFhY_x$n_Yy}4u{gw)?PDJ_e z+8x!J!=BQe>7QMZy}{lO-KRRPTU{cN;#ZuWBbvF0DZF>P0`~&^>9*3=Te=pvsxW^C zO2RiWk3J0-)0Heg)$*#W6pf+Z>lVF|?5q(Umk2V9&GPS;Ta|Qfun8em>zN}ih&XNxRdg z?Z7py zdSWw*jacVVQfH0i5v7h!lz!3@0 zlL>n6CYLuT!&GaWHOK>0U#LocIh1VzO&XKI_WBrN?QxGQ3+tc|=it5)%k`ssF_NWB z<|QAx6*RXJA`hyi^@U>=T)ZHp(Hog5J_OO|njwqX@XOFM1a@9lP=L`+aGU>GeV*A!ncVe(9Hd$Q++eUL2pq!|p&E|qoT8)AAiJX}liX?J*bZKG zAE&A3{)4yRZNhOer~Pq)PvgFhZY^N%yF6`P<4F@7aZS0*S*+fOn2E;WFxPB&u)-p> z&#^`sLs1f~wxtmDHg_+2!BQ|P092NNwSUwqdS!43Px0Sg3axVpd9>C#V^W~@u@zJo zzJ!yhkE)(>k%v29#&~$t3wyA|oIcImtkL1Jz;pkF7+Zxe>Y`Q|`MIn%`63 zoWYnq(__T|hBv_OR6M4evhC)GsN~AOqZc^6kg{j1o}lQQF*=r6{YnQlGpI|AdQ`An z#w_Y1cLt9(-W4k1_Di1jubqSJBV-B}tf$oND@lzNOzeI5Hu*dW&R+ToQ+!YOL#`Y@ zESF`UoR|hUd~jswVm@?CK9~`kEf5ehog40YG*Xc{+%~R$U(R#avE6(2!2ihi)WyQh z@QRlq**rc+F(UMAOnY>k3S-Z=M`RN;!GtNx8h-eW5x<5GQ?@+RfxjxkTmPWmm9yRK zz?Xx!(i|Ael!6DRs*Nxws|B`(JK(5`jE`Mjv_&RP=gnHu;EQcu7#RxA=bGo{W!LdJ zlWi$hv$S1&t_$^617h;u(e@UY4$tknG7dH2)5LEH*&L&)4Tu&6UkZIL+j~%bt?o>} zC;UP8hPhD{7J<1zIK>Xek`qP-fL3o+6x`XtLezQ z+jJ(>ojWEI_#^RKdrzEB>zG=c1k=hzXy~J$#@QhQsk*kpEDKY%F#Qgh$=)S|Xk_C8 zmF^U6k6w;Rdn!n3wrY*Rx$JhSt?AW>&%8tI*_2KcS9BD%_4B(xbHZ_>Dv)8sYLxrjKmBVHU4c65*X+c zq{ltylXoup5jy3fYn^t?>RgZr*2_D)pXmJ=7M1d2ogzEf>%7jlGh)A5hQ6rvj*^PC zl)f5YTmaW+4h6XngwIWnJP&-<@azC@aeI@BU}>6O{}8`7az`b5fgGL%3Tu}|zpPRB z3To<7j?$zLc#}~OHAlEqrL4anKN`e6PWB|&omqdhIt$3=-u>)-T>-aAr4CK^P|FoC z6MaDFG+;aJz1;3nz!}6{lDS6jiILmk($cQzShKc?9wf&nSqEOI%h{^aJ;;ZES>RN} z+g*~cRcj0n1$|@`JHAEiS_O+1VlWxVDk5L zrgjuw9c?gd5$+ZaI9Om^l5}J6l#Ie-EZK7%@Z9g`@+Q|hoTg2q&xt%uIiBV@GMpYP zDOyW6uhfV*=FjiFJSwxfkh`wSe+9_DvU<#DwPfpgRvLYd|3UKV<3*6^dJqK|^(Fk>=ao>1m1AHh7HmmK{}$droo`@`_Uq3X&%95-t7YtW z`}`b|?(RkIS~7(IN)+7ZNnY6PJ{rMN{47E}C`oo37Tl|Zw!|=1IdFT5>_zsZjl{nE zB(78aZ1W4L-Nd+)^eePOb5a`~j!PQX*anWhV^m&z^ktg@F}ss{0z9A58Z){(rA75i za6*+nmOl2)W>FNa8C#LbNSF#Bdo_l)o1%1)K{XJ>cgG!hXMVd==7^p zk{xN(%hOhUcyv}N@Ko?& zsrcS##o$%$h^)rO7-(g2^8QJj?yDN0+4QB}KAB}52?+sNZTH!X#%WCCjBD#B{28J# z60b=Rna>4A_02gt#lRuzGn(~{HyDMqbfJ6aik2|qFhbJN#Gz&#^5O9Lg)B(p`5i>x z?JH8k2j{WJ2N<(sr1avLy@NGnnr3yPkp`l~k+jx?Qt zecL(>IvN$as<^v{a2(bk{)?>bs@&-h;BV^!+zHfbiJqC2Gw#f5m+XBp&GC#6r`Ky+(0aQls z{Nw>x^<=Akds}|D#jTIUV!;I#iZ0Qg`Z$&n$5E(bk(4KW_1cx>7*p~^_HvL)=M2Dz zggs|Ql-@UYHNX+m7Hcu~o!(=J`;b|L>;aR>9NXhBT&UY6GKVyF3OOmeoUc|t6wLu= zufR+v9G6#4rA1SwA$(Wtgwc9Zpr<6^#og)K`Oi}E3p}3myzY2TT-+M6&yor(Ptdu1 z!D;eDtp8@8;sGIWs5>3V=Z4RBg}2f&mDr|$$@uEvQ+#wqe(Z${D)TJ3`a_D|VmTbj zy>N5*@Mbxu^N_=$OJ42g01xPZsmg3$hgSn6{iH?Jo1GP;0+?1=ot=FRl2s&_ET*P@F9-Bu>NX|$s3Y09>u zMp=tsx{y=7@f%n*ENH#cAVRIYa9@-cp8?= z@x_+#DMts13C>5C>cE}K#=XFnM;=0(M{hUh@Fyc#%%_%!Tpn!BXK$OByUwDm{s zc9EI89K32K>$bFUJLIT+Al!=vDLzSW3wq{Kv+73BWJi7IM))yNokIVMl$+e^_xyUR z8zGu%TN_1V@JC0k8LQNnJ$~nL{v3k%LSI4tL3i5Fg#d8#u$vty9kSUMG+GGg5#wPS z6^!jPXQa5!Ys-ynNULUYI3j69z4vUfU<1tfs`=L0>=`P`_-`y0`Y-tbNm@2C`l_kYd*GeQ+{1ID1Hf56s%Ikt8uWxb&AK9&%B%?9p z#{qVWbO#kz%ihbSs#9uiI_^g{)77T9IjV1C>CR6zI>jM3_kyR|<2o;lUh?BdNV}Wi z#V%F^AEj$dw}IfVZ=ndRf-oLA<~7&8(v>4AE>Mxx)NHMq8ga5Bvkj1TBbb?Zo%?y0 z=}33o^!;t!oH&9KIf0jM*Eb6fRvR56>Ivdqz_-8&UTm}G9+{M-?z+j{0Vb;-wnDQy zHBC2PjXdaHShkW{)ulgLop(JfHV+fxtifGiQgER+4@?>}#azIT?P3=*U{rsk=EFZh zkw1iO^tw#E*_kZ(=t7zKU`RGXD5U(f*8jjq+OtgQRA#p21lV00U}sJpH`c20;6^`H zSwGB`?too>vXw99T1-~JRFp{}(zyM+1-!D8k~%G4{j*hu-YXX9tuU9P7q6@`4;(54ENv=xWA_!+DvqAoo7@wJp# z&T64{ctiHe(bA<_g{Vy2ZFM@|U95!;W*3ZBve#m0J$Sht>yW6095+lNqp(?B8G%$R zfeNj0$+HPPU`8-S6aD-nfpIz^0C1&ZK75S3?xL^=tqyxJ2d{*gBV9JtX=)?9sxb0i zLdoLZ`Sc>ojO3idBWfwP@?Nvbqkx$wuWZ7Cf-Z9h&G5SKY_7owNpPjHYT?qjXFK7m zRj|+4$tRW`TvMcT9X;v^@yln$t7>MWP$%_!W?t6&3KFJ2Wimer$~Fk90Iq?e?9SC_eXbkreX z-Q5(A+2^?Nk*(HP$A;OEt9FzOAy7gS_;G=FPDN|+Mm%LnV-eeB8MXG@!?ogDmtpxF z6aywwYa+BNBTfeo9-ZPoJ26*f3u`XNv zNVe{5_7>#5vNgVRl#9}xNwRjxIZu~ci~#U)*&e8^N(W|V-^K;P7U_|1)Ds~&q zIZtLE*kV>Uu(LXnawE1(WTu}Blm;19vLu>hpB(ZH9k3B@b|+SkhV8}boPSM;jdc0Y zDgcmhS3K2e8GRq*KGv$IDjggdR4PD@L!jbYCe2-lQW>qTp5(Ro#!mj%Zl`21?rypx8&KirhA=)q$|}mUZEeHJe*^|hSfhH+>hxl zL7=!|7aH8K^$b?}njsW0g?Ik(;Mm?yGTPBy2hkQzs_xajzF?raUPt~m?(#VNa+me4 zr|a>a!_yHQ&xIrLlm~O)o+esJm3cR39cmlzl_ivz^;NCHh)r*>xo7!1D9WV@>D;N&u0K+6 z->udW?tHY)l{ik+KohzqNFk<@QtwSeGtIQe#cUqC!l?l3$hlqSa9JX4oEc`_-5eLr zLUNzB(s#zUj`ZgD|I-Xp>Xw;;PYsx~VE` zH@DP*8M#Iesl|8t6Yab{Y0>kW0|VOvC!C=6YCV|??}MDk`mfosNlJ)nt=YpGR)V-e z)Eb?jsm&+xTGA(2vzJXF;iqqn&=$gcWWpnAP>R<+d#Ly2={%f0LJ;%B`qyRIU>Mz7 zK%~o_CyEzN+>G+OVmR704W&IXJtf{qa7~%mc$M%-q~bOi@ZEt1BhEQ{8@rYdmts$- zti&+qoMcJ{+>)|Rf5X=jvs%xzM|k9Xnl#42>0ca49l z6t=055FBI|G9Tu#Y*?Z%@Q4l&jJI^9zfVp~Hkwvez!i*zSype90o#e4y*dcS>{ z*D-;bilEyC8M1gA6y)J|t;GjlFT+wOA78G!J|!SANn?^};Z@Suh%nT0V$hCq96Cs< z)uwW5wL~cIeR%J_y7Dndysko@IkOUrZ#mH$`SMk0JyO>KcNnjwr&BjbGh&AbL;57v z{Gm5yGdD(-_K3a1$-NQwoVgsoFdxRX$$}ckjnU09jSPE*0c?6Rh8*G0u&@!CIc|ET zDzGV3=52x&wHaP~_CsP{&@(6zZKiUntlZU+vqgvx=!2@ZYNggiH}J@HB8a2mAW?}9 z_s%QPr2qJB2cCFemHmQ9Vj1r-ZQSF`;96m$^_?0K2my}jmmmNkz%A}wA9K120Tdx0 zof+iBouXLqHh#m~t*Z!gr0!rHG1j6dmA1GC?p6@o93|o+OgOTLd+hJtW4w*+Qi^>D z>EV8UHC0rQu#2}(VSPshXCql4KV(3Rd;cjSgY51X=ws;~6o0?`SMm)Bp@f)k8tH?` z(y@*!p9%?60}GQZjS{>%+vRcUC<*N9&HL|x(19*OL%KTYm_-TA6W;HAz(`L#3{KKxDZ@)4b(FBZ=fcPoC$mpZ z+Y@{>a$(!1XEhdNO-oT96A^d3Oo+c|T^)+3Xa^({I_a^k)=!9KFv<1+3JVKuh%cL# zdc6dEH|LFX>((~}HQ^n1V@~tamTlE!H#bb#BIP{LNqg(+`*5l~uzCSk~E zxD~$cXWL>VY265o(*db#Z7NseXja^^_S18k%JW!oa}M^2vEP_!BaI=mJB{<~M~fk| z@;>vL^a%Spj8IoOws&(|^@9u@9&V7w(~VrnJnfGV+_;OiUELqo zvF)$ItBcE}BR^dn1INp%h{CY+q>v}p;^^>s}ZObQq2@ua3*&ss8?tSJF*)a z&52t$L1XC|OOfo**{g*>>E|U>B=zmseSz4M1Qp>(HhOevCbs1nGw5X5=`Z)61$lM-xpOdv19yo zt_yZcN{m4M&FBMm#Q0?0G>70t6dKt0tb04|7-1hXb)B!OYg1`MbWq+>%ux_lpH#r2 zqSW7-jDC14`X*n~$d+-b1!mB+QJu&+-XXgdRlA;jSGa-vvH77&R{f57Wa$ z^R=XePa4Z`mpK&+r6Ej5XRq)H1vQ8(&MYw}%FQ)U&YtO{(@yy9%z}^k)jBIJUksLO z(E{uP$_L3${3VS^H6Fza%SnAXy>AP<5%-9kOtj%GumJH3TZ6J7$hbYDNc=Oi+ZHiKFJ=w&28;sZ(SH8d78 zMu~^!;+NHLURb_S`=sqT)h*JEfl+K|^zDygE#eZzUIV>Ng#URT&#mntVdm!O>SOfh zt;f}!n*L4VxtAi0HI}7MV)FWFR*MVxm>&8zHxzW;Q&m;XRLY`Jjfvg|lBYLk+R?cA zDkalAjd{N>A=VphYH8ayOkkD+&9RK#*Cfww&MKjCGZ{}CgO%2dEy|SuEmJIW2!mkF zlfAeR&NcyGh9^qa>bpK43%cEDkr4^>=6?~4HhMex6AH4RuU=x}b}|4(Iq@!@dgH-$ zs1w@^StOf*Z~0{~Z^LJnDi;UZSt;v(EYpgrbXo0Q76?JMxk4av_^8rKyakEe@IdLyJft8%tmpB!Sog0>%RG2+2L%H%& za0?GZAupl*U=fZmXM7%2SJ*ihF=`QIt>Pltuv)C&uArQUuk5m9w*CCwIZv^}GQhZR zLORY)mo~H51xAH+fX#y-Slxon6K2#^1ltVBjh#iGH>A}p!79v~)XhD5wOkF)^U|k> zw3Z+2gKT;BkLZ&Jg@v{xN%2mm#V*kN6H>Z5#4n2mnAe(19bI6_h#%4=*tshoyx4T4>DYd!b^5^5nJJU&!)oigdI5(RpKJ{MfpCZ z>Ua_7GOYkJ9G{{`?zXpLIYey_beJ69tZVGq<2WY}>NGjQ6>3dbI(zlu-MuFmZlcJ| zi2|aLpU}mZv9a9S%x~Bmcqm@iRm305+vfe8gy!Aj9(ad%oy4Bo*jlYZTC!}|T~QpFV?EaU@vdctk7n3#S^2PhGYqvJ38mWH>m5^5&G@3gp~3 zr*FD^+~Bm|&6p2Mr3%t7Y(O|#u`X+@ zry|UHrxsxdhO~mN++>?|OXrv>Xt!bP6n88oVQ_`&yZa<0?R?^rguzhSFkMW@m$X9Y zaUroUX)isl$QkI4#Bj0>ln$jLt=H4PO1RG;kJAAK@&)^8;c;l&>FofYwplkSd3#cc z+AJ@sl;lRTry^goDaj}Kkk7Gd1&R&i2g=6fW2GX$cTi&5eolIgBb;)PkfL{vZ8(;> z$uitoX{d%9AcQ8XH40BD5>?ru_ zlAKYBk~B+wfJAyR!3mZ_7XG0aXx|$=uU*fCDTG(Q~U}|eT3Lu=x z2t)ewruZu7%Q>GgM-vOv=RWRE($oM~N&5O}@!o0!VmF;(%W_kYBI@{dBkRZuUH)gK zBlDeJUAAd4jQVWEen$=b#MJx=BZEh_FbVzst(8ZoX;T|cST94XcH=v1Xg;=@?+Q4H zO@+C<(MB4eQT4vEqexFk z(KR&IyL70PZDr3ru7hr#qbyTwd0pJz1}1t3Yj8(!H9fkiX3Oci`ldiXhc2IbUQA!}K&=Es~A z48yy>+RO6|o#o%|m!`&S z!uz((SOGYb=?j3)ggG1E+bt83bp{qzyb|DGOw5cd(X7|2x>;({n+df~b>VfGA2)ZY zHj^Fci(FjvfhVLOvC$~a-FcFEfcXu@V;z;d_g`uSIU~0{C5gbpf8nZcC^CQVsUIjR zY(OLg)Ngo5Y#hj~`t00)hoaK?OaV^BcJG3L@}Q4Y>toHd@=I}PQ|JR>;X$#1hMS%C zcp8**mJM}VMK!_FF4Arp1E~D*Q9N*CkhhIrMj}O0hdiT5HA zx&D&PQ9O|IzHJQiE3DcNjPy74{ge&2Q_*a<7wEQTy(SYc-#oO?gN{arM?9#h=9FIu ztq7)_wMeqx&Ep#pWHPejlP=|~)aIWMrRN~f=ytI}eXi|5u|sh{u}5)~yj@RS62Skt zW}C6T=yneU==3I2ZRjicq$ zfHw~j4Av;-oQ@e<_Njvg&qzkfGx?cD?W&;R!4_sh;F9so1jHWtaFT<<*|J!3d7#a^ z;mZ%CYi2B?s(aoSlB`hch>k5hSVa%NsmIq*+P8`j}jP0cRuJE!{v1bXA9yN3fJ?Y*SIl0EuxH2T+_b! z+zI#kBKg;My*uwlcizqHyj$6Mx3dE>TGA&0#Mhx;F6SmKi!~ZMnkrtb6-G4V29J%# zIn$TP%6*KHjVaw=#@(X&87`m7OZ|-5@ zdKqr4zYzU=%4~0c4VidS6mPfw!jM1gO8lt|$iu;+VrOXk4RuIELH>u@pTJbV=O2)e zBNBPZ9AXPFgqT2K0u)b6>V`1fGVdr3CAt{(_oM5&Fj!ZBc${#|0$`Oay zfDuqDTPPd`_$JrD5DB#4pa^Ee{Jo9LV7MgcH=D?5b6cLy4QSs)e=OC(U2 z4N@c9FZrxY;kIxaQ~2K)`;Y1$CV%g}AAJcz`iu_=!)0U#h9Cq*zzDdZfh|BBYGVtv zuz(-{Jj|@jTmU+CC=$ib!G<36&Bwpf|K2&@<;X#t9N-8e8$m@Y2n?X|9ZU@QW^C(# z#2MoUD5|LN1LToMz`(%BPhhx}lLf@smLFhcXJ`Qh0}PBTp)jOkgu%BU;0H+9eM1K` zfc?<<57~d#>By%qVQB!hK#nRi2-xXMTVi-Mu380R&5MYmlibE>-F8Z$` z{z-v|ovkSx@m&qlF;*4^PEQOhA%fBp^5Otd1B9Ug!T>-gE+ZxL?1z7l0{@2PUx^m{ zO%(Gtxxx@IGvr(Q4Ulg?t^9+~zx3cYjYiUNurl(59f*UN6W9VGBQD6shJ;|`0P=8h z{zJk)6aCdBF}S6b0nAD8r}qKkP6jYbDA)#|BJ!JJ?msy_%p4jaDj+)m=i%m4}Iy{TUlOA1jcLosEGNdBF-|`3vP=WQZXF z7m?_YKV*Ctz{JY?jE#+ti97s$j3Hs)gD1@f>naT)OPu(2Dn8}dScKjr>`_n+jd z{TLcZvApz$NtgwIzB7^K(-bP;ineTh)pBaCe@hepzBQT#a9ARl-D`;S4 zh0F`zv%$Bt17i8T>US02=xw347LdQF{bmdD=&3)K{%-EKsf5o0i8Lud;mBkJF*dNX zu%!@`P*su#*fS&ZIx%D(u^Isq`a%{+3z)f9?OZt>0P!{x+`Qh+nke&)xf{76>YtLM@Nk)LAfg#iYhD@EvNe0IFtr6eTv8@vTZVW)CRyzcMPE1)w9GOl3^N9MTec!CM z`E6$SJ!F1HfvADWpWE`Y0vW@3&X+|CaEVv**9%`kCH; zP57hq&;I!y-2d!M&_C`i_)H-NMixk@;6FR{dx-yw0AcwH`7iUu-z)iH>L2bNzl-}( zio7>P-WId`p}<-SYH-+@&|}+-zF)dG_$b>Mg3Y^oQf4C^rbMEo);e{=%An{rkmVtl?z) zwxId9=2$4}!pB+fo-?>+``eYSrzsTSHudppdqc8&izwnk6++EX+Hw>pZ%C3Pr*k;5rIEuEdj zEo3w}OzP<_>RHX6Dg%tDINv^QxE1o~uzd~}Hnc{cBb_pM-j}z632VTO&g1+bqRORd z<}JjkFry$7c8)v_myLtnF1k%7q80~#N;OWB#dxer8Lia`o6wuIxxim>bP08|ePjId z>RVPlS?&5&GJXU*U+B`^<#wTrhH35XmDpFa5=+^rBTZ1o84g_6m`Y9pY6slAeobx# zNJJst?7X&t&sH($y2hV7(@- z*LwdKUC;gd4>IBqFx&{DWPq@=0fAV=;1+h4FdG0T3-Vj45fqI4S_WWa7Zm)}$;vX~ z+8{Ajadrtdb}k7K9v)6Hy)xDF~9_;ouY#5fS5N6#+_eaf$+Y zC3)DzxOgSSxg;edxWoi?S!7_wZ~(_Qo0NYppL21sD8u2l01n>or7BQo2!Mm<$6(M< zh8VM`YqJ5^0c^UzFK(Mo85&-@PSecz@Fva8v=m{l0vzc#>(?LLbpsSk3y+Qw-;4CL yzJ}hs#QLz`<6a~R;L|mo=eLx9;T`1q2g2470WrqB!^Xyqd54->LQxX)&i@0^#+VWS literal 0 HcmV?d00001 diff --git a/eupl1.1.-licence-en_0.pdf b/eupl1.1.-licence-en_0.pdf new file mode 100644 index 0000000000000000000000000000000000000000..80b5007ed2349dda087424659c75eb30b7321641 GIT binary patch literal 34271 zcma&tQ+F<0ur}z}cAnU_ZQHhOuGqG1+gPz}Co8t?q~G2BjXv22-A8kd`UO>2-9@G- zDo)2l&k92}b`@U)!$!zR=wNIO!^_JcYi4iZYDvicUqzWg+{)I~%$blu+}6m|Ow`Q8 z!PJbOAI8Pi+04ie#xti?SIY&r9og@+e!cywtfhE>2}s;a+E$I4O6|Sq8j=eUjVvjO zGOA2;;sb9Yf4ooD#L*2jcBc@rTDNt304P?Rz<^*qpW~45$o%PdG5oQe;_&45?)P|N zj>2B?r|%f}$LGiWv3)6dYk}tWAzN?dPSkBZr^9XuUSm%{Q)QoEWd-y`y%&Aah9)#EpIpYvpEds zx}AQ?FgbR~jG6Vqn;c%|oz@ObV&S=qswv(JlMINY;+iuE!@E=(NcrXdJ{>gUu;ekr zw_r+W283*2Q9z)R?uYx(gP0n%o#~4DS(P#H?IPs4vqASk?|U$9)~>&cv~PZtPo)bF zrGGkM#pbQ(jyHQY_BbMpR?XnYt)Kye?gMsCo-%(Q6(mY<-k8Q?7re)vo7V`vl>3=% z6v;ISW^EY8fU($t`^c`uu-fMlZ{sHIej&pRV{+3)*n^$DnYR#mIMz4hmtHPX^OrQC zaV0ECoGfKv(R5Qp9Zd{Tuig?SxLI9&^&_&(gBJMwm4NjCrMeb*c&@ zv29W!nU*X}?f2(2-%GfCoStr}^c7%9^z^`3)W_iWn0uW&5w6{0vSVFhpf^138A|OH z!#nIS5vqb>fY2N1A*8R7gBUg-jXd&X3`FDIhn(#~?2B|Q=sf8S2Fny}U z9KXlpju?m^3L|0Z%&SvjdRwCQ2{*QgkvjJjxc$NJ6tBQ>YIVQX3NZ6RrdC60j46hg z$h!;1511ZJPhGU+nKsttHBshZ{m5GTtDv=A#Ot-VV0!E|HjZ}OQFSJTV+WU%IgZl) zC6y%dnQ;uobZk~?jBe62{uVmw)#w3h5+p?};{AwJDV>b8f3izbjNL>1*bA>+Ibukj z*^*SX!<>8_#TcWAady7A;au(7Dc))@E#10?2m+#+Mk&2}3rzy#7x8oDXccPE-NTGz;T zEDDdRA)(8^`J+pzCk+?O7+)LJ@q+%CVfZk-AFmR)k1QHPoMgJ(j+z~Bl_yD+0 zWnV%vgEie1u4VO)s){&z0An(xp|Rp3HWSreycAbWk?U{iW-U^F$w1G@;G($#_-Y5B zz1-GBz)ktdgUEDZ@??G!RQqP$uTj)4|2EB;X4J8}phMUvL3~BjObM~Q!SXvB3frE= zu5dQ@&(V>K-#ThS@liME=F4s;5r&Gz@QIU_$IY@e~fqZ>P%=;$_p3QWQFT>lSC> zjiKe&M{F(z^QO*Q=UX93SvKzK$K(pT{9e3a!3*&-%WjGy$w~SQ>a7NYJSN z5=1ycR@e9`C4z)xJS5aMi0X2lGQuXSj+LRc3#*9`D)a}yet~YQn~Eb7th1c@pef-bcD-&ZF9-DH!X#laLK?XV;cj9msoEHI z^Z?+CUK(tRmJdRQ^o1^N70TbI_@*wey>a|k5g0n$Px{Ss2v$1%U7T)$W@+L~k|4#x zr={O6)3)-+qK;yobAq9a=w8W>HuXkZeujVKC~01h{tS(Eqxg4YBx-!KRWjjO4EkrE z*Pc7GzWERJ)=biA<}yQjg%qN^LyJqVemk13-#B=>Q9Aht|umOiVCAZFkrV0H4Scu{J zFkN+>&DDlJriEbk3!I!nX+-@m+ncRgsmY-OcUDLm8Cpn8B)&Nl#qUdDAgz@}a=6N( zq2dr({wHBP6MEkSH=*O3;b8*YKp}Cw5%IuyjrMbuHFHgcTNMHfZHS1{u$4Yz%0a5KVAYa>ct;-e0zj<%a@@1nnN(D{*H6?Ih6rN7IE^!~M)gOV= z%E*@$dpju~N}a3PES5tN1M5;``U3kXnY9*t>p+cvPjnS9wWhJz0{mn4RE~U$x>}dQ zy%_ciLnl7}N+2g;jqFX8aaO0dM-~b=_drRG@(xpNdsHoEOX0OUa!M?7AJkz#fe(mMiprOUyFW$hRHRA9)ThW2$jQ zOpFs(Q6^m^&8);Eqb2-BOTHo@l|bg40;s{+j-+t@p-2wEpnGgqVyhHK2O?vWv&9&X zRR)fXz)v6>=DRT*%z(sznUj0~w-Ph8jcz|T;b#vYD$#U4ry)&MNh&l&4FTu3)RwC= zhJ#`7;hC0=dkMIw+mV$KTFUHwHY-kJznvJiizHUY74HP4ANx1VY^Gz!UA$J7#vP1s zJ^R?Ro6*m-hkPR|DMMJVV+hXeKu<3mZ*bR8_>j`r3|uSL6^KH&Yo-AhfH2fjtC{3$ zla!)5R3~CQB?aSr6=Slxt;o99(}HXCeLfm9a2eCq$XfVArf!;{+(gxPr-pu%5vlFJ zH)pfaiWWFd20D%z33};7&*-GA&PJayBVVJNWC&gS#Edt_C=06fZMaH#e~o!uY6;p( zte!6A{C97agIK(n=-=7Ni!!VeteyGc>Ar~^Mu$&H1E({l#b$E)`*@YN{xZ^C(6z4q zGyI7OoD$6|e>^Bu>d8)bZrDWOzv@bUG1VU#QS@@(4(&~}$}ikmIAwXr3pK0MB8BQK z0@yf1W;|S%O!-1X-ei71LQCe}rKyYE$Fe80MsLeUh6*~umJy;;i8iMidu1H%7FXg3 zG}T3yqS?xr;v8{0P`K8IO?yg1zoohimxJ8#KcJTNdmdZbg8}g!U0zgCy9(wWG zv^dHt>_)?r4of<7g1U9Zd!CL#Zj_R(555s(Wim)u3$YA!r1Uu_Xs3yzd`f_b^tO<2 z*BfNY@fR{u=L9dWu|wE0bdS1C(JF7oNtFUOtPZz<5p~kq(br0dS(9hIP$Fl_N!M%e z*Ugs_q~}5%ne)IL(vn!ybEk+dAM$Qtx!OAC8Kf~Xe6a|p$ua_Vh=^}R3v4DgmY??A zTa8rMk3wj1h!jzAyM+Subech>?omSsBUqT4{oCz&@5qt5WEawtW^{1O>vF(Y9be6i zg?F-?*4Oq3WWACfu^GrFz(9-4xwyF={})#2d84iJidIxAQB+DKVOa2hM;-tHp8wL`{|>Mr?v^Z4U>8idlt`umEz7P5)#hBgkTw4E3BibdQX)E0RqkS$roFp?o*s8356#FR$oNvd9m4TG zuF6F&kxEo2?AuLaAZz(G==FN}da|zBYrSWwq<&3?TuKqwyS{Ddes;K-$$8S__Q}rO zGw%k&sB<6ap37ZdR;@`jT#{ruqMArDBO6DXVFm3;cZf0P#VGbt?Od`saM=b@f9`n2htW)B~Mi(Rr=q>$~s zG`#){^sriJW%O}yXpEI&#n%>v;0pv&s7c;s)GTAhOJ=NKsGuN!Br_Kv^LYV>%zC3| zi_S6FV>5<$wTn(C10+`HsTD~D*Jtd2-C*`*=gNmdyu7b7@{42dwGFN^$P*$3fW*7% z_UJ=_pbj%$&bub*^SDy5A#nX0&X2p==NMQ5k&tr)PQtgM&=h>-;2Cwdyaf=4{M3{Ch)%Gu2l&qr5BS9p(A{RW9Uswq_n{~1 zySlt-1(hVM>8=Jo@|zwXIIkD6F#$?srRs6(UlF}nxfVzva)Y=R!<+0wi7pKEwIWeM z9t4bCD)hi0MV;=Jd`zw|>ERdD7p*xZyv5Xowpv_u!fF-;$BB8V=!jZk!T5(zI>t4tC zyNvR896c$SkbKa_l}e4zhFk8nIyAQb2AJt8-`s&-0^MaDG`f18r6|4EW_VP^k1suz zXu@B2mSl+)q?CDDf8?{>2iqVW@FJ+)(0g?rmq^oPZ$M^3sRcxCuIA+z`AeNs#e(`nM>bUO!4CFQ)%`j=VSFN5sQ?x#6S?gq_T|lKO_6pVtsa8=Y!5R zT<34Sj*Z46XCksuSm2{98!B^S!x>7eWr}; zm`k-X=0%=o5D!H`vI6q1_XW}Iwucwo#ZfO4RyCFW$%%YoC##~@GLuI?v$`W$i~E-1 zvDbD1jVZ@EIe}4j5J(6(9*~&+ z@nH2$9opK}rj{4K@W%lo4lPr3t(B$XomhCBVr>0rH!e3K{~b5!Y(0)&Sb{qHxvrhv zA^-YWKfdgg++Quf{``7-5C3k?&jX3}$H}gV%$<>`W_d4xb9-Iwvz>Vlp0z*6qdOzl z8}Ii1_xa5t4m@UG4LqDXxAb*OYh}%oaxP=7^(W^o2s=GAOIo>QOI<3#8=+38xkuQ2 zZwrU=o0J5rk2@P92LAj4SG9|E&Ayi2^Q%>}vZt7Tignjj{C~Xt&r~@!8*jrNF+YQq zmN+egas1P9`ZjFL#wM~KtKTiyT)#tjdPai3+XSAhaZmj!CZZRe+BrmCZ=?N}|G29^ z2RA2D3*_(a3*El%yBZTTnbmvWZUgzAtl$)19nJktSO04Kh=%#V+3}!Ln$G=Z>-a=g zjMg`CRJHThKO`v{VK)%W30$3my7UY^$*NxXp2dCgT_NgYuXXzAt5`Yr)riIO+53RT zhBLPW$F6)UMrZqkoGI!bo9U+=Qqd{BdU%4c3TMLN@o2$suoaEo?wfG_T>zpP#qbb* z-11V(eMZ&4U>K~L*Ph3x#|ewY%N%lkBAxe4$9;x!c?^GW_IIozD7PoSoK`o7I$WQY zh=A4~bK3b!R-Utj?{rAiP{}-cjUcsoI{SOX#XECOx8Jz?s(=0uf?h1;@0^c<+d#pU zA*~__f=tObH@*Axw|FcB4GFJlT@{-`;da6F(N)AIdk#CIH_z?J=0Ev;!>0=IT1eJp zQ4l}6FXUNxH^l&NMQ?uQA|UKk@wanpu2#Q>j?k#?>@wG{{ZIQWH@CM>hAs3?rsS-- z_e-A$hslW_jFUTgg`oVx1|V)9vTtk(v&f3+b57My!Au)b|>M# zMOgWPl!g0K^$|n=iWwF%8s?voBkfGY9??Q<(7|wW9ak)o&D3?D1suf{!JKN}V+vvE=OJwd`v|KSKg>?6OeYndK6n2uTlEwO-6AT3Ire6!v)2VuSuPe-HWuNrep!u*>G4@EXhlv1Sw|AqXoE_^ zuxrUAPbj~ZK%t&bWM)j#8_FAU;hY@K&;k3!FbMA?9S@fwW1n}E4i#Te#Ml_B*F^21 z3Y>g=tty`4u}i`WYo{Dx52`|m@IC~!3m&Ms(euQl5_5)e=eKK ziSRuX(Sx(jH-RNW`py~mYUHGL4y*&Z*6!ghARXc?mi6_VPsk6Zwanl-T_Bc3$f)B$ z$t|5u%=Gx<&vo%S4)F|UF~e{5PdnEWZ?ts2Vd8(e>yHZ?h?&R+!e9QOrMTy}9+8cN zG^9Y{Mr%SDrgMs9L55xF0zu86ilDPcjeij3lW$mFpLKJsLerCMQdZCH9vo5?2K3%q zK)jtCT!phZJk`C5^udyB`CpL0Mk=Z z%w-S1V1}aDu}&af+HKDUW&2l2=(ZFq0@l&r1g-=c1a=jF#z}X#fFU0`>EX4-1HO0Vp25K2wTDAkpf95Tei#_SpucRG zOa2ko@@xoM?H#Fce4U2HG``4{`cmX2|2ckP!D))*!WGt8;_y2$-9>ELjQLy{s-5xS z5AdDJtw1eW>|I3X8P}IlxmPYS25Fa(X8)*MQB-EnIee09WGP`I7XV$%&;pMub}}O| zR8BQugBAg2SagUuK@<#^jfnFMVWKLpb)C~b3Ywq=XBkPCcGx0B&YqJZsiop;ulByG zBepOL3oOO^s?w^t5W?UHK2|lXdU%fbk8io`X{IQuY;h)5ZGx*5Y>hICC+t)voU9@R zyPrfVWc37$7pxdG5k4c@l!~-+1ZVyXlrC&`3a_`bBGf^)p&*%O{w&bv<1U4J76gie z$Q`9q*WN(A&=U6N4_c?$(@jS4U!QL?a|ke@VEinniW?^G24IvBC5@(8nICW5z+-%H zzG35$VxxfNoAP%Vhwi#aP}4u(uxgZo885npBQcP%JR`S45XHDE`=+GeqRcX`F=p~j zK&IL2q;E}@02x7}@Tkth#uD&I?kMV*)8^6a<<{Jg7{DgngewLKb}@PnAPkyCt-()a z-Vs7dZf`7fi$94-OYkiY6J=00pvE?J@(6yPsjO~xK9{;_E@~Iq>;)-e9h-;v?_DIN zI>;Rx-#QXF3$FO4VZG$U&ZFYd^S4eDd0%0WJN0yeCI04c8p0x ztj2M@-dQliOEN|nvns}Q*~LG8!a&s&D97!p;RP!NN-a!vdBgu08kNkGf_Y~m^oib# z(5nM3{cxoNmAC)iLapoEv#11^w*5XGy!tC7;w>1%tHw8Gnph)15Lcz)I;2(^9gL43 zEY`?H9UHa`XUMaj6jRpeBsuYhGhw%GBTd=HK4^zTN0%?I0zNn<$+Q>y^*5QwP}MS* z;8b$kJRzu-R*iFT=o=U)Hq$U>(_8CG=$5oYP=x}{@dWTT!77~s4wH%=Wa;_8=lf}I(zUIt3^6e}KDrE5(|D3_$R(}L+o%-w}` zbW+O{cU~HL7KI-9>_dNt*jT{!1QYrEUq{2D z_i%_>Hq*Plgc=E{A&erX(5_cqWx+j!jl-W|rYLqfB&ko_>Z<3Ue7SU5EDksdd#-ex zXa5-&Ws;IhCk`&g{$)kJ1(}@yUhEN4mw!q-VEV^la_AlR!iIl_gWYO@$_p`mv?mtS zSBQm0+0{Sd{+>@Hrqj_TgbUMd7`FyxrVDQbW^rVt1py-u{2&=inZ4V#xT}x`QF+%L z!uzzgCR_<#!}^K-LCzs4`JnocF=>?!-q+&_38FSt{^yPbipoN5HROkbQg(5wg2aY1 zE<9y8Rh^9^G(1r3gfJmS$!vQc`RT`u~LFkSHqv)wEK`xp!tI;NJ8=&S` zK+!z0U)Lr|GcqP{s(OV(Kr$4a8l|ff*AxeHf4pVbxQKqH9Z_j@WTLZV)N(F4LyM|T z?U1}T6Vt@6UcjhiA%9adiQEgwR^$9kfoUO@m-MY)zwlD6w3cgSnrqPn3utt`GJRrOIB1ly ztavy@-l#DkQ`eW1jup)Iyc0GodYiJAYf&+y$$`-yVNp4X8FnD*vh7qfrHYM|KvNE0}tOi7e zhBGCyrADZS^^^7wrIsW+Twt*>!KF$NrjN3Buvp&Qx%(pE-_)!H2$b%$7&W^I9#7kq zzl)Wdm%(9V=?V~9!mkjj%6ZU9b+Naepx*W;^iot8JzX1GM8T_AnBmV~%p?lcNkIpc z_I&8%W9aQAnS7bWlk?|U6JH)Cx(hqyN`1NS;VRXi5jYBB81~(k(l;&7USd5rT`yph zR{9t|AYkaB!l;_`)kSq2j^)ItF9AVciq*WtH!Ori+iNqWQXQZWByJgKQ-schEPsHo zOizUUCn4*ms^EXSHzf|M5xR)*==uT~N83F^Jx-FD1+ha{A>dOSGL zEhEnJIEH$Us3sy1px$_T6kW|+plkeGLm^F=m;R1lp?uoR`}ibgO#>6c$h?kF2y z?M{w+`$NG})X(_^9itN6+eD(a=G5|H%nb~C!{DnyxR8UltFZl<-M(jnxnqPVbxekGgGHOy{EMgSKda^k; zaxtnRbhNKMb6!qsief=qb!|7O@r@;wHxdw|di?QI;UMMlSm){u z%sNXQFn3Z8aRm@cJ0>Tlz@WvaJnj$C8^XH&ObWGQ54Pbm1sbmgNiGg_1PMbBdDrht ze&Qn80u^IRR(wU?p0x(Xef%>cJ*kZElWechP_QNY2jeBnu1bmqX6p-1Hp!3`8Vx== z)E|+~VHt6P?ro|vMw#S9DLj2VOOF#+;O8n;3WMlTPy*y)W*m4SYHbsa^#D#YhwelL zj_gCL^a8{siZkMzW0`H(k$i;soSeh)rpJ0D2`Q+9Ei^ocTn#b|Hc>YK1~&*fdqj!h zB?=9(e`6+d?zbFk)}7Xx?WZBeQM+SnjzKMO6eVAv?~$Ur`i~r-@rrs-c~(nbnpgY6 zDfKn2z6Z6lJ`G^~%0Ue1$Gv=BK5S^s%Sx1ImsJJFPFJhx9Z%PW5*PFqi>VzFXd`?8 z6N2$anTEgpl0)rn4IeC+A!0{Yr1g#CLkv6<>}Ph_IXomV7EHG5CZ zt3+ZAqJJXIh9xH2- z8B6hsiAI<0l1M|cF%7w42TQAkL&&S)cU3otK?0v1Gaa(q$be*uz$lkH8odK1IN#a7 z|ITJFv_(trf;;=%AdcZOnNY?ViX>)BQ^oNK*LxRERUjd{SZW!mXWZtDP+WKh=69BX<;IZ41-g1W#V2TipM8yH z>jdT>TbF@PerBs+e2`;1=^3wDwo;VadfiC+dHS~vGI0&>p{xom$RLw&7z%poI0*RS z8t-5F#EmKPHqo(3ig&*GLqaOdJtT?~=uqZ+(FPYE;#GTTi}aeY2LV04gE4}#erNZx3ApZ}#y|2H4~Kg`6$`TsH#$Ny#~&i`M`q-Wr=-G=PHR^HM6&&`Jw zIsq)Bo!yG(+Vp~aK^8xWKZjJtopvJ1nrn8=@-q+?Mbi*-oy?1(e8yG0cgX*0cbERf z@87aoG5hly!IzVlhq}Yh!H>Vo%VTEk(_7y;?APn@&Pn{K1w$WOcH!bEpp~nSzq;C} zL2gHZE=z806=w+{^_=6>!(waZ=BvEfxB5~Ve#hJ)%iYBfonO7RE(VWh*7xZ_r`rm! z=F4csy%TbtRaVE{$`p&h^3!1F@!3`NXQ$n*z1ID+Zc%T^!=ClGym71EJmF+zU2V;SD$J$@7 zpY~>f+E<^^OjB3D)LOPl3G6f z*0+u0@7KrG*QkAX zKscv5yV%uhy`BcV`1!eUc$f8Op-0GyYaHWM07fwV(S&XXWXnOU)UeC@3E!xc zptCP_x~ItUNAwsdnxXCN7WgM3(jVYx8QIu|yVphnBfLEm8=O5;AAC#Sr8aafIQ6aw z*rQ|Wt5n8O1|Eu7o$IoWx&QhE0l8dx_PUtlMJ1oX1C0q7&7JB3!j-{Xq!CUjOzE{~limBLS0R#kjk}&*O956Ci=L^7PPUaj*=TI>E87(}x8`qq*4`E-JoaJ` zIXLY?mv=IIKx7Fkc~>mv4MGZAmx7Gdx$-=iB^`JrGK7dpp$-u-ozfZ6^Ar|?@>Y!N zqn(Q~cKZ1bQ@PhD_^ODwT!H$(rBGJi6%M&ETPNH=3R>2HarwM(IkU@(*o6*LWGfVJ zrG^7%moxImOqX6NFI2iiLQP6UE9V+wK}(Wkgjh%9dc29PlKO>~_MC$boZqTvLl!su zUUAf!gqXFV$eE$FQb@mpfbgOc5gac_L_x0biP2=7m65=lv`!bonENG7z7%ZF*twLOeHNRI|oW zD{Tau%r-ljTyBZ;n~+?fUtW2=a_SX(?_=|v`}3Sjj&=Lhavl$Uxs~NsS@ZLNqQjSn zT*uWw;Uk)kjm(NyI#UwZ$iev+XfFstFMb+ibE=_!?p0n#2KM^PX{!+>htjC9Ux*I2rEUR<&}- zROqN_w8L4V#nEju7FP*+9;LG;-F=M~>^W;ib(2jaE(&_jSrZ#OQ|HrHR!~N+rl)@{ ztIKqK7e}aQSOpG(va1Zj+2l1wh|!V|$gq6zAD4*G$+F)3-_IE@>FPUl9aDO1;QaYfh z^!w+MW<+d1K}R_ZS9Oe``k`*I-fBd%DbuW4MP+N`{2~>J-T?4!h65Qo+h|DKq(KRY zqN{H<+jCZiF|RRpyN6`_rr|o)5o>v7D^S6j1=#zR^ZlZ^NF20A-A-2LpCW~DF!}U^ zZe0uw0e%$}oW|mA=SoFxs%|^c+EM7asyd8qL*pd&k!5(IOKcjhvVS6t3(hyyBnX1{ z18o-q8_C)m=Bd8tlXg%nkORtWN@W3 zIUz2n;&G2OiQ1=>9}9Sd>q`$t778NRRH1TWN$&Tb)2tej`@q6Z6H=Cqz-e=|k#5gs z_0}lGT29i0x2etEfS`tEl^>26R+yB+Zg#OS#@VZCY>d8meyE8b7=u+>wRQudZI37h zrb(62yQGK&Uh?u$UYtg3wRMUgh|` z&!j2I0J-V#jLRevM5NNm=HF`dvgnWzhXl5O9dEsIyqGy!) zt)pWhQkadfAmo5i6moj(K<=~3qU{5S-~e%`ER<1bexeAvks?1)#Gi)kla;86pp*74 zib~88n8?b@AD{P2LThgEpaAC*un#GbzG+m%zg*KFYPPg_E6}_x$QH!0&4|Xk#uqO>)6Yuu!0KSL;u!t2qFLLI;Xx4Q0Wac=w^f62 z_=E$to77QA{FN9uFhHu=XTGEu183z6#SZpIXou`l#X~V_lnw zSg3fZ8>v~#Q1Eu4BPu;Mu0<4}TMOvZcgI2?9o?ILTdN|LTDMwkXX045Wp(n`7cdo}Sdz5JFczqR(buaa3rsY) zKPik^E}oq@Q9CX|M@!P9CU#Sa`E|J(FVTjArGKf}+VgBSCyfU-ni zO=VUUzfxSn8Fx;w-f&>v1J}ovDqh}HF;ER2|BVlog7h;{fSU`L_E;Q{-d)WL;ULA!$z-mtLxx}^{?Ovaidm%NDjE=bCE{*~fzVXmb zXn;@3vnX2E&vTGatiC9PQ;TGvf>aXvaEe;$A%BXI*A%^5td(3ryFpam7e)9~W<&pM zt_415s^8RRZ{?Jx`{?0l6*jaETPq6e56l#$Fp*18hUOp>;n7g*sISPN?T>hanh&eR zdx^rVu}>VWqVLM3xjKg9Iv68sIkO!w!WH;Ruf5~J*eF`_GcqwS=c&1e%cGatR-aGYC?$oM=1zEH=%7o?= z+*@HyDJ6D*YxT{KCoA z!8Z!z+!%)(nLJ_<=}L;EXOE^NG4r4p*f#Qv%Wbw|7(?WL12{%>yF&l>@>}Ge%NCG- z4i|M*a`CJjM~G5){62m3?VRRSU*_WSr^ZQOaS*i)2bc@K>@$|=e;G0%_Lrpv ztH9nA=Hd;9oO}S1nuw!mz+fvoyindP{EF}|sz{c*oi-{jT=AJKwY*{oEz78_lrv!& zal_bR&YUDM9KUU)<6X~>&JS$j0`$x@G^ImK48bl0Xb@Vl1&=$&8dkQt{>(`f&^#s& zgxfG>`@2)8Na4t+zamWzyJ~ir$NJ{6VU_X?g2HpD4*iaw%xVlgLt@Pd05YSalx`VE z_Z$O*65<%{!Aw?j)tG&3D?9Y%uqr$;2UiWA9eRF}V_04*yE2YJkEh|r004t7&r09L zPFl(uQnNQH=?zSy9yW;xC~ibupu(}B5F8?}Vi;|HDq0t$@GtWAT+zDGjUnD~W7whd zo|RK560^2sMIs>qLG|FpYf1$2yieI5zK(gGXSDo;**>r6NIhAnfYagaZ+f%PFD_<` zHPw=RF-BW=@UK8^1H#ZLBIuuWCS1ukT8~Gs7WhRKXLYm!@I;~^Mwz@2i363(WFNp8 zKPX2LcARlqrA~8L0=_zq{5y{_eJ6F9vgAn1=7Nj8#$^1F5z#yFpHCN8F~XzO_#0g- zvLgQKAl!M*NOaP?U9)dzlq-lz`GrMd&A_lOrBm9Ti9mGTYsHr#%eWQ&zQD-Opg#{N zvv4@wzfy1g&Ltm4M1|}Wc^RJST_DAvNYr+@Z|hmU)h`CzA{CV*fWEc&@Z!h0QPczN z%p!87AE%X0uo!+c*Mn_$cNOfp28E{v6LfXwDYGf#YpPS{!8)IRbxQ;A#Z)fGp2N_x z3KdY;$87GJZOLwQoBI^FiNS12R+vJZWzN}8WtQYx729vt2aeV#d=W-dlp%*>vVW$F zFUaBw1IvJ{%CmGR_ryX5RYv!X&^u>h2;aHgq{i0CY2Z^LO}O`cQCVlk)D<8e>-O}Ik{bdO+|K;-mwxi&}27Q2MU ze#Yk5um_Jf&V7c37WcBy6_j!ZKt_BF%G0pfn2LQNw6td{m33rrJm-o&1`8cO?=qPJOy33AK>tf&s))0;BfvS;2 zj=+KTntxLM)TYfIyO+m2Gj0YjHO}s+@3;#JYi5&d!+TVCTU^hrDDA$=t+*9_jYQNG zI@>jBZwVZzcv142G@(7kM+eEp-4uVY$11Nx?96n`AkBXuyhqrg!LE$&hXMiMirh2* zb_<7o7OfGE4!G()V_!|#c}qh<)X&1V~j}Oo@`)JttvGvh$ht{+RL$zMpO0|&9GZc zn^;l;9`pNB0B)EMq%X1%OI@5sutA*SWyUZ}MM|ZB7uKC%sR)T>XCYqGt`EojDs|l zO8;F2V9-)T!ua4$)e5f^2=;!J1Q16 zPAnI@3aldbD+htQc!T?CS4LJZ^JV4kd|W*(d+9X!>HW}9kLw(@z(IFfWA6L#B}HNQ z3CAG*QQu~j+5RcNg~8pV1bm#T&M;UQwrWh834+vibwOVgKpJ_*AgDNP6v>ym?E4?! zA9PRTx6#un(ktYFbDgK$88N~KZ1vUc89N`@Hup)7Oz-$g#M@nu*3hB zb0Y-cML&B4Xl4#*D7cMXPu-ep1=TTCP7A7xH9pBY(%-^`*~;ky|i z`=Qw~eMO?w0^*XVc=Zzz44bPWD;~iKU8!lPO08lsNbp(2|E>@x5}XxWE-$mWicPbo zjVsc&m6ECLns7I&?$c4DrL&J-6$ky&gsl1C)OO~LE6UYmC)lrC*utvRkRgiSjrP&z z`}0Mm!gSM=8;kWM-SK5dvpNK#e+G1hL+(8-4h9%!+`N%GuPDPKPx5tUx!-(8q^%A$Cf=@F@NmOdr>5#F=wyPDTwtF)h! zD|cWn{PmJzr=l_Vtv_u0_WMyzrfqaaNz$NZvO{l{|9c_ zSlRzyxMgPgUvSI(Kas)z&qLJ_JuL@Z4ivw&`i%BmEzeE4rtXdJj%-r~EKO1*L$Cq< z@h*<0l6Cc_AG`@fefMNB^+dI+v};1>AW0V9C{sKGyb{07*4V{v_a96@PX9lPH$OK& z-Td$8rnb{xa}VG@?_bw9_2m={ea+pKOzth-9^P(?>`Saw?wZ;adb_PWTXL~@?&S^& zFWp1u{t9uqRycgm-NkAFZSO;+5+gRau8a`hhSnT&^s zpZkus+3t)YYcHQCj#oHK)rB6#?*!K!zZZ_Gavbp9{x0*eZZ_vH79CsHLTu((6g_?M z@Vfph-5bw8z#Ki?8R8qOu8e`~pMAsk`Npq^je1$ha8Eam%QKu$(zm+=ae6FpF3HQB z3s0VBpZfykzufO9Pe)tHQr&v_Pc7LkJQ+LPQ_0ov7MfeCwJU$KtyK0&X>VLrZ_?{< zJS0l9iSRN$?VS61xw*a{5kHrPoGewhW?_YFdSq?fy2~V1WQ!--L0>Lk9CHA6I$PKE zh>u=dfw)K)M+ghZJMB9iebIC$*yIvLJQ*60I~>ITNZvCqVu*#~Q8yPZkQg2Wj~n$^ z?INVw{dCW#h_q$qZ&&fMEQxD!7VI$!L5_nU&Oq_qT`-vL0eEv9n{Ip$3f}``hdbR3 zCe5|hB~JnTMYQ;iLkWL`+}zGp_8rDHK<m#yX$S$(Q ziwR-FdZk1EvH#L{br5lxxT)YAR!dYgE-F~R&wrJYt+s%%hE$TxKbZ_h2gCWYxLZ@P z|GQ0@ozm{1Rp zYc8lGwSn17fszt(oJX=W;}j?$G4yvJ-yiJ2I9PUvwCLY#ww0{K==YSGdQdZy)1DHq z)iUwDbXNKWRz^V1<7Y4+ONh$;J2B8OjuaA6E4ngwWMN*6OW1X(t71<`IuOEWl4r?8yysKW-;lxpT#}7haVRwq*P{ZOGKFlI|5BsMroyw5eG~&^v8T6fw848`o3do>jPmm-=_bwE(f%2mePmM@_jTH@qu2vGHH3m4t z;KwKf0J@Ju`HkWj`FynHEjk5+dwVF9uw}kVLg-8Ux8J@^)~8OS!&-8}QQcx#+rK}q zC@EEAY`dbmpw*2SbMDC)?A5*$>*DY%H>?N(=sNBSbJR(5c+d+z3c_PB^3WwnIUuUA zoMcj7d%(!W)4}cHyzR}d5ZD+IcpQ@3&U@EQJ9kux9LSZbOoH&0U>AygAo-54K;*E1 zl??L>x(d~0=<_4UnEyrqbbRsN1B@iB%(JDpZ#oA?W>{6+vc|cUKG{Q|nM3SLN#Gao zM=lznQrXCUpI{ql{Qf)6EL#XO^`}44%fQP0`wwhKfg--JcbsG;pXz|Z ziB*VDoNa+^6Kwc53tFMD2n)SV$GpbbjZ!h|$G zQ)mct3bbNYH_YrH=1r*G(eri+*yTK1o|+;(!2t!Vy*jg<_x-4Gb9Uj^%;R{gC=< z=;mY^(<_{@hW@2cP*BxiU4^KG`4e3EX|dOYYz#Ok^9(wC@_l}+PpI~HVY7(Th<5u5 z9V4jaa>o^To4NJ4W0`dg(!vrL`Mt) zv`8uLuEmQx6xX6HTHLLLQd;i7{&shF_jk_y{y68}c^>lQO;*;qMa+(D>bD(g%?S;7C zb4vPPwz2pV<+$|xxyljV=mXd$%41?=F8tmk$Rhler!eN5+UeXJjaCUh5xYTRbl0)G zZ4rvAdK%*T3JtSDJ$HFy+FqX2BD-1A=_o^*-c9C*e0hAmp`Ek!c>3d>6RCCjt-9L` zNMB<0%s0PZ93nYUd|g%TiX5skFfpd9paS1hR?DYZr~J%O|1wa-BKT^C zn7Wobt;1vDBIm0IR=ej>!X+D5I*C}sGUY>*WN|l)0ZZS!lE=Fe9AC%Mv`67chfiq_ zi~%3Bk7B_z&*iyJzO2D9!S-GBI8e17R*FAp4#IbOZ=vg=)ge>Hr2Bla+808C+y%MRK^gS8m( z-%g#M;8Dg`p~Q)|i#E@D`di|tWas-1lZJI}pf_Wp7Z)68B`cwQEO{1v5?ooPVD`=y zDX8l;)q;9oDb7$Rx8VvZ&$$G|$Rv!Uo}ds}mF<2}f{&2JiM8I#s24$Yvw)!^dmC9^ z8)1SCr-O`%jRv~NrHU&l2B#Vv5*#Y^xH@{TUWfOl((ZwWqG(-8Y9doHJt!d30^%Ts zecH4EY$7}=!rcpXaUln+K!ZYFwN3**JMFYpUqA zf!L}rYKqX{iz=cn1HTBE;8r-ySc=V9A1X;pOeO8g! zo(L{l#8)EvA48a(yXSCfN%NmI<0PKKmA6$55Q=l6r3TL=lNHiAlF}6QvhuK#qy|#0 zsJoqV^i%Zpv1`+!$kU++&4aroJlgHHlOoW$qJiCRK_{IO!R_1LKwLK3isY)()An1~b>E_i+ZiW^VTWSsKy4GtdaMn!W zqR~Z4CLa%zzfI&iAvaHeeI1sl(dlh+gQUbWnrt}Lrzya1_=MxJDr@n;Tuuu66FHik zW}B9DdP)CIYETg~WSC8bs8*d5Td&e1cqP}ypN2P}C7Y1*OVWAL*3gz=FH#r1@u6cF zabc7UP&BNZtY-Tq)m&PQa{3`7d*~BXzWtYq6Qn!gkl5(zn7&!`k_~QEnu%%9GsEZ` z5EI^B(xgMrqFf1sPx`m_c?FrM`U1g%IfnrR#Cm8I_(RTf%ja>$sba(ieFVyh)NIr1 zhBQ?j%BT^cRTPNUwD`Gd!xTJEC5n~Q-VrwWD%-C1m zRnui*$LHY_Rp@vu?>Gv}`!>A=Q@h8@2xvSTNN4xz|#Q<+BmMlyRu^*UF z0;n@pD^P~K=j+l&y^^+HONiLC4&`P^^SMSeK|@W|au;}VDPl;<6H%P1@|6K>{km%` zaG_gT^J)0OBTP%{Ai(Dl@f`#YUm?P>urVfm^pJ;V6-5z{z%F)h&w24@4)!pXNCOK> zPeH}+lpLaSO$ z*N#JDUqc)Bm^{iAMgSTGS&Lte3M?0jzYM)%s8xE)Y4#*4%^-hVAZQE3I+EBmpW0C@ zeeUNmqRvA5NRpl@^DD9OlwyJX@s9dPpkER_d zMNABx4EFu4(Fga&He8#kE%2Bv)h4k*OYv8c4*M1C6yj=t3d|=$&uRBqHItCYy@J1Y zQy`?$QM(85Zp2>irjG)T#?ox6xZwdva1M0g9y~$8)`nZhKMStzrw7c%vVHeaOrFmcv?Iv!a``W zx8J(zpe_E&Fa~>VucwvDQLt#2vdz-zY@E!3jm*&$M9^6jSc-`xhf_dXTs}^qg$;np zK<7G9|dk1M&|_4#tf{AW;7;H*YUmO`8>Ra^XvxO zZ8VBp4TmB8tqqqOO@6}%{?|%~@csiePqD4bD`L|e`cy|*MVKTWNI3VQGK3>)$?=oq zGTNhcNV0`*4Ffw1$k{RYFBReZ3SJ$h)Kt9Zrbs_mZIfCX&`YmHQy+ak9J6yVYQ_PN zXMmErk6eZO>}-v)#!(a+xs>g*!GYOuPTbE}^0su|iA{{M z&!}v${7kWvt)}DJlr`dCAkC(l?rrAnmU^^6?N2Pkv7|)A1Nize$zgig57+QY-ei1sda;$^ zxHrUtB$mNDzMaFdIVZ)*L{0V@d0~e+u$EZ9t9~Vv zTz)V6e)UIguP*z>4#b-kCE9O4;7AE%tv(?&)xXPyW)rb?ERV&PCNRy5j>$vf>AAMg z%((Y(>Io&FJXkNvSg}jA`!&oq)k%cpEM^li#P(-*SZW&eR#!% z<}~0P=h&p@FuL1y-i5}xybF>g_+YwOw8!{+GpZ(8$KB7Z|Q=Tny=~yC?N$Q`-9@2W@~!zEiK{R6Y9{v0Q<$$EIJ>Z z#|1CeW>hMCrXgfUbL(917Rd@OMAX@8EdR_S+}~|DdU2LU{nafRTR?+5T$3vmYrO)y z{c(~b)mQ!t&vuP7G<0J|G2S=Ug^T_P=V?x*WV<{3=ferA^;A2W{U)gF_j9N-Z>kqh zv(WZ(!#f|y4%iMndDr(VZ49fnVyv&29+_d7v;R>v#H;5D=v zIr}04V>GuJLZ8va`_f?vtjx#Y6&LH|bD2hr$NBp-6#+uTA@{`)Rkqm}vkOGhq8)If z(`Yo_w)xgbv6g1$cKa~t9LbNOM?PHAa2!)OejhqdLx-FYBXZ0XKg@C8OC0r@f8bE% zOnpr=xUrg9z?AU~FX{9Gg0d>!E;AZ{>7HDw6^HRDZbC%lz-eR6`?zM8Aq6g1uPnBN z8&dfBxYAkilZ1$}ilm{;E@{UTR35XyD^$#YMcPhtV_QI2NDy9occ~KPB!yit!79yT zQI;02S4EX0*6vR|8jq#s+jeBs$-G}o``W!3DN1G(5kzUsbbZLdE_ql|FivUMR+)=BGJ3tqBr6M%N(8@-G; zuJ~>~YK!oKo3!G^wk>mt4#pC%ApbH%KYu{}i_-_=3I>3^C89{d$JfsudF8oT33K}Y zz`|h~$V>>R>pe5Q5Bc_}bCR`Ht?ozU%iOMpgJdS@AKUHyqDLtxIYg~B!JVq!es+rj zgD4u=&FDIjDb+7@@C~C!dl_VOk-AFs@qC@EfNkWh0MmiC-mtSxC{#p#7SCqZCYJ& zMr{sd2ZhJfBQe3hJsn$ zUck-^^jU^2!ww^6X!8!(xjm+Rfz78HIA;W8+ZWomqS2j!(s=6Bw>|p0QklxEefynu z%@Q>#hhVpLRYhUs`!K@~E@rc%HzAl~OGdTL?V&fVEtd}LJZ?Y6P7%Ie?;jkFB@(K% zeQqwU@;x_m|GJiATVT1r0JfYs@R~hrNhNu_&$?{1?;bXk33@SFADw))IN@&R?it_4 zQ*L|B=Wg(Aci`=(iaqg*fln^=3RA9zAI{CT1XInf!OP??Y^NT$IoP(_?>rn?o#U;<1=Hm|q_KOXEeyEvc7vFyzK205c zJ9F@6`|{{P!ASs{9RXr_S z$(#Bfe1g0h*jaOM9&Pk9{&BEYI`rJ0l@|AfJ%qzXmEkfqXch*s*VUa8Y6jI>(Y)Ce%z2qI;k>Gc`q7@rqJ*Z%>wc zYsY*2f{9VGa6LlNG0|N``CQr%@#_!g3HnKE+L7V;0!+z0@AveN_m|f#AyjPzL8}oR zRn6cV4qSic?bmKhawup`IbHPjibn4XCT@Zb+Kq_XNYgk3?N&|gl8X*ovmf3xr|I*K zJW6g$zwd<@Il49XBDPF%X=#NPVppWX=a#1P0hIb-nNU2;@1y@5j`Ec(H$8d)nHD2@ zSN-^93J+L3svE1bHAJc<-Ou)!1afKZw<^-0wUo+U(H-%t?ypM-?$XfdOhDw&9n zLOGp{geO3X$UzQC1-3HvZDS>pE3=BUPrt}kfcL92o<}G58A-Z0MDLJ?(@!BfXtobD z<DT^KN%Yk{`@W-5U!@#C!+e71qOM9;{Yen&bx}9X=N-Bq6n*L=8Uza1xSVO2PFwno zF*UF-_K-^)84c5P+FYz!XT+1T^Ud~*W{lT`l2Zv1pA-RBZ7T3U`I?DwAUFu!xoLHa z%bB5loP_a{Sr;D75~QGA;nhqtFupy?u9+{f8>x`ZXuLs~W9GXT@uWN2n&enrL48SJ zspI5krWCL^xL_2#uccz))5p@ML6Je0y#3PYoe*ms9mlihc5AIAE7_#?r

n4HJGS zO&B}R*yltwiv&Nz?dRK(sekcF`^8F2po)vBB3wrg_e+f zJ#UED+>5lzIV)&g-t3H1OT!R6X=gI~+)7^1uQuf7%kLs4QLUKlK1^)SHpbz2DWt>_V`!LBmKzDe?hj+x> z#DZ%KoD7q1X=M*sKlCU~#L}L)2a3Fxp&U}Wr;f5m@ov2<^yXW4pt^UW+Ur4;l?SiM zsIzrK>H4lyKt&~)3t!7bZ}6@z+s0oMJJNhVgIq?)wG-eBki>_mLY`$tY=l22ifpLV zMk-W;C$g!H)CWlk=1olj+X7e$@wqmHj^3$^So8(Pg87jMW$1=^yOVM$#IeN%&n@~V zT`!RI!dpgn?hWA|e3a3Ap=NWS$*0byQ67;SuR+DjM$T{hzSmlusVIXpF} z0{JlU);DDlFui_74$@imtTLvf*{CN`BOg6Z8_&dB_Q}k12dYjTvc4d5O0{ zKo~E69f8tCQ@{zDazpZ^s(E)YUkGJ6OxQ$vAS;S$5wI=d#i#jQAxC9hg}_l`gsYKE zFM4DJImYS`nl+q&>M(E+SQg5_#k@v&j|~k3216EL=b>?Qt86bjMTb4$n82=!WRMVC z%NG0ejJ21I>~T0{aMCBf0KkbHMnMsI2!96B=OM1+R@q0?n#abD5#8&JaFPQaO6_KS z^XO>#R(o9Z@xt*R{kOQ%$`A^45@ZLmW2-biDUR84Y59BN%fY{MG=%4VAMGH>)o8MX zzWjp*ZS%}&xh9G#$WNW^ek^LHg`2qOlT}flLQ01>U!@<3@Q1EsBP5GpELxtjyt&Un zIr#W6{t48^vpfUm%qgC*DnqUaQ}OMv%IgV~RBqSWZxcB6x#ptKyocN+{>n+>)!}cV z;3SgK@2O&4&}>y6o*oTD{SQN8{P*QKIAR=%s@1C7SZnPQJd6T2HdR!W{f2asq~KYl{fgiq_iE1kKLLQ_?1|R zvzmChy)*sV#xg6MZM#VaCVaK=0S!yT(}ds8WkcZYtU?Ws3^IQF zuk0y?PNK5X71OpPzX_0(2oD}QaoV?q8;wy#eTA+D96vS)8N=nBEt%wBW|2kDhH4^) zsAv_DN*ap$n-?g^HX#Pdn^i42+CDPfjl*>jpy=7ar%=!Y1(4bn8}fWDPcY*l<>quT zV1oMd1Ur#EB)|qx2%;kUH(Iix%;PKMqV7B!%&HWOIGuPTx__Sh>M%U7kX2wVXDOiq z?h>I#ZhQ|h$=Tn7+hG#H+k%tF_m$6ffWZ9n7W$JN;JA$$$-z6HhB?tIb6O;_E;l^hq4L{L*5eS7+=5B&OV6p z0bOk-0ZnUKW2ceSSZF+YUu)irmf|1SUuvXrK<62#Kd78m3*txR#^0f$OEpi!F(>Rw z2$vwYNB$g`mTiGQ;Dm&_6R7qDJQ?}Q`tq10kWTjFOO(Mc>Fm)M*zIIWQ;I599va*5 z6Qi_9Vkb2%(ywb4ZOk|N8CgsG)x*alR|&oajX#V~mYxYxw+P_TVXBm1tUAp@wMHs_ zD1%j#SG!4!)BZ)7+Z%!b*<1$9OSB>lX*wN5ZbmZmpXP5xsNu&PDh@cl;$aU(eL>f! zIq>Y|F2G0LXt$b~Kml|{jHhlhJHLw8qScDo zj^aZ3wI#g98<{3NQ;g_k!UEhI`cuWyd<;Y;NDZNWgqoXO2_Dw_FO87T9f*dU{28Mk z*NSkyAwS0oTE`?iz>dO#Hqti+OpiAjGE{^_yqR)B=e?#P=jI|L!L>9f&CL!KlXN+! z3IbQAb2L5~?FZ24SbkTZ9G5TGG^~03jqq!}2Z?6RTd_&Qqd?!2n73qg0OE|dW!BU9 zW2a~U&kloBL>`YOvE`^TAs-DV@lp?2w~wLW9kMfO;Yw#I8zf^U>(I`O+2x4oioB{M zHJbL^t7W#iI##}}`bQxqBhWI7qZBb53n2*m#*;zi=#8p0EG!p|{9euhVQ|%Vr7CkZ zv=bz1tB7&gZ2uY)qh_I)zS2*dNtkW1>)ewmp@lgJt-VDYFV8cE{56R*3)JGhtybv> zR&CViEGy4y??326d-lBd1GYR)SO$K-jkIaBf1h+MA%SlM5VtoWf|PhPi`-*`**{&U zi0C0=ugLKi?FhyfHIu_1#%pOOGYwd~mlX_VzSzyVQd-Ts9(?+^LXryJiZVtqGbG~G zDlh5ZFCaBA{J?TQPKJjrugnBvm$@XrULLoW0AavDrOGG-Tb#ZXU^-x99k5nVs-X%m zp?+zymf)>xQZs%eEni*jfk_+BT6?@pTbzOh)U~EH*=RuFkx;d8Ca$FR2UCTtGt$z8 z=mFURt8_}_gRLB@+Ma)1sAR2Ay6}K6eHvttQs?ryM#UuXqZ{+%R12l_N$dHnSm%bg zT{piqlUnZlrkxT3A~-3|G{$l$t%4@bhZ3*~rE1b_LI)=alNpgpf3{V~psSfUt4H_& zMfp?uS_>z>ow9tw+Ly zeX5u=Ed=KFM{aABO^h`YvD9&jO^24Av0XU31mM zT0F}N8swe@E7PQOSU=jbNITp>gXw(nB&N7k@9g3+et ze1OgK+1dNtTqIvuNe^u?V^ALeKHPil&&cgl>3sr6`hJbeL-K9bih&0=<>PxaBB?Rf zl;as4LZb=0myTJNA-y<9rL0g8|;O_>*X|>Y`6tXA-fBF`gO^q9ZBqDpg$O7ew-917%+N zju?wlpizTIvEM3U$JyyWm6$-Y+1<)&8*+D}Mx`N09_;ZYTF_V!p&ylwFb}zroSceP zYguEKOO$^oS-AUAXJD`ty#^k5oM9AAX?&&D^^Q(nFM!|+-G!QBU3wdle@Iwix*!Jm z5Jf@X0-L{f5_3Vcg&K}G&fX5a#Y_HLi|!!4Cv*$PyKzyry4qjr$g;)g%FCE=0zJ81 z=_{Ub4}o8h+QSZ0=m}3s_aR@b^4IhAM3>Y%*T-?$D3iGnRDZ^u`nLELf|2z%#`X{w zkTb}1_oCN7$^JD>G*J^MaVDM|1p@&;bB$6lB!6mbr1htRiI#S~ zBi=Nt*w+lm6OkDe@9~#u@W^pb;0)?L9253OWtff7bE6t7g4U60ba9W4pRfEV*dufe zEl;qzhtyd%U5}N1Z+eM+nF2NHC8dbPGn#Q#p=#W>n`9D}I041whMunkW~7L@VFL@f zI2rn>L6x|?EzYc`rS;m=t?=iIOxLFmuDRiOsnKK2q? zKeRQdSQL&Ibt_1TE)4P1I7g?KySg<)d`sX5fUu30eH?7wIGRe>_g^@uJ(setDO0}* z#T4(@z#R~D!&dV$QXF{?0BvvEGxpO6vLQZGk?@|SB_IEqh}xkSK{!+=Af+n;tW7X} zj>A%CFsR;oFJs0d+O7fbtio_Bz+3nl_x?UZ-^R)E^mVp2;

k;M+T~2ED$_Ps+67 z!!@!SbfgC)-(=vc5r8w%QoIY?ohQYa^D|cc9Zxa{C4xALDJ`wmHwKHZaOd?rt#cT` zdd_8tmkM0ADdCA$^}c&V@Ft4sPoWq{t!tF@FWn0J^A5B@^;<07s?OOa)$8&8>e+J@ zUKCf({73J(vCO}Dxo&cxROYSedJ#Dln2t5%RijQ(e;YYXX_pu>+vPE|C1*^ua$1B9=Lt5 z+-*=-tJt!34Q0ByLQ$+e*m#9_E-OW&WmD^a)P#Zk5DrfPzZ-c{jfj4g84c;tbg@Y( zT%0`rzyt9UudtZ$n3Yak`sRX%LvdNU#p#MnXww>0+ce_w0%tDRIIyMjTn79kBU@f5Lu;3>WQ27V`U={dW9e3dSX#&z^$9sGuZDazC}NBD~E6? zw}KFd=x-Vy^oO%9PIOK_DvG~R*G%_SF%Q&j{JyJT^-&x4Vg9t)vcWZG;m1Mde=aQk zdHy3eJ3G&xALeuL|6_kAC)0EK*&AQESiI6alb2+49Sdw_uZv37bdGX-N%?S{R(9-L8+_+ZR<|yX z+%7(hawhyx#9aCOx%;@=@nYnS9jYS#(Yk|o-6v*5P(!+Qc&CL3A9+$u_X+sz!@f_c zi(*Wp5Vr2l`@&=@8};jv{Fz2V?t7`m-+_|e!VP`O5naltqfumOsv)0-T8S{z981eK z&EBkH@=%~_NSDZZHsOpO8ApZ5c01;OdDtf_&2RqYXicdXeb(2H7f?`!>2&kWV|-X| z>+=sT`Ay9xtCw%y1wocbC@+p=dL^loDIUVdzYuKo3(jdFkFL7mZ=3t>^6`5$8V5Wf zgsAMj1kuvszLBkJx2{18n509!%H=pch%-uoH7qD*n5X~K8!SJg{TIfqAN?L?KNWbA z8H_A=iRdRO6rzny*4~~82lvbm({83zUv>p2PzLLXzf6rUZ@2GMKU}$w{}DI;GB1yk zK+h*PMYBz2$QLUiz5?47z_FBx2z8UZIq>!_tb-?LbD7FX7 zbIP)%$DMTtUzBQgWkoqgTYLL)MTNxmWKn^erpU6*s0!6q#X0w-!ZNQfX$zbB%DJnfZ&AXKbk?d=ldsI9x18`#Yjn_% z5YEw}U-mY~`M=;VjLzYG^x5-KoIILNcz6tfLU%ml11-V``U_pk$I&Q5h+Alx+eD|h z*KqzV2KNP;^U&FIdPKi{@=#@8f=DnbpU1e3MW2sL$d<{?eal|iE3(~5Pr26g>I~Ew zZ$sm^=Bn!8pGQnIwsUw%^Tq;=nNeZ8JfzG&IDE!()kgOz;7iLoE|(ZSXK^jnp#;gM z6zX{Ipmg5+CV#w)AjNX!bgtl`=G${ebY=mukIy8IG8*9D7>rCJ#CDv~j-$}iIaTov zRc`0YrW(o!i(U%XmX%bWf9(CHXJVWRafx)Rg`db?9}AvJmF1`eUS-S_gE|(Pb$6{9 zgFlw9>oBMWMRE74hLIN*fD*dPrkO3M1m1X4jDK^Cqb;H*qFS7%hi~Wb4xV`adIzc8 z!=E6e_cLP}2@)REzHWN_PKpLj)q7?l3hjIR=MC;CYVb?V}Yy!1L&G1L1x?lN6aii&g<2xl@ z2~ZwjF%(f;%h_m_Kbc+DJz>#!!$^3?OnR=SWT=~2=Z(O31hB;SPy(H;_;f2~=m-rR zV<{dVSX&G?BE?vhMjx&I83(>S+jY`Vt>k9|cQUmPtj7xlfp1tVh|rjo&55wz?jaqL z&pH5Sd%*d4D5P}QYG@G{DUuyItVaP~Sn*t+evdPrG*2y>#E5$x)raY7JwE6(K!{Q9 zs$Jd?P(C&MnWByjsxidhE$<#l6g2&#IH%Go0|$xgeOd_c#5?%F{Mog)D#Ht8;1%dq zrPk@aP)x32QDHxSBMh#~a^&Hnn|5yNNsiKB;n^|oXO+JjHZR0+>Ywmp|XCGIIJelT)Q_x zDa29R1qnRIjT}ZUZo2FaGZFIQ>$9lt0@IUvRv*9zKn9t!%jL(f%{g~8qk!&pT16JvhhpOAbB?fQ|UWthm#>P?pZYP&;TpBQ;-P5G;_qe;Sx zA1x1~9Cz9eW3|G~$_fhT?4;dB%k*zh4u1>TRZd!Ms&MQc3V&&+6g&tBMJ4Q^KiDw?)>1!5e!oqq<_cVn~^O z-pK;`=dJq~@#;~}Uo0%tw)t-j7ij8VT<8|6A=x}}H>3_J0t znpv zw)uWfWF+-!Rq*jLUd!3!Ln=*qkN|@JHx`TXh|M1?w)+0}22r{BM)r^7v00O>KQTQG zxJlc>b5+-SU7r>IJoG!5aV9vFCN$DD1ml5UP;g|-j&;FIl`&N?yh?h?fudG~6Na+Z zTki<5ms*BTjW(;U7gUR$Ok&u&BD@-bZ%i?=%d)SNzMrw9$kn>=s$g%S-KYNvae2F1 zcqm#ECLaWmJUYC*1gn+RBzD3%9sCHGlY)`LL#w-1kaa~8R$}xrBZ%B4Nk#ak4jze} zOk2l{Xrg2t7oA>I&&f`--VTX;uBkK7`@=_DE#i+crR|2&;d8uB#VdUDeW3$q?QDGK zC}pBAznWr9TZ^_;VpX5P4)?=RYsjwKBeBzh>%7{`nw)r?vQ+zZi8yoUO}YXP(~}-^ z8F-V?xiJyW#4!O4d#TaBv4s9l$ZUd1v7YIsA$4*DRJ|)|^0DN6Yo_MU_87;CjF1%X zg;UeXm_EhaE`9OQ+u^OkPh;KMu>2L-r;8fZk2;8MJ(tz;xz{il`yuTOJ~Q zcq=vzf4azLH1PnmNVmN2o|#FiDC=8Z(-F>ZH5Ofv<1NbyM|vcyJ`~L?b~qn-*7FTN zOgvq9vqdCCjK^!|N|2EXu#e}6t@GiP{?RD|z3g@EQW`mH>Au=f^h^o*Uj6OiQ@NvUnlVM`I_K=-EFV^Dno0SU4gy1TA=9D2R z5gZ`fvzV#yxnzA1p7uVU-NKE<!>J=%ttZPzpC}DR719!#eoTxKP%hm6(q2r%B$QEgMQy$m z-_IKRPD?{;qRi4a!X8clU)=vuW?e)aAfRTZX4f9Ca0trJqr;(z$=Y_G(ec4z@Iv7) zCE^$lUH4OjA7j|wGDGj6VBcZl3m@gqRT|6TitY)( zB2S^}m9d%c{wtESfpNe5LTG>9sXwdYAslpA12{b5;PbG|?K`AcdOPcOgxaB)BxZLCc{ zn~KNkMS7*|-eBRe6d4u*p*bzNXegrkI5TA&9eu&FzsdHd)nM{@9ltW>@XQH)-63W` z5x4WVYQR~&N+nYd^TAX|jOGD_{)fjOF|TG~x@@eG7VQQPS{HLKTgKzpv~I%k{FXEu zck7IGjAmT<+^BziOxe+y4)QL(Usbm2?K7Msl5_3zj*ihU>{IwRzQc3NoP+Gk71ne- z0GYj8=k2?LfYP@z$Y%+UO4>!lwc0Qc$tv)Wo3y}S12XO1yx4zEL$?@kj@y{ff`)hhS zKtZZtb0|m&3b6!>+j{`@Z)c1=u&V@k`8f;$Aa#2e*xi*tHW-x?)XoKXcb@<(8R~5B z>SPLahA~lfvNzR$y6A&omtq1znotiHSUG8X7%}_Z>=OXH^AaQhdxn*F2HtLdQNK-s zVZ8;${>w*+|J!hW%Ki(^-JJ=4!}I6N`u{=o+|JkUWdDNa_Heyl@o?PM;lFz8|Ku$~ES3I1*T*C3c z++Ljj#-8t9?tjDc(|lmQ33GSgoxeK+K?;@-XP`dd?#=&)7peo;@Ajbq?rcsS!0|^O z$KBzMfIG$10fsOsi>=x?(o$St-{v`X4@ZL4|pC$f-s`*7daj-M=wnqkOTG~RLS>#}CVrlwE?|+*N z07==wE+@9MGY4|sSpZ1G)x_oY-EeC$?0?nvUGBdi+!^>8v^!9(1+DM zFUZl=-USLVv9W=gx%_-^vNX4FG5m|)T_gRA-%pIPP#ZU>i=`>}SAhQopkM0xXZb&L zRJXSU+x^PnPpCl7KWyp0;CPq&pQy<5hgNqg{Z#Ou+QNK_A7pCp;Q71mzw}$bK^M2T zf&7a6ze4&8hkso=4Xb!A5s`o06eRT6!vSn+4RryUK+P@fgz1i&+v$Lo5MeqUZY4G) z2MMT!rL30|RKrVI)6~n_lpjL(SOg&S7z}$aLs>k)4lHJO_g6auq0XklR1)fn(p0wu z9s(Y=4z@54271`q*f|S$2-Dpa5rD;SlRo22C}pA(7{+bLCgfyVWWpX7{Z=}=`37a90WihcXxMIcTQG&Cvy-x zKR-W+jRVBN!2+XTarU%x0ei66IX{BQy8~kRBpc2+jfpMoS| z>Nr_CxLDfT0dEC^P3&D=gz0#p>^vO&CVVVVK7KA1c6JDag$?#&;Wss7V~22hK3&OvNf@Q#Dal)zq788Jtj3BP2P$ywJ8C6A~jD@|k%O9eFe5~ADbfAA}{yzW+ z)>{8wcUZ0b>&}HhzZVb!-B$B2TK38Q9gvg-*nPla~kf z$7=wyV(4w)dK)<01~#_=@HTM!8M}LA|M~b=jOAyV)NM@p=cCze%}2Zw&1&3D&p3 zY!39Z=Mn(`|A&76bBzA#u3+eiDK8Fyx_{gdR5*!lMH&T#MO?yUEYj_o#s^)}%D8RLXO z5p#91_zlUOz~2oNMydeyguU4y&h)T2?8R*bHFcq95E1#;Y}dmHY6f5fa`FRUGXC`g z@^Ev*#&%}FKhikaxcFd0x7#1k?$0z1zJI2%aq|97$HC9Z3manptt{+;9X3e*dm0z- z?a2FYX`Jlrx5vHzmc|J~@mm@jAK&l%IN7+lVdL$;mF3{){&#Y@f0xU~&&Bn78Yerf z5`N>y$;QJ5GymUY@NoZL2QYp-zo)Tt!&>yWvM_!;unFC7{P@^@uLD@w+hx?h)3I^i z&XN9>#>vL}4;dW%{9Lfn=HKWz*tog=A(xx?b|U#VIvBryQwMDP|E4VW?^rnaIB#cd zzw_h!9SbKL$L*}*Z*;J-T)$(1*%kNi^}@-^c=%iv_^IAkenn?Ux-gioKISe zLxPJ-oRgi8TS8nCHag@M=aFQW;N{^LmlFBEQyAVZLOHvDom}pmk`wmh1<=q)sYnC< E56rUt9RL6T literal 0 HcmV?d00001 From a9cdf5b11ff8b7aa45c4793bcddccbbc9a21772f Mon Sep 17 00:00:00 2001 From: Jean-Marc Petit Date: Wed, 25 Feb 2026 14:56:07 +0100 Subject: [PATCH 2/2] Added Commons disrectory and a few lines to .gitignore. 2026-02-25T13:65 UT. --- .gitignore | 12 +- Commons/Driver.f95 | 181 ++++++++++++ Commons/Makefile | 70 +++++ Commons/datadec.f95 | 36 +++ Commons/elemutils.f95 | 437 +++++++++++++++++++++++++++ Commons/ioutils.f95 | 141 +++++++++ Commons/kind_map | 18 ++ Commons/modelutils.f95 | 654 +++++++++++++++++++++++++++++++++++++++++ Commons/rot.f95 | 277 +++++++++++++++++ 9 files changed, 1825 insertions(+), 1 deletion(-) create mode 100644 Commons/Driver.f95 create mode 100644 Commons/Makefile create mode 100644 Commons/datadec.f95 create mode 100644 Commons/elemutils.f95 create mode 100644 Commons/ioutils.f95 create mode 100755 Commons/kind_map create mode 100644 Commons/modelutils.f95 create mode 100644 Commons/rot.f95 diff --git a/.gitignore b/.gitignore index e911088..9e36487 100644 --- a/.gitignore +++ b/.gitignore @@ -5,4 +5,14 @@ ModelUsed.dat OSSOS-model.out *.log *.ipynb -*.so \ No newline at end of file +*.so +*.o +*.mod +Classical/Classical +Commons/Driver +Commons/GiMeObj.f95 +Detached/Detached +Inner/Inner +Plutinos/Plutinos +Scattering/Scattering +Twotinos/Twotinos \ No newline at end of file diff --git a/Commons/Driver.f95 b/Commons/Driver.f95 new file mode 100644 index 0000000..495c15d --- /dev/null +++ b/Commons/Driver.f95 @@ -0,0 +1,181 @@ +program Driver +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* +! +! File Driver.f95 +! +! J.-M. Petit Observatoire de Besac +! Version 1: May 2013 +! +! The purpose of this program is to compare models of the Kuiper Belt to +! reality by actually comparing what the surveys (say CFEPS or OSSOS) +! would have found if the model was a good representation of the real +! world, to the objects that were actually found. +! +! The workflow of the survey simulator driver is as follows: +! +! Loop (some condition): +! call GiMeObj(arg_list_1) +! Check model ended: +! set exit condition +! call Detos1(arg_list_2) +! Check detection and tracking: +! store results +! +! The GiMeObj routine is in charge of providing a single new object at +! each call. The Detos1 routine determines is the proposed object would +! have been detected by the survey. +! +! The survey simulator expects orbital elements with respect to ecliptic +! reference frame. +! +! Logical unit numbers 7 to 19 are reserved to use by Driver.f and +! SurveySubs.f and should not be used by GiMeObj or any other routine +! that you may add to the driver. Please use logical unit numbers +! starting from 20. +! +! As currently written, the driver includes a file 'GiMeObj.f' +! containing the definition of the model. The correct way to use this +! feature is to have one's GiMeObj routine in a file and +! create a symobolic link: +! +! ln -s GiMeObj.f +! +! For more information, please refer to README files in the source directory. +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* + + use gimeobjut + + implicit none + + integer, parameter :: screen = 6, keybd = 5, verbose = 9 + integer :: lun_h + type(t_orb_m) :: o_m +! color array NEEDS to be length 10 or more! + real (kind=8) :: h, epoch, hmax, rn_iter + integer :: ierr, seed, n_iter, nmax, nchar, values(8), seed_mod + character(80) :: det_outfile, line + character(100) :: comments + character(10) :: time + character(8) :: date + character(5) :: zone + logical :: keep_going, rec + + lun_h = 10 + keep_going = .true. + +! Get arguments +! Seed for random number generator + read (keybd, *, err=9999) seed + seed_mod = seed +! -maximum number of trials (<0) or calibrated number of objects + read (keybd, *, err=9999) nmax +! Maximum value of H to draw + read (keybd, *, err=9999) hmax +! Do we want to record the objects inside GiMeObj ? + read (keybd, *, err=9999) rec +! Name for the model objects outfile + read (keybd, '(a)', err=9999) det_outfile + det_outfile = strip_comment(det_outfile) + +! print *, n_track_max + +! Open output files and write header + open (unit=lun_h, file=det_outfile, status='new', err=9500) + + call date_and_time(date, time, zone, values) + write (lun_h, '(a17,a23,2x,a5)') '# Creation time: ', & + date(1:4)//'-'//date(5:6)//'-'//date(7:8)//'T' & + //time(1:2)//':'//time(3:4)//':'//time(5:10), zone + write (lun_h, '(''#'')') + write (lun_h, '(''# Seed: '', i10)') seed + write (lun_h, '(''#'')') + write (lun_h, '(a,a)') & + '# a e i Omega omega M', & + ' H epoch comment ' + +! Initialize counters + n_iter = 0 + rn_iter = 0.d0 + +! Open and read in object distribution +100 continue + +! Main loop: loop on objects + if (keep_going) then + +! Select object. + nchar = 0 + call GiMeObj (seed_mod, nmax, hmax, rec, & + o_m, epoch, h, comments, nchar, ierr) + + if (ierr .eq. -10) then ! Something wrong with this object, go to next one + goto 100 + else if (ierr .eq. -20) then ! Something very wrong happend, stop + write (screen, '(a)') 'GiMeObj returned -20, stopping.' + goto 2010 + else if (ierr .eq. 100) then ! Reached end of model, prepare to stop + keep_going = .false. + end if + +! Count number of iterations; may exceed 2**31, so use a real*8 helper + n_iter = n_iter + 1 + if (n_iter .gt. 2000000000) then + rn_iter = rn_iter + dble(n_iter) + n_iter = 0 + if (rn_iter .ge. 9.d9) goto 2000 + end if + + write (lun_h, 9000) o_m%a, o_m%e, o_m%inc/drad, o_m%node/drad, & + o_m%peri/drad, o_m%m/drad, h, epoch, comments(1:nchar) + + goto 100 +! end of the if( keep_going ) loop + end if + +2000 continue + write (lun_h, '(''#'')') + write (lun_h, '(''# Total number of objects: '', f11.0)') & + rn_iter + dble(n_iter) + close (lun_h) + + call exit (0) + +9000 format (f9.4,1x,5(f8.4,1x),f6.2,1x,f13.5,1x,a) + +9500 continue + write (screen, *) 'File --- "', det_outfile, '" already exists. ' + goto 9502 + +9502 continue + write (screen, *) 'Make sure "', det_outfile, & + '" do not exist and restart ModelDriver.' + call exit(-1) + +9999 continue + write (screen, *) 'Usage: ModelDriver < input' + write (screen, *) + write (screen, *) 'This will read the following from the keyboard:' + write (screen, *) '' + write (screen, *) '' + write (screen, *) '' + write (screen, *) 'rec' + write (screen, *) '' + write (screen, *) 'where:' + write (screen, *) & + ': integer used as seed for the random number generator' + write (screen, *) ': -maximum number of trials if < 0' + write (screen, *) ' Calibrated number of object if >= 0' + write (screen, *) ': maximum value of H up to which to draw' + write (screen, *) & + ': do we want to record the model objects inside GiMeObj' + write (screen, *) ': output file for model objects' + + call exit(0) + + 2010 continue + write (lun_h,*) "# Failed while generating model object." + call exit(-20) + + +end program Driver diff --git a/Commons/Makefile b/Commons/Makefile new file mode 100644 index 0000000..a0722ca --- /dev/null +++ b/Commons/Makefile @@ -0,0 +1,70 @@ +U = datadec ioutils elemutils rot +OBJU = $(addsuffix .o,${U}) +SRCU = $(addsuffix .f95,${U}) +FPPU = $(addsuffix .fpp,${U}) +MODU = $(addsuffix .mod,${U}) + +M = modelutils GiMeObj +OBJM = $(addsuffix .o,${M}) +SRCM = $(addsuffix .f95,${M}) +FPPM = $(addsuffix .fpp,${M}) +MODM = $(addsuffix .mod,${M}) + + +GOBJ = GiMeObjF95 + +FC = gfortran +FPP = gfortran -E -x f95-cpp-input -fPIC -DPYTHON + +FFLAGS = -O3 -fPIC -x f95-cpp-input +FFLAGP = -O3 +CFLAGS = +CPPFLAGS = + +TOPTARGETS := clean test +GIMEOBJ ?= OSSOS-classical-model + +all: Driver + +%.fpp : %.f95 + $(FPP) $< -o $@ + +%.o : %.f95 + $(FC) $(FFLAGS) -c $< + +.PHONY: clean + +clean: rmtmp + rm -f _$(GOBJ)*.so Driver GiMeObj.f95 + +rmtmp: + rm -f *~ *.o *.mod *.fpp f90wrap*f90 + rm -f $(GOBJ).py + rm -f *.log LOG .f2py_f2cmap + +.PHONY: link + +GiMeObj.f95: link + +link: + \rm -f GiMeObj.f95 + \rm -f GiMeObj.o + ln -s $(GIMEOBJ).f95 GiMeObj.f95 + +Driver: link Driver.f95 $(OBJU) ${SRCM} Makefile + $(FC) $(FFLAGP) -o Driver $(OBJU) ${SRCM} Driver.f95 + +Modules: _$(GOBJ).so Makefile + echo "Modules have been built" + +_$(GOBJ).so: link $(FPPU) $(OBJU) $(GIMEOBJ).f95 Makefile + \rm -f _$(GOBJ).so + f90wrap -k kind_map -m $(GOBJ) $(FPPU) 1>f90wrap_GiMeObj.log 2>&1 + f2py-f90wrap --fcompiler=$(FC) -c -m _$(GOBJ) $(OBJU) f90wrap_*.f90 1>f2py_GiMeObj.log 2>&1 + \rm f90wrap_*.f90 + \rm *.fpp *.mod *.o + sed -e 's/import _$(GOBJ)/from . import _$(GOBJ)/' -i.bck $(GOBJ).py + rm $(GOBJ).py.bck + rm .f2py_f2cmap + mv _$(GOBJ)*.so _$(GOBJ).so + diff --git a/Commons/datadec.f95 b/Commons/datadec.f95 new file mode 100644 index 0000000..a8c0c32 --- /dev/null +++ b/Commons/datadec.f95 @@ -0,0 +1,36 @@ +module datadec + + ! define length of array parameters + integer, parameter :: nw_max = 20, nparmax = 20 + + ! define some useful constants + real (kind=8), parameter :: Pi = 3.141592653589793238d0, drad = Pi/180.0D0, & + TwoHours = 2.d0/24.d0, TwoPi = 2.0d0*Pi, eps = 1.d-14, & + Flatten = 0.003352813d0, & ! flattening of earth, 1/298.257 + Equat_Rad = 6378137.0d0 ! equatorial radius of earth, meters + + real (kind=8), parameter :: km2AU = 149597870.700d0 + + + real (kind=8), parameter :: gmb = 1.d0+1.d0/6023600.0d0+1.d0/408523.71d0 & + +1.d0/328900.56d0+1.d0/3098708.0d0+1.d0/1047.3486d0+1.d0/3497.898d0 & + +1.d0/22902.98d0+1.d0/19412.24d0+1.d0/1.35d8 + + ! Internal variables + real (kind=8) :: om_lim_low, om_lim_high + common /om_lim_com/ om_lim_low, om_lim_high + + ! define data type to represent survey efficiency and pointings, and objects + type t_orb_m + real (kind=8) :: a, e, inc, node, peri, m + end type t_orb_m + + type t_orb_p + real (kind=8) :: a, e, inc, node, peri, tperi + end type t_orb_p + + type t_v3d + real (kind=8) :: x, y, z + end type t_v3d + +end module datadec diff --git a/Commons/elemutils.f95 b/Commons/elemutils.f95 new file mode 100644 index 0000000..1cbb0f7 --- /dev/null +++ b/Commons/elemutils.f95 @@ -0,0 +1,437 @@ +module elemutils + + use datadec + +contains +! \subroutine{coord\_cart} + + subroutine coord_cart (mu, o_m, p, v) + +! This routine transforms delaunay variables into cartisian +! variables. +! +! ANGLES ARE GIVEN IN RADIAN !!!! +! +! \subsection{Arguments} +! \subsubsection{Definitions} +! \begin{verse} +! \verb|a, e, inc, node, peri, m| = osculating elements \\ +! \verb|x, y, z, vx, vy, vz| = cartesian variables: $X, Y, Z, Px, Py, Pz$ +! \end{verse} +! +! Angles are in [radian] +! +! \subsubsection{Declarations} +! +!f2py intent(in) mu +!f2py intent(in) o_m +!f2py intent(out) p +!f2py intent(out) v + implicit none + + type(t_orb_m), intent(in) :: o_m + type(t_v3d), intent(out) :: p, v + real (kind=8), intent(in) :: mu + +! \subsection{Variables} +! \subsubsection{Definitions} +! \begin{verse} +! \verb|cos_e, sin_e, cos_i, sin_i| = sinus and cosines of $E$ and +! $i$ \\ +! \verb|delau| = Delaunay variables: $l, \cos(g), \sin(g), \cos(h), +! \sin(h), L, G, H$ \\ +! \verb|e| = eccentric anomaly \\ +! \verb|mat| = rotation matrix \\ +! \verb|q_vec, qp| = $q$ and $\dot q$ \\ +! \verb|tmp| = temporary variable +! \end{verse} +! +! \subsubsection{Declarations} + integer :: i + real (kind=8) :: delau(8), cos_e, cos_i, e, mat(3,3), q_vec(2), qp(2), & + sin_e, sin_i, tmp, signe, f, de, fp, fpp, fppp + +! Computation of sinus and cosines of angles. + + signe = 1.d0 + if (o_m%a .lt. 0.d0) signe = -1.d0 + cos_i = dcos(o_m%inc) + sin_i = dsqrt(1.d0 - cos_i**2) + delau(2) = dcos(o_m%peri) + delau(3) = dsin(o_m%peri) + delau(4) = dcos(o_m%node) + delau(5) = dsin(o_m%node) + delau(1) = o_m%m - int(o_m%m/TwoPi)*TwoPi + delau(6) = signe*dsqrt(mu*o_m%a*signe) + delau(7) = abs(delau(6))*dsqrt((1.d0 - o_m%e**2)*signe) + +! Rotation matrix. +! The rotation matrix is the composition of 3 matrices: $R_{xq} = +! R_3(-h) \cdot R_1(-i) \cdot R_3(-g)$: +! \begin{displaymath} +! R_{xq} = \left(\matrix{ +! \cos(h)\cos(g)-\frac{H}{G}\sin(h)\sin(g)& +! -\cos(h)\sin(g)-\frac{H}{G}\sin(h)\cos(g)& +! \sqrt{1-\frac{H^2}{G^2}}\sin(h) \cr +! \sin(h)\cos(g)+\frac{H}{G}\cos(h)\sin(g)& +! -\sin(h)\sin(g)+\frac{H}{G}\cos(h)\cos(g)& +! -\sqrt{1-\frac{H^2}{G^2}}\cos(h) \cr +! \sqrt{1-\frac{H^2}{G^2}}\sin(g)& +! \sqrt{1-\frac{H^2}{G^2}}\cos(g)& +! \frac{H}{G} \cr}\right); +! \end{displaymath} + + mat(1,1) = delau(4)*delau(2) - cos_i*delau(5)*delau(3) + mat(1,2) = -delau(4)*delau(3) - cos_i*delau(5)*delau(2) + mat(2,1) = delau(5)*delau(2) + cos_i*delau(4)*delau(3) + mat(2,2) = -delau(5)*delau(3) + cos_i*delau(4)*delau(2) + mat(3,1) = sin_i*delau(3) + mat(3,2) = sin_i*delau(2) + +! Eccentric anomaly. +! We solve iteratively the equation: +! \begin{displaymath} +! E - e \sin(E) = l +! \end{displaymath} +! using the accelerated Newton's method (see Danby). + + e = delau(1) + sign(.85d0, dsin(delau(1)))*o_m%e + i = 0 +1000 continue + sin_e = o_m%e*dsin(e) + f = e - sin_e - delau(1) + if (dabs(f) .gt. 1.d-14) then + cos_e = o_m%e*dcos(e) + fp = 1.d0 - cos_e + fpp = sin_e + fppp = cos_e + de = -f/fp + de = -f/(fp + de*fpp/2.d0) + de = -f/(fp + de*fpp/2.d0 + de*de*fppp/6.d0) + e = e + de + i = i + 1 + if (i .lt. 20) goto 1000 + write (6, *) 'COORD_CART: No convergence: ', i, f + write (6, *) mu, e, de + write (6, *) o_m%a, o_m%e, o_m%inc + write (6, *) o_m%node, o_m%peri, o_m%m + stop + end if +1100 continue + +! Coordinates relative to the orbit. +! The cartisian coordinate are given by $\vec X = R_{xq} \vec q$ +! and $\vec P = R_{xq} \dot{\vec q}$, where: +! \begin{eqnarray*} +! \vec q & = & \left(\frac{L^2}{\mu}(\cos(E) - e), +! \frac{GL}{\mu}\sin(E), 0\right), \\ +! \dot{\vec q} & = & \frac{\mu}{L(1 - e\cos(E))} +! \left(-\sin(E), \frac{G}{L}\cos(E), 0\right) +! \end{eqnarray*} + + cos_e = dcos(e) + sin_e = dsin(e) + q_vec(1) = delau(6)**2*(cos_e - o_m%e)/mu + q_vec(2) = delau(7)*delau(6)*sin_e/mu + tmp = mu/(delau(6)*(1.d0 - o_m%e*cos_e)) + qp(1) = -sin_e*tmp + qp(2) = delau(7)*cos_e*tmp/delau(6) + +! Cartisian coordinates + + p%x = mat(1,1)*q_vec(1) + mat(1,2)*q_vec(2) + p%y = mat(2,1)*q_vec(1) + mat(2,2)*q_vec(2) + p%z = mat(3,1)*q_vec(1) + mat(3,2)*q_vec(2) + v%x = mat(1,1)*qp(1) + mat(1,2)*qp(2) + v%y = mat(2,1)*qp(1) + mat(2,2)*qp(2) + v%z = mat(3,1)*qp(1) + mat(3,2)*qp(2) + + return + end subroutine coord_cart + +! \subroutine{pos\_cart} + + subroutine pos_cart (o_m, p) + +! This routine transforms delaunay variables into cartisian +! variables, positions only. +! +! \subsection{Arguments} +! \subsubsection{Definitions} +! \begin{verse} +! \verb|a, e, inc, node, peri, m| = osculating elements \\ +! \verb|x, y, z| = cartesian variables: $X, Y, Z$ +! \end{verse} +! +! Angles are in [radian] +! +! \subsubsection{Declarations} +! +!f2py intent(in) o_m +!f2py intent(out) p + implicit none + + type(t_orb_m), intent(in) :: o_m + type(t_v3d), intent(out) :: p + +! \subsection{Variables} +! \subsubsection{Definitions} +! \begin{verse} +! \verb|cos_e, sin_e, cos_i, sin_i| = sinus and cosines of $E$ and +! $i$ \\ +! \verb|delau| = Delaunay variables: $l, \cos(g), \sin(g), \cos(h), +! \sin(h), L, G, H$ \\ +! \verb|e| = eccentric anomaly \\ +! \verb|mat| = rotation matrix \\ +! \verb|q_vec| = $q$ \\ +! \end{verse} +! +! \subsubsection{Declarations} + integer :: i + real (kind=8) :: delau(8), cos_e, cos_i, e, mat(3,3), q_vec(2), & + sin_e, sin_i, signe, f, de, fp, fpp, fppp + +! Computation of sinus and cosines of angles. + + signe = 1.d0 + if (o_m%a .lt. 0.d0) signe = -1.d0 + cos_i = dcos(o_m%inc) + sin_i = dsqrt(1.d0 - cos_i**2) + delau(2) = dcos(o_m%peri) + delau(3) = dsin(o_m%peri) + delau(4) = dcos(o_m%node) + delau(5) = dsin(o_m%node) + delau(1) = o_m%m - int(o_m%m/TwoPi)*TwoPi + delau(6) = signe*dsqrt(o_m%a*signe) + delau(7) = abs(delau(6))*dsqrt((1.d0 - o_m%e**2)*signe) + +! Rotation matrix. +! The rotation matrix is the composition of 3 matrices: $R_{xq} = +! R_3(-h) \cdot R_1(-i) \cdot R_3(-g)$: +! \begin{displaymath} +! R_{xq} = \left(\matrix{ +! \cos(h)\cos(g)-\frac{H}{G}\sin(h)\sin(g)& +! -\cos(h)\sin(g)-\frac{H}{G}\sin(h)\cos(g)& +! \sqrt{1-\frac{H^2}{G^2}}\sin(h) \cr +! \sin(h)\cos(g)+\frac{H}{G}\cos(h)\sin(g)& +! -\sin(h)\sin(g)+\frac{H}{G}\cos(h)\cos(g)& +! -\sqrt{1-\frac{H^2}{G^2}}\cos(h) \cr +! \sqrt{1-\frac{H^2}{G^2}}\sin(g)& +! \sqrt{1-\frac{H^2}{G^2}}\cos(g)& +! \frac{H}{G} \cr}\right); +! \end{displaymath} + + mat(1,1) = delau(4)*delau(2) - cos_i*delau(5)*delau(3) + mat(1,2) = -delau(4)*delau(3) - cos_i*delau(5)*delau(2) + mat(2,1) = delau(5)*delau(2) + cos_i*delau(4)*delau(3) + mat(2,2) = -delau(5)*delau(3) + cos_i*delau(4)*delau(2) + mat(3,1) = sin_i*delau(3) + mat(3,2) = sin_i*delau(2) + +! Eccentric anomaly. +! We solve iteratively the equation: +! \begin{displaymath} +! E - e \sin(E) = l +! \end{displaymath} +! using the accelerated Newton's method (see Danby). + + e = delau(1) + sign(.85d0, dsin(delau(1)))*o_m%e + i = 0 +1000 continue + sin_e = o_m%e*dsin(e) + f = e - sin_e - delau(1) + if (dabs(f) .gt. 1.d-14) then + cos_e = o_m%e*dcos(e) + fp = 1.d0 - cos_e + fpp = sin_e + fppp = cos_e + de = -f/fp + de = -f/(fp + de*fpp/2.d0) + de = -f/(fp + de*fpp/2.d0 + de*de*fppp/6.d0) + e = e + de + i = i + 1 + if (i .lt. 20) goto 1000 + write (6, *) 'POS_CART: No convergence: ', i, f + write (6, *) e, de + write (6, *) o_m%a, o_m%e, o_m%inc + write (6, *) o_m%node, o_m%peri, o_m%m + stop + end if +1100 continue + +! Coordinates relative to the orbit. +! The cartisian coordinate are given by $\vec X = R_{xq} \vec q$ +! and $\vec P = R_{xq} \dot{\vec q}$, where: +! \begin{eqnarray*} +! \vec q & = & \left(\frac{L^2}{\mu}(\cos(E) - e), +! \frac{GL}{\mu}\sin(E), 0\right), \\ +! \dot{\vec q} & = & \frac{\mu}{L(1 - e\cos(E))} +! \left(-\sin(E), \frac{G}{L}\cos(E), 0\right) +! \end{eqnarray*} + + cos_e = dcos(e) + sin_e = dsin(e) + q_vec(1) = delau(6)**2*(cos_e - o_m%e) + q_vec(2) = delau(7)*delau(6)*sin_e + +! Cartisian coordinates + + p%x = mat(1,1)*q_vec(1) + mat(1,2)*q_vec(2) + p%y = mat(2,1)*q_vec(1) + mat(2,2)*q_vec(2) + p%z = mat(3,1)*q_vec(1) + mat(3,2)*q_vec(2) + + return + end subroutine pos_cart + +! \subroutine{osc\_el} + + subroutine osc_el (mu, p, v, o_m) + +! This routine transforms cartisian variables into delaunay +! variables. +! +! \subsection{Arguments} +! \subsubsection{Definitions} +! \begin{verse} +! \verb|a, e, inc, node, peri, m| = osculating elements \\ +! \verb|cart| = cartesian variables: $X, Y, Z, Px, Py, Pz$ +! \end{verse} +! +! \subsubsection{Declarations} +!f2py intent(in) mu +!f2py intent(in) p +!f2py intent(in) v +!f2py intent(out) o_m + implicit none + + type(t_v3d), intent(in) :: p, v + type(t_orb_m), intent(out) :: o_m + real (kind=8), intent(in) :: mu + +! \subsection{Variables} +! \subsubsection{Definitions} +! \begin{verse} +! \verb|cos_i, sin_i| = sinus and cosines of $i$ \\ +! \verb|delau| = Delaunay variables: $l, \cos(g), \sin(g), \cos(h), +! \sin(h), L, G, H$ \\ +! \verb|e| = eccentric anomaly \\ +! \verb|f| = $f$ true anomaly \\ +! \verb|g| = $g$ argument of pericenter \\ +! \verb|h_vec| = $\vec h = \vec X \times \vec P$ \\ +! \verb|p_vec| = $\vec p = -\mu \frac{\vec X}{r} +! - \vec h \times \vec P$ \\ +! \verb|r| = radial distance \\ +! \verb|tmp1, tmp2| = temporary variables \\ +! \verb|v2| = velocity squared +! \end{verse} +! +! \subsubsection{Declarations} + real (kind=8) :: delau(8), e, f, h_vec(3), p_vec(3), & + r, tmp1, tmp2, v2, signe + +! Computation of angular momentum and eccentricity vector. +! \begin{eqnarray*} +! \vec h & = & \vec X \times \vec P, \\ +! \vec p & = & -\mu \frac{\vec X}{|\vec X|} +! - \vec h \times \vec P +! \end{eqnarray*} + + h_vec(1) = p%y*v%z - p%z*v%y + h_vec(2) = p%z*v%x - p%x*v%z + h_vec(3) = p%x*v%y - p%y*v%x + r = 1.d0/dsqrt(p%x**2 + p%y**2 + p%z**2) + p_vec(1) = -mu*p%x*r - h_vec(2)*v%z + h_vec(3)*v%y + p_vec(2) = -mu*p%y*r - h_vec(3)*v%x + h_vec(1)*v%z + p_vec(3) = -mu*p%z*r - h_vec(1)*v%y + h_vec(2)*v%x + +! Computation of momenta. +! \begin{eqnarray*} +! L & = & \mu\sqrt{\frac{1} +! {\frac{2\mu}{|\vec X|}-|\vec P|^2}}, \\ +! G & = & |\vec h|, \\ +! H & = & h_z +! \end{eqnarray*} + + v2 = v%x**2 + v%y**2 + v%z**2 + tmp1 = 2.d0*mu*r - v2 + signe = 1.d0 + if (tmp1 .lt. 0.d0) signe = -1.d0 + delau(6) = signe*mu/dsqrt(signe*tmp1) + delau(7) = h_vec(1)**2 + h_vec(2)**2 + h_vec(3)**2 + delau(8) = h_vec(3) + o_m%e = dsqrt(dmax1(mu**2 - delau(7)*tmp1, 0.d0))/mu + delau(7) = dsqrt(delau(7)) + + if ((p%z .eq. 0.d0) .and. (v%z .eq. 0.d0)) then + delau(4) = 1.d0 + delau(5) = 0.d0 + tmp1 = 1.d0/dsqrt(p_vec(1)**2 + p_vec(2)**2) + delau(2) = p_vec(1)*tmp1 + delau(3) = p_vec(2)*tmp1 + else + +! Longitude of node. +! \begin{eqnarray*} +! \cos(h) & = & -\frac{h_y}{\sqrt{h_x^2 + h_y^2}}, \\ +! \sin(h) & = & \frac{h_x}{\sqrt{h_x^2 + h_y^2}} +! \end{eqnarray*} + + tmp1 = 1.d0/dsqrt(h_vec(1)**2 + h_vec(2)**2) + delau(4) = -h_vec(2)*tmp1 + delau(5) = h_vec(1)*tmp1 + +! Argument of pericenter. +! Let us call $\vec N$ the vector derived from $\vec h$, pointing +! to the ascending node: +! \begin{displaymath} +! \vec N = \left(-h_y, h_x, 0\right) +! \end{displaymath} +! From this, we get: +! \begin{eqnarray*} +! \cos(g) & = & \frac{\vec N \cdot \vec p}{|\vec N||\vec p|}, \\ +! \sin(g) & = & \frac{\vec N \times \vec p}{|\vec N||\vec p|} +! \cdot \frac{\vec h}{|\vec h|} +! \end{eqnarray*} + + tmp2 = 1.d0/dsqrt(p_vec(1)**2 + p_vec(2)**2 + p_vec(3)**2) + delau(2) = (h_vec(1)*p_vec(2) - h_vec(2)*p_vec(1))*tmp1*tmp2 + delau(3) = ((h_vec(1)**2 + h_vec(2)**2)*p_vec(3) & + - h_vec(3)*(h_vec(1)*p_vec(1) + h_vec(2)*p_vec(2))) & + *tmp1*tmp2/delau(7) + end if + +! Mean anomaly +! We define $\vec X_{orb} = R_1(i) \cdot R_3(h) \vec X$. It turns +! out that $\vec X_{orb} = (r \cos(g+f), r \sin(g+f), 0)$. Hence: +! \begin{eqnarray*} +! \cos(g+f) & = & \cos(h) X + \sin(h) Y, \\ +! \sin(g+f) & = & \cos(i) \left(\cos(h) Y - \sin(h) X\right) +! + \sin(i) Z +! \end{eqnarray*} +! Furthermore, we have the relation: +! \begin{displaymath} +! \tan(\frac{E}{2}) = \sqrt{\frac{1 - e}{1 + e}} +! \tan(\frac{f}{2}) +! \end{displaymath} +! and finally: +! \begin{displaymath} +! l = E - e \sin(E) +! \end{displaymath} + + tmp1 = (p%x*p_vec(1) + p%y*p_vec(2) + p%z*p_vec(3)) + tmp2 = ((p%z*p_vec(2) - p%y*p_vec(3))*h_vec(1) & + + (p%x*p_vec(3) - p%z*p_vec(1))*h_vec(2) & + + (p%y*p_vec(1) - p%x*p_vec(2))*h_vec(3)) & + /delau(7) + f = datan2(tmp2, tmp1) + e = 2.d0*datan(dsqrt(signe*(1.d0 - o_m%e)/(1.d0 + o_m%e))*dtan(f/2.d0)) + o_m%m = e - o_m%e*dsin(e) + o_m%node = datan2(delau(5), delau(4)) + o_m%peri = datan2(delau(3), delau(2)) + o_m%a = signe*delau(6)**2/mu + o_m%inc = dacos(dmax1(dmin1(delau(8)/delau(7),1.d0),-1.d0)) + + return + end subroutine osc_el + +end module elemutils diff --git a/Commons/ioutils.f95 b/Commons/ioutils.f95 new file mode 100644 index 0000000..1d568e5 --- /dev/null +++ b/Commons/ioutils.f95 @@ -0,0 +1,141 @@ +module ioutils + + use datadec + +contains + + subroutine trim (base_name, i1, i2, finished, len) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine trims the input character string 'base_name' and returns +! the +! +! Author: Jean-Marc Petit +! version 1: December 2019 +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! base_name: Input string with file name in it (CH) +! len : Length of input string (I4) +! +! OUPUT +! i1 : Index of first character of file name (I4) +! i2 : Index of last character of file name (I4) +! finished: Whether it reached the end of the string when searching +! for first character +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) base_name +!f2py intent(out) i1 +!f2py intent(out) i2 +!f2py intent(out) finished +!f2py intent(in) len + implicit none + + integer, intent(in) :: len + integer, intent(out) :: i1, i2 + character(*), intent(in) :: base_name + logical, intent(out) :: finished + + finished = .false. + i1 = 1 +100 continue + if ((base_name(i1:i1) .eq. char(0)) & + .or. (base_name(i1:i1) .eq. char(9)) & + .or. (base_name(i1:i1) .eq. ' ')) then + i1 = i1 + 1 + if (i1 .eq. len) then + finished = .true. + return + end if + goto 100 + end if +101 continue + + i2 = len +110 continue + if ((base_name(i2:i2) .eq. char(0)) & + .or. (base_name(i2:i2) .eq. char(9)) & + .or. (base_name(i2:i2) .eq. ' ')) then + if (i2 .eq. i1) goto 111 + i2 = i2 - 1 + goto 110 + end if +111 continue + + return + end subroutine trim + + subroutine read_file_name (base_name, i1, i2, finished, len) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine trims the input character string 'base_name' and returns +! the +! +! Author: Jean-Marc Petit +! version 1: September 2017 +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! base_name: Input string with file name in it (CH) +! len : Length of input string (I4) +! +! OUPUT +! i1 : Index of first character of file name (I4) +! i2 : Index of last character of file name (I4) +! finished: Whether it reached the end of the string when searching +! for first character +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) base_name +!f2py intent(out) i1 +!f2py intent(out) i2 +!f2py intent(out) finished +!f2py intent(in) len + implicit none + + integer, intent(in) :: len + integer, intent(out) :: i1, i2 + character(*), intent(in) :: base_name + logical,intent(out) :: finished + + finished = .false. + i1 = 1 +100 continue + if ((base_name(i1:i1) .eq. char(0)) & + .or. (base_name(i1:i1) .eq. char(9)) & + .or. (base_name(i1:i1) .eq. ' ')) then + i1 = i1 + 1 + if (i1 .eq. len) then + finished = .true. + return + end if + goto 100 + end if +101 continue + + i2 = i1 + 1 +110 continue + if ((base_name(i2:i2) .ne. char(0)) & + .and. (base_name(i2:i2) .ne. char(9)) & + .and. (base_name(i2:i2) .ne. ' ')) then + if (i2 .eq. len) goto 111 + i2 = i2 + 1 + goto 110 + end if + i2 = i2 - 1 +111 continue + + return + end subroutine read_file_name + + character(100) function strip_comment(str) + character (*), intent(in) :: str + integer c_idx, i1, i2 + logical finished + strip_comment = str + c_idx = index(strip_comment, '!') + if (c_idx /= 0) then + strip_comment = strip_comment(1:c_idx-1) + end if + call trim(strip_comment, i1, i2, finished, len(strip_comment)) + strip_comment = strip_comment(i1:i2) + + end function strip_comment + +end module ioutils + diff --git a/Commons/kind_map b/Commons/kind_map new file mode 100755 index 0000000..1724b5b --- /dev/null +++ b/Commons/kind_map @@ -0,0 +1,18 @@ +{ + 'real': { '' : 'float', + '4' : 'float', + 'isp' : 'float', + '8' : 'double', + 'dp' : 'double', + 'idp' : 'double', + 'sng_ad' : 'float', + 'dbl_ad' : 'double'}, + 'complex' : { '' : 'complex_float', + '8' : 'complex_double', + '16' : 'complex_long_double', + 'dp' : 'complex_double'}, + 'integer' : { '2' : 'int', + '4' : 'int', + '8' : 'long_long', + 'dp' : 'long_long' }, +} diff --git a/Commons/modelutils.f95 b/Commons/modelutils.f95 new file mode 100644 index 0000000..2c756a7 --- /dev/null +++ b/Commons/modelutils.f95 @@ -0,0 +1,654 @@ +module modelutils + + use datadec + + interface + function func(nparam, param, inc) + real (kind=8) :: func + integer (kind=4), intent(in) :: nparam + real (kind=8), intent(in) :: param(*), inc + end function func + end interface + + type func_holder + procedure(func), pointer, nopass :: f_ptr => null() + end type func_holder + +contains + subroutine incdism (seed, nparam, param, incmin, incmax, inc, & + dist, ierr, func) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine draws randomly a variable according to probability +! density \verb|func| with parameters \verb|param|. Same as previous +! routine, but can remember up to 10 different distributions at a time. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! +! J-M. Petit Observatoire de Besancon +! Version 1 : October 2006 +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! seed : Random number generator seed (I4) +! nparam: Number of parameters (I4) +! param : Parameters (n*R8) +! incmin: Minimum inclination (R8) +! incmax: Maximum inclination (R8) +! dist : index of the selected distribution (I4) +! func : probability density function +! +! OUTPUT +! inc : Inclination (R8) +! ierr : Error code +! 0 : nominal run +! 10 : wrong input data +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in,out) seed +!f2py intent(in) nparam +!f2py intent(in), depend(nparam) :: param +!f2py intent(in) incmin +!f2py intent(in) incmax +!f2py intent(in) dist +!f2py intent(out) inc +!f2py intent(out) ierr + implicit none + + integer (kind=4), intent(in) :: nparam, dist + integer (kind=4), intent(inout) :: seed + integer (kind=4), intent(out) :: ierr + real (kind=8), intent(in) :: param(nparmax), incmin, incmax + real (kind=8), intent(out) :: inc + type(func_holder), intent(in) :: func + integer :: i, ilo, ihi, di + real (kind=8) :: random + integer (kind=4), parameter :: np = 16384, nd = 10 + real (kind=8), save :: proba(0:np,nd), inctab(0:np,nd) + logical, save :: first(nd) + + data first /.true.,.true.,.true.,.true.,.true., & + .true.,.true.,.true.,.true.,.true./ + + ierr = 0 + di = min(nd, dist) + if (first(di)) then + inctab(0,di) = incmin + proba(0,di) = 0.d0 + do i = 1, np + inctab(i,di) = incmin + dfloat(i)*(incmax-incmin)/dfloat(np) + proba(i,di) = func%f_ptr(nparam, param(1:nparam), inctab(i,di)) + proba(i-1,di) + end do + do i = 1, np + proba(i,di) = proba(i,di)/proba(np,di) + end do + first(di) = .false. + end if + + random = ran_3(seed) + inc = interp(proba(0,di), inctab(0,di), random, np+1) + + return + end subroutine incdism + + real (kind=8) function Variably_tapered(h, params) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine returns the number of objects brighter or equal to H +! following an exponentially tapered exponential with parameters in +! params. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! +! J-M. Petit Observatoire de Besancon +! Version 1 : October 2021 +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! seed : Random number generator seed (I4) +! params: parameters for the distribution (4*R8) +! +! OUTPUT +! Variably_tapered : Random value of H (R8) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) h +!f2py intent(in) params + implicit none + + real (kind=8), intent(in) :: params(4), h + + Variably_tapered = 10.d0**(params(3)*3.d0*(h-params(1))/5.d0) & + *exp(-10.d0**(-params(4)*3.d0*(h-params(2))/5.d0)) + + return + end function Variably_tapered + + real (kind=8) function Variably_tapered_diff(h, params) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine returns the differential number of objects at mag H per unit +! mag following an exponentially tapered exponential with parameters in +! params. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! +! J-M. Petit Observatoire de Besancon +! Version 1 : August 2023 +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! seed : Random number generator seed (I4) +! params: parameters for the distribution (4*R8) +! +! OUTPUT +! Variably_tapered_diff : Random value of H (R8) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) h +!f2py intent(in) params + implicit none + + real (kind=8), intent(in) :: params(4), h + + Variably_tapered_diff = (log(10.0d0)*3.0d0/5.0d0) & + *10.0d0**(params(3)*3.0d0*(h-params(1))/5.0d0) & + *(params(3) + params(4)*10.0d0**(-params(4)*3.0d0*(h-params(2))/5.0d0)) & + *exp(-10.0d0**(-params(4)*3.0d0*(h-params(2))/5.0d0)) + + return + end function Variably_tapered_diff + + real (kind=8) function offgau (nparam, param, inc) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine returns the unnormalized inclination "probability" +! density as a non-zero centered gaussian. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! +! J-M. Petit Observatoire de Besancon +! Version 1 : August 2014 +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! nparam: Number of parameters (I4) +! param : Parameters (n*R8) +! inc : Inclination [rad] (R8) +! +! OUPUT +! offgau: Value of the probability (R8) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) nparam +!f2py intent(in), depend(nparam) :: param +!f2py intent(in) inc + implicit none + + integer (kind=4), intent(in) :: nparam + real (kind=8), intent(in) :: param(*), inc + real (kind=8), parameter :: Pi = 3.141592653589793238d0, TwoPi = 2.0d0*Pi + real (kind=8) :: fe, s1, angle, t0, t1, t3 + + if (nparam .ne. 2) stop + t0 = param(1) + s1 = param(2) + t1 = 2.*s1**2 + angle = mod(inc, TwoPi) + if (angle .gt. Pi) angle = angle - TwoPi + t3 = -(t0-angle)**2 + if (t3 .lt. -300.d0*t1) then + fe = 0.d0 + else + fe = exp(t3/t1) + end if + offgau = dsin(angle)*fe + + return + end function offgau + + real (kind=8) function cold_inc_sq (nparam, param, inc) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine returns the unnormalized inclination "probability" +! density for cold objects. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! +! J-M. Petit Observatoire de Besancon +! Version 1 : April 2023, 7th +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! nparam: Number of parameters (I4) +! param : Parameters (n*R8) +! inc : Inclination [rad] (R8) +! +! OUPUT +! cold_inc_sq: Value of the probability (R8) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) nparam +!f2py intent(in), depend(nparam) :: param +!f2py intent(in) inc + implicit none + + integer (kind=4), intent(in) :: nparam + real (kind=8), intent(in) :: param(*), inc + real (kind=8), parameter :: Pi = 3.141592653589793238d0, & + TwoPi = 2.0d0*Pi, drad = Pi/180.0d0 + real (kind=8), save :: angle, a0, p + logical, save :: first + + data a0 /5.0d0/, p /1.0d0/ + data first /.true./ + + if (first) then + a0 = a0*drad + if (nparam .ge. 1) then + a0 = param(1) + end if + if (nparam .ge. 2) then + p = param(2) + end if + first = .false. + end if + + angle = mod(inc, TwoPi) + if ((angle .gt. 0.0d0) .and. (angle .lt. a0)) then + cold_inc_sq = angle**p*(a0**p - angle**p) + else + cold_inc_sq = 0.0d0 + end if + + return + end function cold_inc_sq + + real (kind=8) function warm_inc_sq (nparam, param, inc) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine returns the unnormalized inclination "probability" +! density for warm objects. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! +! J-M. Petit Observatoire de Besancon +! Version 1 : April 2023, 7th +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! nparam: Number of parameters (I4) +! param : Parameters (n*R8) +! inc : Inclination [rad] (R8) +! +! OUPUT +! warm_inc_sq: Value of the probability (R8) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) nparam +!f2py intent(in), depend(nparam) :: param +!f2py intent(in) inc + implicit none + + integer (kind=4), intent(in) :: nparam + real (kind=8), intent(in) :: param(*), inc + real (kind=8), parameter :: Pi = 3.141592653589793238d0, & + TwoPi = 2.0d0*Pi, drad = Pi/180.0d0 + real (kind=8), save :: angle, a0, p + logical, save :: first + + data a0 /12.0d0/, p /1.0d0/ + data first /.true./ + + if (first) then + a0 = a0*drad + if (nparam .ge. 1) then + a0 = param(1) + end if + if (nparam .ge. 2) then + p = param(2) + end if + first = .false. + end if + + angle = mod(inc, TwoPi) + if ((angle .gt. 0.0d0) .and. (angle .lt. a0)) then + warm_inc_sq = angle**p*(a0**p - angle**p) + else + warm_inc_sq = 0.0d0 + end if + + return + end function warm_inc_sq + + real (kind=8) function H_dist_cold_2(seed, nparam, hparam) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine draws randomly a number according to the cold belt H_r +! distribution, represented by an exponentially tapered exponential, +! with parameters I've fitted on the OSSOS cold belt data. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! +! J-M. Petit Observatoire de Besancon +! Version 1 : October 2021 +! Version 2 : December 2021- updated parameter values. +! Version 3 : January 2022 - forcing slope at small sizes. +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! seed : Random number generator seed (I4) +! nparam: Number of parameters (I4) +! hparam: Parameters for asymptotic slope(s) (n*R8) +! hparam(1): start of asymptote +! hparam(2): contrast at start of asymptote +! hparam(3): slope of asymptote +! hparam(4): end of asymptote +! +! OUTPUT +! H_dist_cold : Random value of H (R8) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in,out) seed + implicit none + + integer (kind=4), parameter :: np = 16384 + integer (kind=4), intent(inout) :: seed + integer (kind=4), intent(in) :: nparam + real (kind=8), intent(in) :: hparam(4) + integer (kind=4) :: i + real (kind=8) :: params(4), random, h_min, h_max, n, c + real (kind=8), save :: proba(0:np), htab(0:np) + logical, save :: first + + data first /.true./ + data params /-2.6d0, 8.1d0, 0.666d0, 0.42d0/ + data h_min /4.6d0/ + + if (first) then + h_max = hparam(nparam) + htab(0) = h_min + proba(0) = 1.d-10 + n = Variably_tapered(hparam(1), params) + c = hparam(2) + do i = 1, np + htab(i) = h_min + dfloat(i)*(h_max-h_min)/dfloat(np) + if (htab(i) .lt. hparam(1)) then + proba(i) = Variably_tapered(htab(i), params) + else + proba(i) = n & + + n*c*(10.0d0**(hparam(3)*(htab(i)-hparam(1))) - 1.0d0) + end if + end do + do i = 0, np + proba(i) = proba(i)/proba(np) + end do + first = .false. + end if + + random = ran_3(seed) + H_dist_cold_2 = interp(proba, htab, random, np+1) + + return + end function H_dist_cold_2 + + subroutine H_diff_hot_5(nparam, hparam, np, hs, dist) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine computes the differential distribution of the hot belt H_r, +! represented by an exponentially tapered exponential, with parameters +! fitted on the OSSOS cold belt data, then scaled to hot. +! This version provides a continuous differential function, except for the divot. +! +! This version has modified parameters to obtain a cumulative distribution +! that looks like, and has the same normalisation as H_dist_hot_3. +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! +! J-M. Petit Observatoire de Besancon +! Version 1 : June 2024 - From H_dist_hot_4 +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! nparam: Number of parameters (I4) +! hparam: Parameters for asymptotic slope(s) (n*R8) +! hparam(1): start of asymptote +! hparam(2): contrast at start of asymptote +! hparam(3): slope of asymptote +! hparam(4): end of asymptote +! np : Number of points to return +! hs : Values of H at which we want the distribution +! +! OUTPUT +! hs(0) : H_min +! dist : Values of the cumulative distribution +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) nparam +!f2py intent(in), depend(nparam) :: hparam +!f2py intent(in) np +!f2py intent(in, out) hs +!f2py intent(out) dist + implicit none + + integer (kind=4), intent(in) :: nparam, np + real (kind=8), intent(in) :: hparam(*) + real (kind=8), intent(inout) :: hs(0:np) + real (kind=8), intent(out) :: dist(0:np) + integer (kind=4) :: i + real (kind=8), save :: params(4), h_min, h_max + real (kind=8), save :: a1, a2, a3, sl1, sl2, sl3, h1, h2, h3, scale + + data params /-2.6d0, 8.1d0, 0.666d0, 0.42d0/ + data h_min /-1.d0/ + data sl1 /0.13d0/, sl2 /0.7d0/ + data h1 /2.5d0/, h2 /5.8d0/ + data scale /2.2d0/ + + h_max = hparam(nparam) + h3 = hparam(1) + a2 = scale*Variably_tapered_diff(h2, params) + a3 = hparam(2)*scale*Variably_tapered_diff(h3, params) + sl3 = hparam(3) + a1 = a2*10.0d0**(sl2*(h1-h2)) +! The normalisation is done with the exponentially tapered exponential, +! as fitted on the OSSOS cold component, then scaled by 2. This +! determines the normalisation of the exponential between H = h2 +! and H = h1 +! N( from reference frame to ecliptic ; +! -1 ==> from ecliptic to reference frame +! v_in : input vector v_in(3) in the initial frame +! eps : Epsilon, inclination of ref frame [rad] +! om : Omega, longitude of node [rad] +! +! OUTPUT +! v_out : output vector in the final frame v_out(3) +! ierr : ierr = 0 :: normal exit +! ierr = 100 :: error ieqec neither 1 or -1 . +! +!*************************************************************** +!f2py intent(in) ieqec +!f2py intent(in) v_in +!f2py intent(in) eps +!f2py intent(in) om +!f2py intent(out) v_out +!f2py intent(out) ierr + implicit none + + type(t_v3d), intent(in) :: v_in + type(t_v3d), intent(out) :: v_out + integer, intent(in) :: ieqec + integer, intent(out) :: ierr + real (kind=8), intent(in) :: eps, om +! obliquity at J20000 in arcsec + real (kind=8), parameter :: Pi = 3.141592653589793238d0, & + secrad = Pi/180.d0/3600.d0 + real (kind=8) :: ww(3), coseps, sineps, cosom, sinom + + coseps = dcos(eps) + sineps = dsin(eps) + cosom = dcos(om) + sinom = dsin(om) + +! regular exit + ierr = 0 + +! to allow a call like :: call ref_ecl(ieqec,vv,vv,ierr) + ww(1) = v_in%x + ww(2) = v_in%y + ww(3) = v_in%z + + if (ieqec .eq. 1) then + v_out%x = cosom*ww(1) - sinom*(coseps*ww(2) - sineps*ww(3)) + v_out%y = sinom*ww(1) + cosom*(coseps*ww(2) - sineps*ww(3)) + v_out%z = sineps*ww(2) + coseps*ww(3) + else if (ieqec .eq. -1) then + v_out%x = cosom*ww(1) + sinom*ww(2) + v_out%y = coseps*(-sinom*ww(1) + cosom*ww(2)) + sineps*ww(3) + v_out%z = - sineps*(-sinom*ww(1) + cosom*ww(2)) + coseps*ww(3) + else +! anomalous exit ieqec not allowed + ierr = 100 + end if + + return + end subroutine ref_ecl + + subroutine ref_ecl_osc(ieqec, o_mi, o_mo, eps, om, ierr) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine converts osculating elements back and forth between +! given reference frame and ecliptic plane. +! Uses ref_ecl to do the work. +! Reference frame is given by Omega (om), the node longitude (rotation +! around Z axis of ecliptic) and Epsilon (eps), the inclination of the +! reference frame (rotation around node axis of ecliptic). +! +! ANGLES ARE GIVEN IN RADIAN !!!! +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) ieqec +!f2py intent(in) o_mi +!f2py intent(in) eps +!f2py intent(in) om +!f2py intent(out) o_mo +!f2py intent(out) ierr + implicit none + + type(t_orb_m), intent(in) :: o_mi + type(t_orb_m), intent(out) :: o_mo + integer, intent(in) :: ieqec + integer, intent(out) :: ierr + real (kind=8), intent(in) :: eps, om + real (kind=8), parameter :: Pi = 3.141592653589793238d0, TwoPi = 2.d0*Pi, & + drad = Pi/180.d0, mu = TwoPi**2 + type(t_orb_m) :: o_md + type(t_v3d) :: posi, poso, veli, velo + + o_md%a = o_mi%a + o_md%e = o_mi%e + o_md%inc = o_mi%inc + o_md%node = o_mi%node + o_md%peri = o_mi%peri + o_md%m = o_mi%m + call coord_cart (mu, o_md, posi, veli) + call ref_ecl (ieqec, posi, poso, eps, om, ierr) + call ref_ecl (ieqec, veli, velo, eps, om, ierr) + call osc_el (mu, poso, velo, o_mo) + call ztopi (o_mo%inc) + if (o_mo%inc .gt. Pi) o_mo%inc = o_mo%inc - Pi + call ztopi (o_mo%node) + call ztopi (o_mo%peri) + call ztopi (o_mo%m) + + return + end subroutine ref_ecl_osc + + subroutine forced_plane_damp(a, inc, ifd, Omfd) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* +! This routine returns the inclination iforced and node Omforced of the +! forced plane as a function of semimajor-axis in the region of the +! main classical Kuiper belt, according to second order theory with all +! 8 planets. +! +! This version damps the forced plane to the invariable plane according +! to inclination. +! +! The fit is valid only between 35 au and 47.74 au (2:1 MMR). +! Returns 0 for a < 35 au, and the invariable plane for a > 47.74 au +! +! Returns degrees. +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* +! +! J-M. Petit Observatoire de Besancon +! Version 1 : November 2021 +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* +! INPUT +! a : semimajor-axis [au] (R8) +! inc : Inclination [deg] (R8) +! +! OUTPUT +! ifd : Inclination of forced plane [deg] (R8) +! Omfd : Node of forced plane [deg] (R8) +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* +! +! Set of F2PY directives to create a Python module +! +!f2py intent(in) a +!f2py intent(in) inc +!f2py intent(out) ifd +!f2py intent(out) Omfd + implicit none + +! Calling arguments + real (kind=8), intent(in) :: a, inc + real (kind=8), intent(out) :: ifd, Omfd + +! Internal variables + real (kind=8) :: ci0(4), ci1(3), ci2(3), ci3(3), co0(3), co1(3), co2(5), & + co3(5), alpha, beta, gamma, delta, om_lim_low, om_lim_high, damp, y + real (kind=8), parameter :: epsilon = 5713.86d0/3600.d0, & + omega = 387390.8d0/3600.d0, ce_damp = 16.d0, wi_damp = 6.d0 + + data & + ci0 /-0.115657583d0,34.8097343d0,7.79198557d-02,-1.06252408d0/, & + ci1 /-0.392018467d0, 40.5282974d0, 0.137285933/, & + ci2 /-0.391531527d0, 40.6837921d0, 0.190758109/, & + ci3 /0.391087621d0, 40.2941360d0, 0.146276891/, & + co0 /971.35006612732775d0, -50.061061665365997d0, & + 0.74715152013989472d0/, & + co1 /4694.6567730091410d0, -249.62252879980477d0, & + 3.4214148236506192d0/, & + co2 /11.7058372d0, -294.139587d0, 30.8862801d0, & + -1049.41772d0, 15.2022839d0/, & + co3 /30.4611244d0, -1203.70593d0, 2.47961617d0, & + -16.6130295d0, 302.805359d0/, & + alpha /1.d0/, & + om_lim_low /200.d0/, & + om_lim_high /20.d0/ + + common /om_lim_com/ om_lim_low, om_lim_high + + damp(y) = 1.d0*(1.d0-tanh((y-ce_damp)/wi_damp))/2.d0 + + if (a .lt. 35.d0) then + ifd = 0.d0 + Omfd = 0.d0 + else if (a .lt. 37.d0) then + ifd = (ci0(1)/(a - ci0(2)) + ci0(3)*a + ci0(4)) + Omfd = co0(1) + co0(2)*a + co0(3)*a*a + else if (a .lt. 39.0d0) then + ifd = 10.d0**(ci1(1)/(a - ci1(2)) + ci1(3)) + Omfd = co1(1) + co1(2)*a + co1(3)*a*a + else if (a .lt. 40.0d0) then + ifd = 2.477d0 + Omfd = (om_lim_low-163.187d0)/(40.0d0-39.d0)*(a-40.0d0) & + + om_lim_low + else if (a .lt. 41.9d0) then + ifd = 2.477d0 + Omfd = (61.791d0-om_lim_high)/(41.9d0-40.0d0)*(a-40.0d0) & + + om_lim_high + else if (a .lt. 47.74d0) then + ifd = 10.d0**(ci3(1)/(a - ci3(2)) + ci3(3)) + beta = -((co3(1) + co3(3))*a + co3(2) + co3(4)) + gamma = (co3(1)*a + co3(2))*(co3(3)*a + co3(4)) - co3(5) + delta = max(0., beta**2 - 4.d0*alpha*gamma) + Omfd = (-beta - sqrt(delta))/(2.d0*alpha) + else + ifd = epsilon + Omfd = omega + end if + ifd = min(ifd, 40.d0) + ifd = epsilon + (ifd-epsilon)*damp(inc) + Omfd = omega + (Omfd-omega)*damp(inc) + + return + + end subroutine forced_plane_damp + + subroutine ztopi (var) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This subroutine resets var to be between 0 and 2Pi. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in,out) var + implicit none + + real (kind=8), intent(inout) :: var + real (kind=8), parameter :: Pi = 3.141592653589793238d0, TwoPi = 2.d0*Pi + +1000 continue + if (var .gt. TwoPi) then + var = var - TwoPi + goto 1000 + end if +1100 continue + if (var .lt. 0.d0) then + var = var + TwoPi + goto 1100 + end if + + return + end subroutine ztopi + +end module rot