-- $Id: ion.occ 1795 2004-06-28 11:19:13Z Martin vonLoewis $ #INCLUDE "mathvals.inc" #INCLUDE "filelib.inc" #USE "dblmath.lib" #USE "course.lib" #USE "file.lib" VAL NUM.WORKERS IS 2: VAL NUM.IONS IS 1000000: VAL IONS.PER.WORKER IS NUM.IONS/NUM.WORKERS: DATA TYPE MEAS RECORD REAL64 abstand: REAL64 bfeld: : VAL []MEAS m IS [[ 490.0, 2.765 ], [ 488.0, 2.935 ], [ 486.0, 3.030 ], [ 484.0, 3.165 ], [ 482.0, 3.275 ], [ 480.0, 3.331 ], [ 478.0, 3.364 ], [ 476.0, 3.313 ], [ 474.0, 3.298 ], [ 472.0, 3.276 ], [ 470.0, 3.121 ], [ 468.0, 3.047 ], [ 466.0, 2.924 ], [ 464.0, 2.770 ], [ 462.0, 2.587 ], [ 460.0, 2.400 ], [ 458.0, 2.198 ], [ 456.0, 1.980 ], [ 454.0, 1.783 ], [ 452.0, 1.582 ], [ 450.0, 1.375 ], [ 448.0, 1.182 ], [ 446.0, 1.011 ], [ 444.0, 0.848 ], [ 442.0, 0.703 ], [ 440.0, 0.583 ], [ 438.0, 0.481 ], [ 436.0, 0.390 ], [ 434.0, 0.315 ], [ 432.0, 0.249 ], [ 430.0, 0.187 ], [ 428.0, 0.143 ], [ 426.0, 0.104 ], [ 424.0, 0.072 ], [ 422.0, 0.045 ], [ 420.0, 0.027 ], [ 418.0, 0.009 ]]: VAL MAX.VALUE IS 10000: REAL64 FUNCTION bfeldprofil(VAL REAL64 r) VAL INT n IS INT TRUNC (r*500.0): REAL64 res: VALOF IF n <= 30 VAL REAL64 p1 IS m[n+6][bfeld] / 3.364: VAL REAL64 p2 IS m[n+7][bfeld] / 3.364: VAL REAL64 anteil2 IS (r - ((REAL64 ROUND n)/500.0)) * 500.0: res := (anteil2 * p2) + ((1.0 - anteil2) * p1) TRUE res := 0.0 RESULT res : PROC magnetfeld([NUM.WORKERS]CHAN OF REAL64 in, out, CHAN OF BOOL eof) [MAX.VALUE]REAL64 v: BOOL stopped: SEQ SEQ i=0 FOR MAX.VALUE VAL REAL64 r IS (REAL64 ROUND i) * 1.0E-7: v[i] := bfeldprofil(DSQRT(r)) stopped := FALSE WHILE NOT stopped REAL64 r2: ALT eof ? stopped SKIP ALT i=0 FOR NUM.WORKERS in[i] ? r2 VAL INT x IS (INT TRUNC (r2 * 1.0E+7)): IF x >= MAX.VALUE out[i] ! 0.0 TRUE out[i] ! v[x] : PROC rand.r64([3]CHAN OF BOOL query, [3]CHAN OF REAL64 numbers) INT64 seed: REAL64 num: BOOL stopped: SEQ seed := 8 stopped := FALSE WHILE NOT stopped SEQ num, seed := DRAN(seed) ALT i=0 FOR 3 query[i] ? stopped numbers[i] ! num : PROC rand.gauss(CHAN OF BOOL q.rand, CHAN OF REAL64 to.coord, v.rand, CHAN OF BOOL from.coord) REAL64 r1, r2: REAL64 x2pi, g2rad: REAL64 v1, v2: BOOL stopped: SEQ stopped := FALSE WHILE NOT stopped SEQ q.rand ! FALSE v.rand ? r1 q.rand ! FALSE v.rand ? r2 x2pi := r1 * (2.0 * DPI) g2rad := DSQRT((-2.0) * DALOG(1.0-r2)) v1 := DCOS(x2pi) * g2rad v2 := DSIN(x2pi) * g2rad from.coord ? stopped to.coord ! v1 IF NOT stopped SEQ from.coord ? stopped to.coord ! v2 TRUE SKIP : PROTOCOL RNG CASE uniform gauss.energy gauss.initial.y : PROC rand.coord([NUM.WORKERS]CHAN OF RNG query, [NUM.WORKERS]CHAN OF REAL64 values, CHAN OF BOOL eof) [3]CHAN OF BOOL r.query: CHAN OF BOOL g.query.1, g.query.2: [3]CHAN OF REAL64 r.result: CHAN OF REAL64 g.result.1, g.result.2: REAL64 val: BOOL stopped: PAR rand.r64(r.query, r.result) rand.gauss(r.query[1], g.result.1, r.result[1], g.query.1) rand.gauss(r.query[2], g.result.2, r.result[2], g.query.2) SEQ stopped := FALSE WHILE NOT stopped ALT eof ? stopped SEQ g.query.1 ! TRUE g.result.1 ? val g.query.2 ! TRUE g.result.2 ? val r.query[0] ! TRUE r.result[0] ? val ALT i=0 FOR NUM.WORKERS query[i] ? CASE uniform SEQ r.query[0] ! FALSE r.result[0] ? val values[i] ! val gauss.energy SEQ g.query.1 ! FALSE g.result.1 ? val values[i] ! val gauss.initial.y SEQ g.query.2 ! FALSE g.result.2 ? val values[i] ! val : DATA TYPE measpoint RECORD INT num: REAL64 msum: REAL64 csum: : VAL REAL64 r2fenster IS 0.000049: -- Fenster hat 7mm Radius VAL INT WIDTH IS 1024: VAL INT HEIGHT IS 1024: VAL REAL64 SCALE IS 0.000014: VAL REDUCE IS 2: DATA TYPE measarray IS [WIDTH][HEIGHT]measpoint: DATA TYPE whitearray IS [WIDTH][HEIGHT]BOOL: PROTOCOL photon CASE add;REAL64; REAL64; REAL64; INT eof : PROC ppm(CHAN OF photon from.ion, CHAN OF measarray res) measarray points: BOOL stopped: SEQ SEQ x=0 FOR WIDTH SEQ y=0 FOR HEIGHT SEQ points[x][y] := [0, 0.0, 0.0] stopped := FALSE WHILE NOT stopped REAL64 x,y,winkel: INT detect, w, h: from.ion ? CASE add;x;y;winkel;detect SEQ w := (INT TRUNC (x / SCALE)) + (WIDTH/2) h := (INT TRUNC ((-y) / SCALE)) + (HEIGHT/2) IF ((w >= 0) AND (w < WIDTH)) AND ((h >= 0) AND (h < HEIGHT)) measpoint p IS points[w][h]: SEQ p[num] := p[num] + 1 p[msum] := p[msum] + DCOS(winkel * 2.0) p[csum] := p[csum] + DSIN(winkel * 2.0) TRUE SKIP eof stopped := TRUE res ! points : PROC ppm.write(CHAN OF BYTE data, CHAN OF BOOL eof) INT fd: BYTE d: BOOL e, closed: INT res: SEQ file.open3("output.ppm", O.RDWR BITOR O.CREAT, S.I644, fd) closed := FALSE WHILE NOT closed ALT data ? d file.write(fd, [d], res) eof ? e SEQ file.close(fd) closed := TRUE : INT FUNCTION MIN(VAL INT x,y) INT res: VALOF IF x maxval maxval := count TRUE SKIP maxval.r := REAL64 ROUND maxval mgesamt := 0.0 cgesamt := 0.0 igesamt := 0 SEQ j=0 FOR HEIGHT/REDUCE SEQ i=0 FOR WIDTH/REDUCE INT count, rval, gval, bval: BOOL white: REAL64 msum, csum, pp, coswinkel, sinwinkel, tmp: SEQ count := 0 white := FALSE msum := 0.0 csum := 0.0 SEQ kx=0 FOR REDUCE SEQ ky=0 FOR REDUCE measpoint p IS total[(i*REDUCE)+kx][(j*REDUCE)+ky]: SEQ IF weiss[(i*REDUCE)+kx][(j*REDUCE)+ky] white := TRUE TRUE SKIP count := count + p[num] msum := msum + p[msum] csum := csum + p[csum] igesamt := igesamt + p[num] mgesamt := mgesamt + p[msum] cgesamt := cgesamt + p[csum] IF white rval, gval, bval := 255, 255, 255 count = 0 rval, gval, bval := 0, 0, 0 TRUE REAL64 punkt.gesetzt: VAL REAL64 MAGNIFY IS 270.0: SEQ msum := msum / (REAL64 ROUND count) -- also m/i csum := csum / (REAL64 ROUND count) pp := DSQRT((msum*msum) + (csum*csum)) coswinkel := DACOS(msum / pp) sinwinkel := DASIN(csum / pp) IF sinwinkel < 0.0 coswinkel := -coswinkel -- jetzt haben wir den gesamten bereich -pi bis pi TRUE SKIP punkt.gesetzt := REAL64 ROUND count tmp := DCOS(coswinkel) rval := INT TRUNC (((punkt.gesetzt * MAGNIFY) * (tmp * tmp)) / maxval.r) tmp := DCOS((DPI / 3.0) + coswinkel) gval := INT TRUNC (((punkt.gesetzt * MAGNIFY) * (tmp * tmp)) / maxval.r) tmp := DCOS(((2.0 * DPI) / 3.0) + coswinkel) bval := INT TRUNC (((punkt.gesetzt * MAGNIFY) * (tmp * tmp)) / maxval.r) rval := MIN(rval, 255) gval := MIN(gval, 255) bval := MIN(bval, 255) data ! BYTE rval data ! BYTE gval data ! BYTE bval out.string("MGESAMT ", 0, out) out.real64(mgesamt, 0, 0, out) out.string("*nCGESAMT ", 0, out) out.real64(cgesamt, 0, 0, out) out.string("*nIGESAMT ", 0, out) out.int(igesamt, 0, out) out ! '*n' eof1 ! TRUE eof2 ! TRUE eof3 ! TRUE : PROC init.white(CHAN OF whitearray res) whitearray weiss: VAL INT WINDOW.RESOLUTION IS 1000: SEQ SEQ x=0 FOR WIDTH SEQ y=0 FOR HEIGHT weiss[x][y] := FALSE SEQ i = 0 FOR WINDOW.RESOLUTION VAL REAL64 pi.2.i IS (DPI * 2.0) * (REAL64 ROUND i): VAL REAL64 pi.2.i.wr IS pi.2.i / (REAL64 ROUND WINDOW.RESOLUTION): REAL64 x,y: INT w,h: SEQ x := DSQRT(r2fenster) * DSIN(pi.2.i.wr) y := DSQRT(r2fenster) * DCOS(pi.2.i.wr) w := (INT TRUNC (x / SCALE)) + (WIDTH / 2) h := (INT TRUNC ((-y) / SCALE)) + (HEIGHT / 2) IF (((w>=0) AND (w < WIDTH)) AND ((h>=0) AND (h < HEIGHT))) weiss[w][h] := TRUE TRUE SKIP weiss[WIDTH/2][HEIGHT/2] := TRUE res ! weiss : PROC ion(CHAN OF RNG r.query, CHAN OF REAL64 randoms, to.magnetfeld, from.magnetfeld, CHAN OF photon to.ppm) VAL REAL64 q IS 1.6021892E-19: -- e in C VAL REAL64 m IS 1.6605655E-27 * 4.0: -- 4u in kg VAL q.div.m IS q / m: VAL REAL64 dt IS 1.0E-10: -- 100 ps REAL64 v0: VAL REAL64 b IS -0.005: -- 5 mT REAL64 vx, vy: REAL64 rx, ry: REAL64 alpha.sin, alpha.cos: VAL REAL64 alpha IS 85.0: -- 5 Grad Schraeglage der Oberflaeche VAL REAL64 target.y.top IS 0.004: -- Target ist 8mm lang VAL REAL64 target.y.bot IS -0.004: VAL REAL64 target.x IS 0.0: -- 1mm hinter Mitte VAL REAL64 deltapol IS ((DPI / 180.0) * (18.0 / 0.001)) * (dt / 10.0E-9): -- etwa 18 Grad pro Millitesla pro 10ns BOOL reflected, stopped: REAL64 rnd: INT i: SEQ SEQ n=0 FOR IONS.PER.WORKER SEQ r.query ! gauss.energy randoms ? rnd v0 := DSQRT((2.0 * q.div.m) * (12000.0 + (100.0 * rnd))) -- 12keV mit 200eV FWHM vx := -v0 -- Bewegung in "Minus-x"-Richtung vy := 0.0 rx := 0.200 -- 20cm ausserhalb des Mittelpunktes r.query ! gauss.initial.y randoms ? rnd ry := 0.000 + (0.00025 * rnd) -- 0mm ueber dem Nullpunkt, 1/2mm FWHM Strahldurchmesser alpha.sin := DSIN((DPI / 180.0) * alpha) alpha.cos := DCOS((DPI / 180.0) * alpha) reflected := FALSE stopped := FALSE i := 0 WHILE NOT stopped REAL64 b.r, rot.x, rot.y, r2, mag: SEQ rx := rx + (vx * dt) ry := ry + (vy * dt) r2 := (rx*rx) + (ry*ry) to.magnetfeld ! r2 from.magnetfeld ? mag b.r := b * mag vx := vx + (((q.div.m * dt) * b.r) * vy) vy := vy - (((q.div.m * dt) * b.r) * vx) rot.x := (alpha.cos * rx) + (alpha.sin * ry) rot.y := (alpha.cos * ry) - (alpha.sin * rx) IF ((rot.x < target.x) AND (rot.y < target.y.top)) AND (rot.y > target.y.bot) SEQ IF (rot.y < (target.y.top - 0.0001)) AND (rot.y > (target.y.bot + 0.0001)) reflected := TRUE TRUE SKIP stopped := TRUE TRUE SKIP IF i < 4000 i := i + 1 TRUE stopped := TRUE IF reflected REAL64 tmp.vx, tmp.vy: REAL64 winkel: REAL64 polwinkel: SEQ winkel := (DPI / 180.0) * alpha polwinkel := 0.0 tmp.vx := (DCOS(winkel) * vx) + (DSIN(winkel) * vy) tmp.vy := (DCOS(winkel) * vy) - (DSIN(winkel) * vx) tmp.vx := tmp.vx * (-1.0) winkel := winkel * (-1.0) vx := (DCOS(winkel) * tmp.vx) + (DSIN(winkel) * tmp.vy) vy := (DCOS(winkel) * tmp.vy) - (DSIN(winkel) * tmp.vx) i:=0 stopped := FALSE WHILE NOT stopped REAL64 b.r, r2, mag: SEQ rx := rx + (vx * dt) ry := ry + (vy * dt) r2 := (rx*rx) + (ry*ry) to.magnetfeld ! r2 from.magnetfeld ? mag b.r := b * mag polwinkel := polwinkel + (b.r * deltapol) r.query ! uniform randoms ? rnd IF rnd < 0.0046 SEQ to.ppm ! add; rx; ry; polwinkel; 0 stopped := TRUE TRUE SKIP IF r2 > r2fenster stopped := TRUE TRUE SKIP i := i + 1 IF i < 1000 i := i+1 TRUE stopped := TRUE TRUE SKIP to.ppm ! eof : PROC main(CHAN OF BYTE keyboard, screen, err) [NUM.WORKERS]CHAN OF RNG q.rand: [NUM.WORKERS]CHAN OF REAL64 v.rand: [NUM.WORKERS]CHAN OF measarray parts: [NUM.WORKERS]CHAN OF REAL64 to.magnetfeld, from.magnetfeld: CHAN OF whitearray whitechan: CHAN OF BYTE data: CHAN OF BOOL eof1, eof2, eof3: PAR PAR i=0 FOR NUM.WORKERS CHAN OF photon to.ppm: PAR ion(q.rand[i], v.rand[i], to.magnetfeld[i], from.magnetfeld[i], to.ppm) ppm(to.ppm, parts[i]) rand.coord(q.rand, v.rand, eof3) magnetfeld(to.magnetfeld, from.magnetfeld, eof2) init.white(whitechan) ppm.write(data, eof1) ppm.exit(parts, screen, whitechan, data, eof1, eof2, eof3) :