From 5b6f642c5693b64cb5beda0d3dc177452cec21ee Mon Sep 17 00:00:00 2001 From: darksnake Date: Tue, 14 Jun 2016 16:48:47 +0300 Subject: [PATCH] Fixed propagation procedure --- build.gradle | 4 + gradle/wrapper/gradle-wrapper.jar | Bin 53556 -> 53319 bytes gradle/wrapper/gradle-wrapper.properties | 4 +- .../numass/trapping/SimulationManager.java | 3 +- .../java/inr/numass/trapping/Simulator.java | 211 +++++++++++------- .../java/inr/numass/trapping/Trapping.java | 4 +- 6 files changed, 136 insertions(+), 90 deletions(-) diff --git a/build.gradle b/build.gradle index db849df..1b609c6 100644 --- a/build.gradle +++ b/build.gradle @@ -17,6 +17,10 @@ tasks.withType(JavaCompile) { options.encoding = 'UTF-8' } +tasks.withType(JavaExec){ + enableAssertions = true +} + repositories { mavenCentral() } diff --git a/gradle/wrapper/gradle-wrapper.jar b/gradle/wrapper/gradle-wrapper.jar index ca78035ef0501d802d4fc55381ef2d5c3ce0ec6e..d3b83982b9b1bccad955349d702be9b884c6e049 100644 GIT binary patch delta 3698 zcmZWrc|4Tc8-ESc=rZ<^waLD3kz7lQEfGo9$dZw5$QCYXF1pGR@|Gl&B5SspC=nWC z%w%h$1sRmBvXpGa@11wT#1Jm-1NnR7^Y*r~VK(Y99X9KQhoCnsQN z*>eRg&j}fe9oix#R}5T)a_7S4!`Og-Ia&~>REnKxx)fovq{%d(G8=M4{^ZA1SS)~52-CqE&HoXZUX=@9srOB8CYSkKt!0} zAwIyrp>A=)r7}p>+rZPp=ZK&bTlr8&ko9P`O1%KbRTYYju`o$^=F+vg$908HF^{Zg zRg7LVF2#ZUR|flzev92Gt^e@$q}A<9@+${tyJK_re#@Wjy4kR?-aE5_2;tdCV;H8H zT&<$owwkA9DUs2l1y0AVezA6OR<_h?CwO++xD@TUQ{xkP%9cToG)uav983w3aMi*; zBTUTID(>4!R+yuS?E7nn*(xcexAW&)5u35_-RImXK^Tuqqp4c5b9O==0v|P=&0(~@ zJ6w4pOd0T5Hy8SlHtMQ%PW%J^gvO!yA;l*O3LcL0ahUzKA^1j8)eZ$~>VxhFDS?D+ zRpg2OI_}szr3UE?MZL3soX3a71f2Ev3OK(=#|HD>7Yf zvsGx8kIb3U!rX?!t}mXV*oNSO&;Lj+OQx)r$sk>OcU+dE)>L2BH6!o!&fcAI(=1V5 zIwLoG%H6K7K*G(DFIVwnS!BdJQrZnV;Zv@c#&n>W$xp6jpCMfzNqo{E8Z2Eav5J%$ zdrTj!TMTkx+|XupPLi) zTygIA-%wNE|YF`sf*H2s;$R-2=`7NOf#q6kjnVAp7_?!o?szBPvk02X%6m`SBdjL zU+&zFz8h^%Mui@usEEnWnV<*eNJ%A{6;i4G?&Hnbdwz08=@_E&of<_c(G=BBXjgKE z*@VKfk|@n@`9`Jrd^oy=Z|eHPyCXjhszYk$RWL36nu#vM(}UwO%IEJg2o1~mDIWrt z+1E=MN0*5S>b!5{yL~#kjr+qcC){35RZJVLPVhU_HZOd`vG3)Z{A%Z_NB*Bzh*5X! zl|BhN(yZ|{nu7Nnsvd4{(J^oTQIVGO^ut^fI%+;QsoeCO)8JX%y$8^H2E7*&dpc5` zytKt>H|tXfe$gidG*eFY7vKo%rxJ>UC)fV(FJgz1BsKx0+(co>pjW1dY?MG`>WC3zM*9Gp{(oAQ=1br^v7u8 zuEYTfE;Yo0CRI`w*%DZA?79dUp zV*#9~@3KQgfrsap$qpI9LjXBDQI0-bwZFzUe>5X5Gwf0#KT=@kXj za_L$<2e8J`eAqihE0)We9eSzp^R2C8uY`tjL@bIuOkBxy{-=oj#l{!@;QUok3_gkU zJS-b8w@Y4A^vzEVqq4j+B`+LE1AO%2qyWy|qNJC5Ox2Sj3#&Q_9?!*izDU>GnV_6b zU@ndw(u=sLXCF#R>Pcm+HBx^}bu3NAys#a!#qFRv*$pwiTA)uFEQZAJF>-Kxq^tAi~^~~(!HOs)sU3&zJ7Pi~*-ODFO zXW3oIJ(q{rO?UfCkHA}Wsirwaa8F*Izv*Zq_D9>Sj@+DtVAL=vSQhQ(I(QRUGg3Jt04mOr_* zUvDJ&u5ce;uZkPI94c`L<^2^`w>Gt+8priyNo$=mM(bk6oGAgda)XMf%`0xYB&U&k zMSqoFHFom-xEychKC2+5hN_64cno{S;_yMeW-xg6{^$_eebim&= za8^_kI)Xkw(g;gD=T-a z(7LQ2Apa9s2~r+WnB$0ml*s_7=7}YA*77{v#Rv2?x12ilRnb_~dC<-6J1TvK|FETh z%&XafoaK(|>2+Pbc4h{z)x-{2=i$Qg;CZgeQ-ug7y`` z`h-d9f#i#KGOp>yJ=m$Orc8E?R10d(vefYM!vSjN4<-Jtdp~Q#| zUY8Ke0XX741p{5U4HzKy_)f6~jR7DVJY|c=a`H1*BM~TE+l>H#_s|1lFIY&F&GOqu zERUFA@vlZ1u=6l!(lN9XSx5l5fd*?5Fj+V*;;)!eBP@H&O(l* zRj~YQq?Zc7Q{{#8ETo`lgf)mSspo*F^efrf;lY9GIo6bR-AR@oka&>A8EEKZA?1(k zSsb<|kR8p;kavUZBm}iZ04!{gm)UIS$N#pG4I#w7`K$vaTXi7Hl?1&TrT@`DmJ}~!^LJPjxJD6z$pO0pGD&Oz z;K&XDdRxff?O}3(2paTo5CQ!=xM4w-x23o!CjeOS0)W&O!MYpBEr^5$uS zh%#HqjvnAG4>@o9 zO$-8a0pQ4={NMHgK+s6=4H{$8N*S{oOlZrf0V}V&E<_mn_;g zu(F#Kw+%uY??L@BV(w~~sN5rt2(N&$Q{P-&eQQm44wOzgv}?Ck^N%X9u;=e;LOOgc z6aeZ%9WT=zycCn?!k?z0s+{8j0K+Xy>8gVhOzLHIyQWv1jSuYVWmP~8*e&F~C*ZF> Kt!+{*O!R*iz94e| delta 4032 zcma)9dpwle8vYD}ouOO?nPG?+cQaHX+EMOuCrN}5!nmd|8^w+yeQdRdNGK7xq#882 zjLAJIaw#E~L`HUW<$Uwa-a2#6`Qyy*H@|m%&-Xshds%DETEi&dt}Wyiw6)~s*#rRm z{6I9}YNntfKV;1K7Io2?`mSP9%VA4lT)5!;bF6?#WJ7@*Eq}82qD4| zFmRS<1_pwWMKDl<<`?CK5mlifHo_`DQCXzcm$fE` ziQ-q|co13NcqI-ix7KpSzs*2A8gh55Q|VvfYuvSX;&s<9?J>&z^|js~ySQ@OXcM}r z%`&-kdn|8ot8$~Pv0d^M-fF1Mr%HsLMrtM7U6OWFP~vr&w#V9yW-Yi8%T21p8T(=t z(la{9eoOW!T-wAWDHf4`VcdZfE0v!~_d7t9rqB4S#bO!EYLv(C6KQs_%E_9HzH^nW zOeedwKT#eMl(8-ED1=zc=@14pPx29lvC52VA}J-~TA9YXW0X?sHwj%hly;c(gdBQg zt2wZ zw!o0Gl$O8@dS1HLit)lK+FA9F2;?_KMMgeWCik`8y-Lr(FNKDavn$5~U3Fr$le6QR z@pC%DF{4_7tA&Yej;rN~PYbed4wKh}Tb?s4It8c(chtoqmo(q_?HIWfd3S4A4yy12 zwrbvJAUeD8c%Z5(H81b!XmPyBJu7QUx=>LRRXpUI0JUBlw_JLd_w&eGvULT!zg|&1e_`sxfCMF|(@V0BmfC?kf_z3sdmm^z zHFxL5841!rX>6^Jgl1rm$oA={IWiTK#UI8SJE@cY$har?+4Qo47MHNU=$B7E2psBL zAUo`BxI4W`A5$TOx0Rf{e$qL`MdDMo>|S9;cvJK6*Xl?I+|kI${-tMFV&go%EBUf_ zPh!Wn-S1JCC{#+`)vg7R>Z0&AvPN*`DF;fM%b_SFNlM!DT*G|yhieyAz8Txw%H&pC ztQ2Zce7*_^eQN&&EfH!*kJ%$b&u~OdGE}a$>_CT~I!t$voILP4*p+U>Oie|qeY6-X zl*DhPHw2e`E}i|3+UsXhtr9O+Ch)m9)#LnMfenRky6XKNy!>n}Msp%0+$6hF!GkLb zmV-wQ5FUjZdqmv$oR>N_xw+x<4dV=xB~vuHIp3*#(g{s&@<7*~GX3M7hQ6m*SV@@# z&!@n`P_9aWiR#3p$SVZ}cc}WVM0ao3u@%2>1x_l{V9AtXj@f9@v5KtGBAYzJD)(b1 z_|cXz$5BRdD$Ua5Lt6dxFfX}uLa@O+j3+NcqNde%FlX01Ugp;46tP*vntAXh?Vj!C zjYZ=MmeT?yk!hZnLvp;3H(!riZ>Ks&tMAx6Q2T%@BUHVVDpu)qjp;8f?BaWL-=aSC z4HvD*{o#OPs9b#jHALFDW)11p@{zzRfHNk)del(E7+k*FI8rv^P?JQyOu@?!s#DH7 zFcb7Lk4l$Lh5j&Hdb)nXH-54`CS!V=_nyWxCeaJYgm=VvV^lGw}f%(13bOU?b1YCGCVa%U|a{M-p85^p~j&W%Zzy z@s`pfakP&CJoC>y)J`0@Np5@Xn^-wt-v8(1TxYz9Qc#*$=wK8VF>UhuYrT$;dXw6T zAE!{|DlURgLgRu*s5yJ1oI3~~s`|pd%erKQI_swNqw$reZ{sE;XHXHct&27B7Pl+U za-qV+xs4MNlx`_gNrX#14e80Zv*XSp6;4d)t`zjwM7$gnQai73mxtV{b@I#3bxRw& zp8Ne*iBhWk7fa5l@MjdC(MIp@Ik$659!T>-&ky|?Tal<(@tnA#xI*$5B;mcu$nKeD z3Y9NvRxxV!w7tzITOV$|nP+##2h5P_@7>00<(dg}TWS^V1oqe;e0@nyzAi@f^GUBa zE$q}FWgKa;d{QMYxF?bDLA68j4+;YRqytECTZ7oQu{Rt#o+(6NkaGtOSI3AwH}EoZyOoVmkRJAt&> zDXGCPVXS_SlF~PkldH_9cyPdQpx@1e!JFQF?)NDJv;xH@ixeA2j^K6=KQCI|*Uoos zYEms~CcXSqIQCajf)1#crxqpAss6c#_=ta)a9;VEzn5S2lg_ob<4f z-q~1Nu2=eG>BZ$kjFjwqv0@$ToX2xb3C>v*K2_I=C?q@((iIng2TI1Po8W=m!B_0F zn9GTeeHIJ2IPnXzj^M$(R~~i%0N8*hJycFO^mN-OcbhwznOk>ub*i;@^-kKaboV-h zs^RUmh)!BsHAJ0SBGEZCN>Ip7i>ReDseP{ELx$w=h<}lT zk-lT4S-=*1S^r<24D(C$Ftt2hm|FM#$IMxpr zLX3Ok!U}59bT)3RU=U4@wY}yTQ+|nl0eHv2i^Gby`f}H}NePL>mqyJVt|eTeiYb7X;z7VQ8&IQI@r* z%m(f&^g}tzO9B8MYRNyJ+5`<AK}ULTLba9r*DO(mLS~U^bQQ-N+Ot1whip=5P|t=`*sAyLVN*;ug4K1G!j zCjo{YLD0u^2-*hz=wwTIY&Lv$1FB6Ks!f|iD(^bXweOOFW19B)&R>QcbD`Qe9W#sm z&Y2g)7NOK^Er0H7RzprHPB?n<0DzPL0O)f>h`0$YBy%`3gKi8$y&Ni1V`EP19^VbW z0tLJcO%vywYFBOqd{iO6F|SHIW~i-DSNvGza0F7Uh6q?ttrG(WbM3((l51cMff|hJ ze reportIf = (res) -> res.state == Simulator.EndState.ACCEPTED; System.out.printf("%nStarting sumulation with initial energy %g and %d electrons.%n%n", initialE, num); + output.printf("%s\t%s\t%s\t%s\t%s\t%s%n", "E", "theta", "theta_start", "colNum", "L", "state"); Stream.generate(() -> getRandomTheta()).limit(num).parallel() .forEach((theta) -> { double initZ = (generator.nextDouble() - 0.5) * Simulator.SOURCE_LENGTH; @@ -144,7 +145,7 @@ public class SimulationManager { } private void printOne(PrintStream out, Simulator.SimulationResult res) { - out.printf("%g\t%g\t%g\t%d\t%s%n", res.E, res.theta * 180 / Math.PI, res.initTheta * 180 / Math.PI, res.collisionNumber, res.state.toString()); + out.printf("%g\t%g\t%g\t%d\t%g\t%s%n", res.E, res.theta * 180 / Math.PI, res.initTheta * 180 / Math.PI, res.collisionNumber, res.l, res.state.toString()); } diff --git a/src/main/java/inr/numass/trapping/Simulator.java b/src/main/java/inr/numass/trapping/Simulator.java index 70a462b..0ad8c25 100644 --- a/src/main/java/inr/numass/trapping/Simulator.java +++ b/src/main/java/inr/numass/trapping/Simulator.java @@ -45,7 +45,7 @@ public class Simulator { setELow(Elow); } - public static enum EndState { + public enum EndState { ACCEPTED,//трэппинговый электрон попал в аксептанс REJECTED,//трэппинговый электрон вылетел через заднюю пробку @@ -64,14 +64,30 @@ public class Simulator { this.thetaPinch = Math.asin(Math.sqrt(Bsource / Bpinch)); } + /** + * Set gas density in 1/m^3 + * + * @param gasDensity + */ public void setGasDensity(double gasDensity) { this.gasDensity = gasDensity; } + /** + * Longitudal magnetic field distribution + * + * @param magneticField + */ public void setFieldFunc(UnivariateFunction magneticField) { this.magneticField = magneticField; } + /** + * Perform scattering in the given position + * + * @param pos + * @return + */ private State scatter(State pos) { //Вычисляем сечения и нормируем их на полное сечение double sigmaIon = scatter.sigmaion(pos.e); @@ -116,13 +132,13 @@ public class Simulator { * @return */ private double freePath(double e) { - //FIXME double cross-section calculation + //FIXME redundant cross-section calculation //All cross-sections are in m^2 - return new ExponentialDistribution(generator, 1 / scatter.sigmaTotal(e) / gasDensity).sample(); + return new ExponentialDistribution(generator, 1d / scatter.sigmaTotal(e) / gasDensity).sample(); } /** - * Calculate propagated position + * Calculate propagated position before scattering * * @param deltaL * @return z shift and reflection counter @@ -131,23 +147,27 @@ public class Simulator { // if magnetic field not defined, consider it to be uniform and equal bSource if (magneticField == null) { double deltaZ = deltaL * cos(pos.theta); // direction already included in cos(theta) + double z0 = pos.z; + pos.addZ(deltaZ); + pos.l += abs(pos.z - z0) / cos(pos.theta); - //if we are crossing source boundary, check for end condition - if(abs(deltaZ + pos.z)>SOURCE_LENGTH/2d) { - checkEndState(pos); - } - - // if track is finished apply boundary position - if (pos.isFinished()) { - // remembering old z to correctly calculate total l - double oldz = pos.z; - pos.z = pos.direction() * SOURCE_LENGTH / 2d; - pos.l += (pos.z - oldz)/cos(pos.theta); - } else { - //else just add z - pos.l += deltaL; - pos.addZ(deltaZ); - } +// //if we are crossing source boundary, check for end condition +// while (abs(deltaZ + pos.z) > SOURCE_LENGTH / 2d && !pos.isFinished()) { +// +// pos.checkEndState(); +// } +// +// // if track is finished apply boundary position +// if (pos.isFinished()) { +// // remembering old z to correctly calculate total l +// double oldz = pos.z; +// pos.z = pos.direction() * SOURCE_LENGTH / 2d; +// pos.l += (pos.z - oldz) / cos(pos.theta); +// } else { +// //else just add z +// pos.l += deltaL; +// pos.addZ(deltaZ); +// } return pos; } else { @@ -161,31 +181,45 @@ public class Simulator { double b = field(pos.z); double root = 1 - sin2 * b / bSource; - // change direction in case of reflection. Loss of precision here? + //preliminary reflection if (root < 0) { - //flip direction pos.flip(); - // check if end state occured - checkEndState(pos); - // finish if it does - if(pos.isFinished()){ - return pos; - } else { - // move in reversed direction - pos.z += pos.direction() * delta * sqrt(-root); - } - } else { - // move - pos.z += pos.direction() * delta * sqrt(root); } + pos.addZ(pos.direction() * delta * sqrt(abs(root))); +// // change direction in case of reflection. Loss of precision here? +// if (root < 0) { +// // check if end state occurred. seem to never happen since it is reflection case +// pos.checkEndState(); +// // finish if it does +// if (pos.isFinished()) { +// //TODO check if it happens +// return pos; +// } else { +// //flip direction +// pos.flip(); +// // move in reversed direction +// pos.z += pos.direction() * delta * sqrt(-root); +// } +// +// } else { +// // move forward +// pos.z += pos.direction() * delta * sqrt(root); +// //check if it is exit case +// if (abs(pos.z) > SOURCE_LENGTH / 2d) { +// // check if electron is out +// pos.checkEndState(); +// // finish if it is +// if (pos.isFinished()) { +// return pos; +// } +// // PENDING no need to apply reflection, it is applied automatically when root < 0 +// pos.z = signum(pos.z) * SOURCE_LENGTH / 2d; +// if (signum(pos.z) == pos.direction()) { +// pos.flip(); +// } +// } +// } - //normalize in case reflection is beyound the source - if (abs(pos.z) > SOURCE_LENGTH / 2d) { - pos.z = signum(pos.z) * SOURCE_LENGTH / 2d; - if(signum(pos.z)== pos.direction()){ - pos.flip(); - } - } curL += delta; pos.l += delta; @@ -194,30 +228,6 @@ public class Simulator { } } - /** - * Check if this position is an end state and apply it if necessary - * - * @param pos - * @return - */ - private State checkEndState(State pos) { - //accepted by spectrometer - if (pos.theta < thetaPinch) { - if (pos.colNum == 0) { - //counting pass electrons - pos.setEndState(EndState.PASS); - } else { - pos.setEndState(EndState.ACCEPTED); - } - } - - //through the rear magnetic pinch - if (pos.theta >= PI - thetaTransport) { - pos.setEndState(EndState.REJECTED); - } - return pos; - } - /** * Magnetic field in the z point * @@ -246,11 +256,12 @@ public class Simulator { while (!pos.isFinished()) { double dl = freePath(pos.e); // path to next scattering // propagate to next scattering position - propagate(pos,dl); + propagate(pos, dl); if (!pos.isFinished()) { // perform scatter scatter(pos); + // increase collision number pos.colNum++; if (pos.e < eLow) { //Если энергия стала слишком маленькой @@ -363,20 +374,53 @@ public class Simulator { */ double addZ(double dZ) { this.z += dZ; - while (abs(this.z) > SOURCE_LENGTH / 2d) { + while (abs(this.z) > SOURCE_LENGTH / 2d && !isFinished()) { + checkEndState(); + if (!isFinished()) { + flip(); + } // reflecting from back wall if (z < 0) { - // reflecting from rear pinch - z = -SOURCE_LENGTH - z; + if (isFinished()) { + z = -SOURCE_LENGTH / 2d; + } else { + // reflecting from rear pinch + z = -SOURCE_LENGTH - z; + } } else { - // reflecting from forward transport magnet - z = SOURCE_LENGTH - z; + if (isFinished()) { + z = SOURCE_LENGTH / 2d; + } else { + // reflecting from forward transport magnet + z = SOURCE_LENGTH - z; + } } - flip(); } return z; } + /** + * Check if this position is an end state and apply it if necessary. Does not check z position. + * + * @return + */ + private void checkEndState() { + //accepted by spectrometer + if (theta < thetaPinch) { + if (colNum == 0) { + //counting pass electrons + setEndState(EndState.PASS); + } else { + setEndState(EndState.ACCEPTED); + } + } + + //through the rear magnetic pinch + if (theta >= PI - thetaTransport) { + setEndState(EndState.REJECTED); + } + } + /** * Reverse electron direction */ @@ -406,9 +450,8 @@ public class Simulator { if (theta > PI / 2) { newTheta = PI - newTheta; } - if (Double.isNaN(newTheta)) { - throw new Error(); - } + + assert !Double.isNaN(newTheta); return newTheta; } } @@ -437,13 +480,13 @@ public class Simulator { double newtheta = acos(result.getZ()); - //следим чтобы угол был от 0 до Pi - if (newtheta < 0) { - newtheta = -newtheta; - } - if (newtheta > Math.PI) { - newtheta = 2 * Math.PI - newtheta; - } +// //следим чтобы угол был от 0 до Pi +// if (newtheta < 0) { +// newtheta = -newtheta; +// } +// if (newtheta > Math.PI) { +// newtheta = 2 * Math.PI - newtheta; +// } //change back to virtual angles if (magneticField == null) { @@ -455,9 +498,7 @@ public class Simulator { } } - if (Double.isNaN(theta)) { - throw new Error(); - } + assert !Double.isNaN(theta); return theta; } diff --git a/src/main/java/inr/numass/trapping/Trapping.java b/src/main/java/inr/numass/trapping/Trapping.java index d236c24..139ef92 100644 --- a/src/main/java/inr/numass/trapping/Trapping.java +++ b/src/main/java/inr/numass/trapping/Trapping.java @@ -16,10 +16,10 @@ public class Trapping { System.out.printf("Starting at %s%n%n", startTime.toString()); new SimulationManager() .withParameters(0.6, 3.7, 4.84, 18000d, 4000) -// .withOutputFile("D:\\Work\\Numass\\trapping\\trap 18, pinch 100A, fields.out") + .withOutputFile("D:\\Work\\Numass\\trapping\\trap 18, pinch 100A.out") // .withFieldMap(z, b) // .withDensity(5e20) - .simulateAll((int) 1e4); + .simulateAll((int) 1e6); Instant finishTime = Instant.now(); System.out.printf("%nFinished at %s%n", finishTime.toString()); System.out.printf("Calculation took %s%n", Duration.between(startTime,finishTime).toString());