1 /*
2 * $Id: rays.i,v 1.1 2005-09-18 22:06:05 dhmunro Exp $
3 * Yorick functions to manipulate 3-D rays on a cylindrical mesh.
4 */
5 /* Copyright (c) 2005, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
11 /*
12 form_rays -- convert from nrays lists to (5, 3, or 6)-by-nrays
13 best_rays -- convert to 5-by-nrays representation
14 dirt_rays -- convert to 3-by-nrays representation
15 internal_rays -- convert to 6-by-nrays representation
16 picture_rays -- return rays at centers of arbitrary "pixel mesh"
17
18 get_s0 -- compute origins for slimits (internally, slimits
19 is based on an s coordinate which is zero at the
20 point of closest approach to the origin, rather
21 than at the x,y,z point selected in the best_rays
22 or internal_rays coordinate systems)
23
24 plray -- plot a ray
25 plray_lims -- set limits for plray
26 */
27
28 /*= SECTION() create and plot ray sets for drat package ====================*/
29
form_rays(rays)30 func form_rays(rays)
31 /* DOCUMENT best= form_rays( [x, y, z, theta, phi] )
32 or dirt= form_rays( [x, y, theta] )
33 or internal= form_rays( [cos, sin, y, z, x, r] )
34 forms 5-by-nrays, 3-by-nrays, or 6-by-nrays ray representation
35 given individual lists of array coordinates. The [...]
36 operator builds an nrays-by-5, nrays-by-3, or nrays-by-6
37 array, which form_rays transposes. The "nrays" may represent
38 zero or more actual dimensions.
39 SEE ALSO: best_rays, dirt_rays, internal_rays, picture_rays
40 */
41 { return transpose(rays, 2); }
42
best_rays(rays)43 func best_rays(rays)
44 /* DOCUMENT best_rays(rays)
45 returns 5-element (x,y,z,theta,phi) representation of RAYS.
46 The first dimension of RAYS may be length 3, 5, or 6 to represent
47 the ray(s) in TDG/DIRT coordinates (x,y,theta), "best" coordinates
48 (x,y,z,theta,phi), or internal coordinates (cos,sin,y,z,x,r),
49 respectively. The first dimension of the result always has length 5.
50
51 The "best" coordinate system is the easiest to visualize:
52 (x,y,z) represents any point on the ray, while (theta,phi)
53 represents the ray direction in standard spherical coordinates
54 relative to the +z-axis. Namely, theta is the angle from the
55 +z-direction to the ray direction (between 0 and pi), and phi is
56 the counterclockwise angle from the +x-axis to the projection of
57 the ray direction into the xy-plane, assuming xyz is a right-handed
58 coordinate system.
59
60 As a specification of a ray, this system is doubly redundant because
61 the point (x,y,z) could be any point on the ray, and the underlying
62 mesh through which the ray propagates is cylindrically symmetric about
63 the z-axis.
64
65 However, the slimits parameter -- used to specify the points along
66 a ray where the transport integration starts and stops -- is
67 measured from the point (x,y,z) specified as a part of the
68 (x,y,z,theta,phi) ray coordinate. Thus, any change in the point
69 (x,y,z) on a ray must be accompanied by a corresponding change in
70 the slimits for that ray.
71
72 SEE ALSO: form_rays, dirt_rays, internal_rays, get_s0, picture_rays
73 */
74 {
75 dummy= use_origins(0);
76
77 dim= dimsof(rays)(2);
78 if (dim==5) return rays;
79
80 if (dim==3) {
81 /* convert from DIRT/TDG (x,y,theta) to (x,y,z,theta,phi) */
82 x= rays(1,..);
83 y= rays(2,..);
84 theta= rays(3,..);
85 return form_rays([x*cos(theta), y, x*sin(theta),
86 abs(theta), (theta>=0.0)*pi]);
87
88 } else if (dim==6) {
89 /* convert from internal (cos,sin,y,z,x,r) to (x,y,z,theta,phi) */
90 theta= atan(rays(2,..), rays(1,..));
91 y= rays(3,..);
92 z= rays(4,..);
93 x= rays(5,..);
94 return form_rays([x, y, z, abs(theta), (theta<0)*pi]);
95 }
96
97 return [];
98 }
99
dirt_rays(rays)100 func dirt_rays(rays)
101 /* DOCUMENT dirt_rays(rays)
102 returns 3-element (x,y,theta) representation of RAYS.
103 The first dimension of RAYS may be length 3, 5, or 6 to represent
104 the ray(s) in TDG/DIRT coordinates (x,y,theta), "best" coordinates
105 (x,y,z,theta,phi), or internal coordinates (cos,sin,y,z,x,r),
106 respectively. The first dimension of the result always has length 3.
107
108 The TDG/DIRT coordinate system is based on the coordinates (x,y)
109 in a plane normal to the ray. Unfortunately, the old TDG and DIRT
110 codes used an angle theta which has the opposite sense from the
111 "best" and internal coordinates. Therefore, conversion from
112 TDG/DIRT coordinates to internal coordinates will reverse the
113 sign of theta. Conversion from TDG/DIRT coordinates to "best"
114 coordinates always results in positive theta, but the angle phi
115 will be pi for positive input theta.
116
117 The slimits parameter -- used to specify the points along
118 a ray where the transport integration starts and stops -- is
119 measured from the point of closest approach of the ray described
120 by (x,y,theta) to the origin x=y=z=0. Therefore, slimits is
121 independent of the TDG/DIRT ray coordinate representation.
122
123 SEE ALSO: form_rays, best_rays, internal_rays, get_s0, picture_rays
124 */
125 {
126 dummy= use_origins(0);
127
128 dim= dimsof(rays)(2);
129 if (dim==3) return rays;
130
131 if (dim==5) {
132 /* convert from best (x,y,z,theta,phi) to DIRT/TDG (x,y,theta) */
133 x= rays(1,..);
134 y= rays(2,..);
135 z= rays(3,..);
136 theta= rays(4,..);
137 phi= rays(5,..);
138 cosp= cos(phi);
139 sinp= sin(phi);
140 sgn= double((cosp>=0.0)*2-1);
141 return form_rays(sgn*[(x*cosp+y*sinp)*cos(theta) - z*sin(theta),
142 (y*cosp-x*sinp), -theta]);
143
144 } else if (dim==6) {
145 /* convert from internal (cos,sin,y,z,x,r) to DIRT/TDG (x,y,theta) */
146 cost= rays(1,..);
147 sint= rays(2,..)
148 y= rays(3,..);
149 z= rays(4,..);
150 x= rays(5,..);
151 /* NOTE SIGN CHANGE IN THETA */
152 return form_rays([x*cost-z*sint, y, -atan(sint, cost)]);
153 }
154
155 return [];
156 }
157
internal_rays(rays)158 func internal_rays(rays)
159 /* DOCUMENT internal_rays(rays)
160 returns 6-element (cos,sin,y,z,x,r) representation of RAYS.
161 The first dimension of RAYS may be length 3, 5, or 6 to represent
162 the ray(s) in TDG/DIRT coordinates (x,y,theta), "best" coordinates
163 (x,y,z,theta,phi), or internal coordinates (cos,sin,y,z,x,r),
164 respectively. The first dimension of the result always has length 6.
165
166 The internal coordinates are what Drat uses internally to
167 describe the ray. The coordinate system is rotated about the
168 z-axis until the ray lies in a plane of constant y (there are at
169 least two ways to do this). The point (x,y,z) can be any point on
170 the ray, and r=sqrt(x^2+y^2) is the corresponding cylindrical radius.
171 The clockwise angle theta from the +z-axis to the ray direction
172 (which always lies in the zx-plane) determines cos=cos(theta) and
173 sin=sin(theta).
174
175 As a specification of a ray, this system is triply redundant because
176 the point (x,y,z) could be any point on the ray, both the sine and
177 cosine of theta appear, and r=sqrt(x^2+y^2).
178
179 However, the slimits parameter -- used to specify the points along
180 a ray where the transport integration starts and stops -- is
181 measured from the point (x,y,z) specified as a part of the
182 (cos,sin,y,z,x,r) ray coordinate. Thus, any change in the point
183 (x,y,z) on a ray must be accompanied by a corresponding change in
184 the slimits for that ray.
185
186 SEE ALSO: form_rays, best_rays, dirt_rays, get_s0, picture_rays
187 */
188 {
189 dummy= use_origins(0);
190
191 dim= dimsof(rays)(2);
192
193 if (dim==6) {
194 /* assure that r is consistent with x and y */
195 rays(6,..)= abs(rays(3,..), rays(5,..));
196 return rays;
197
198 } else if (dim==5) {
199 /* convert from best (x,y,z,theta,phi) to internal (cos,sin,y,z,x,r) */
200 x= rays(1,..);
201 y= rays(2,..);
202 z= rays(3,..);
203 theta= rays(4,..);
204 phi= rays(5,..);
205 cosp= cos(phi);
206 sinp= sin(phi);
207 sgn= double((cosp>=0.0)*2-1);
208 return form_rays([cos(theta), sgn*sin(theta), sgn*(y*cosp-x*sinp), z,
209 sgn*(x*cosp+y*sinp), abs(x, y)]);
210
211 } else if (dim==3) {
212 /* convert from DIRT/TDG (x,y,theta) to internal (cos,sin,y,z,x,r) */
213 x= rays(1,..);
214 y= rays(2,..);
215 theta= rays(3,..);
216 cost= cos(theta);
217 sint= sin(theta);
218 xcost= x*cost;
219 /* NOTE SIGN CHANGE IN THETA */
220 return form_rays([cost, -sint, y, x*sint, xcost, abs(xcost, y)]);
221 }
222
223 return [];
224 }
225
get_s0(rays)226 func get_s0(rays)
227 /* DOCUMENT get_s0(rays)
228 returns the s-coordinate of the point of closest approach of
229 the RAYS to the origin x=y=z=0. The length of the first dimension
230 of RAYS may be either 3, 5, or 6; this first dimension will not
231 be present in the result.
232
233 The s-coordinate represents distance along the ray, increasing in
234 the direction the ray moves. The 5 and 6 component ray coordinates
235 include a reference point (x,y,z) on the ray; s=0 at that point.
236 For the 3 component ray coordinate, get_s0 always returns 0.
237
238 SEE ALSO: best_rays, dirt_rays, internal_rays
239 */
240 {
241 dummy= use_origins(0);
242
243 dim= dimsof(rays)(2);
244
245 if (dim==3) {
246 return array(0.0, dimsof(rays(1,..)));
247
248 } else if (dim==5) {
249 x= rays(1,..);
250 y= rays(2,..);
251 z= rays(3,..);
252 theta= rays(4,..);
253 phi= rays(5,..);
254 cost= cos(theta);
255 sint= sin(theta);
256 x= x*cos(phi)+y*sin(phi);
257
258 } else if (dim==6) {
259 cost= rays(1,..);
260 sint= rays(2,..)
261 z= rays(4,..);
262 x= rays(5,..);
263 }
264
265 return -z*cost-x*sint;
266 }
267
picture_rays(xpict,ypict,ray,theta_up,phi_up)268 func picture_rays(xpict, ypict, ray, theta_up, phi_up)
269 /* DOCUMENT picture_rays(xpict, ypict, ray)
270 or picture_rays(xpict, ypict, ray, theta_up, phi_up)
271 returns 2-D array of rays, one at the center of each zone (which
272 represents a pixel here) of the mesh (XPICT, YPICT). The rays are
273 all parallel to the given RAY (a 3, 5, or 6 element vector). The
274 (XPICT, YPICT) coordinates are in the plane perpendicular to the rays,
275 with the origin XPICT=YPICT=0 at the given RAY.
276
277 If (THETA_UP, PHI_UP) are given, the +YPICT-axis will lie along the
278 projection of the (THETA_UP, PHI_UP) direction into the (XPICT, YPICT)
279 plane. The default (THETA_UP, PHI_UP) is (pi/2, pi/2) -- the +y-axis
280 -- unless (THETA, PHI) is the y-axis, in which case it is (pi/2, 0)
281 -- the +x-axis. This matches the DIRT/TDG ray coordinate convention
282 in the sense that if RAY is [0,0,theta], then
283 (zncen(XPICT),zncen(YPICT)) are the DIRT/TDG (x,y) coordinates for
284 the rays.
285
286 If XPICT and YPICT are imax-by-jmax, the returned array will have
287 dimensions 5-by-(imax-1)-by-(jmax-1). That is, "best" coordinates
288 are returned. The (x,y,z) of all of the returned rays will lie in
289 the plane perpendicular to the ray passing through the given central
290 RAY.
291
292 SEE ALSO: form_rays, best_rays, dirt_rays, internal_rays, area
293 */
294 {
295 dummy= use_origins(0);
296 ray= best_rays(ray);
297 x= ray(1);
298 y= ray(2);
299 z= ray(3);
300 theta= ray(4);
301 phi= ray(5);
302
303 /* form (rx,ry,rz) ray direction cosines */
304 sint= sin(theta);
305 rx= sint*cos(phi);
306 ry= sint*sin(phi);
307 rz= cos(theta);
308
309 /* form (ux,uy,uz) picture up direction cosines */
310 if (!is_void(theta_up)) {
311 sint= sin(theta_up);
312 ux= sint*cos(phi_up);
313 uy= sint*sin(phi_up);
314 uz= cos(theta_up);
315 } else {
316 if (rz || rx) { ux= 0.0; uy= 1.0; }
317 else { ux= 1.0; uy= 0.0; }
318 uz= 0.0;
319 }
320
321 /* form (px,py,pz) picture right direction cosines */
322 px= uy*rz - uz*ry;
323 py= uz*rx - ux*rz;
324 pz= ux*ry - uy*rx;
325 len= abs(px, py, pz);
326 px/= len;
327 py/= len;
328 pz/= len;
329
330 /* project (ux,uy,uz) perpendicular to (rx,ry,rz) */
331 ux= ry*pz - rz*py;
332 uy= rz*px - rx*pz;
333 uz= rx*py - ry*px;
334
335 /* get coordinates of zone centers in picture plane */
336 xp= xpict(zcen, zcen);
337 yp= ypict(zcen, zcen);
338
339 return form_rays([x+xp*px+yp*ux, y+xp*py+yp*uy, z+xp*pz+yp*uz,
340 theta, phi]);
341 }
342
plray(ray)343 func plray(ray)
344 /* DOCUMENT plray, ray
345 where RAY is a vector of length 5 representing [x, y, z, theta, phi]
346 -- a point (x,y,z) on the ray, and the ray direction (theta,phi)
347 relative to the z-axis. The ray hyperbola is plotted with the
348 z-axis horizontal. The portion of the ray plotted is determined
349 by the plray_lims command, which must be issued prior to the first
350 plray.
351 The 3 and 6 component ray formats are also accepted.
352 SEE ALSO: plray_lims, dirt_rays, internal_rays
353 */
354 {
355 dims= dimsof(ray);
356 if (dims(1)!=1) error, "ray must be a vector of length 5 (or 3 or 6)";
357 if (dims(2)!=5) ry= best_rays(ray);
358 else ry= ray;
359
360 dummy= use_origins(0);
361
362 x0= ry(1);
363 y0= ry(2);
364 z0= ry(3);
365 theta= ry(4);
366 phi= ry(5);
367
368 sint= sin(theta);
369 cost= cos(theta);
370 sinp= sin(phi);
371 cosp= cos(phi);
372
373 rtp= abs(x0*sinp - y0*cosp); /* turning point radius */
374
375 /**/ extern plray_zx, plray_zn, plray_rx;
376 /**/ extern plray_n; /* always odd! */
377
378 if (sint < 1.e-4) {
379 /* rays which are nearly parallel to the axis need special treatment */
380 z= span(plray_zn, plray_zx, 5);
381 if (cost<0.0) z= z(::-1);
382 r= (z-z0)*sint/cost;
383 r= abs(x0+r*cosp, y0+r*sinp);
384
385 } else {
386 /* form th == tan(angle)/abs(tan(theta)) where angle is evenly
387 spaced as th varies from -1 to 1 */
388 abstan= abs(sint)/(abs(cost)+1.e-20);
389 absth= atan(abstan);
390 fuzz= 1.e-6;
391 th= span(-absth, absth, plray_n);
392 th(1)= -(1.0-0.5*fuzz*fuzz);
393 th(2:-1)= tan(th(2:-1))/abstan;
394 th(0)= -th(1);
395
396 sqth= sqrt((1.0-th)*(1.0+th));
397 sqth(1)= sqth(0)= fuzz;
398
399 /* rays which nearly hit the axis need special treatment */
400 if (rtp < fuzz*plray_rx) rt= fuzz*plray_rx/sqth;
401 else rt= rtp/sqth;
402
403 zt= z0 + (rt*th - (x0*cosp-y0*sinp)*sign(sint))*sign(cost)/abstan;
404
405 /* clip (zt,rt) to (z,r) with z between plray_zn and plray_zx */
406 list= where(zt>plray_zn & zt<plray_zx);
407 if (min(abs(zt(dif))) > 0.0) {
408 if (cost>=0.0) { /* z increasing */
409 if (min(zt)<=plray_zn) {
410 z= [plray_zn];
411 r= interp(rt, zt, z);
412 } else {
413 z= [];
414 r= [];
415 }
416 grow, z, zt(list);
417 grow, r, rt(list);
418 if (max(zt)>=plray_zx) {
419 grow, z, [plray_zx];
420 grow, r, interp(rt, zt, [plray_zx]);
421 }
422 } else { /* z decreasing */
423 if (max(zt)>=plray_zx) {
424 z= [plray_zx];
425 r= interp(rt, zt, z);
426 } else {
427 z= [];
428 r= [];
429 }
430 grow, z, zt(list);
431 grow, r, rt(list);
432 if (min(zt)<=plray_zn) {
433 grow, z, [plray_zn];
434 grow, r, interp(rt, zt, [plray_zn]);
435 }
436 }
437 } else {
438 z= zt(list);
439 r= rt(list);
440 }
441 }
442
443 /* clip (z,r) to r < rx */
444 if (numberof(r) && numberof((list= where(r<plray_rx)))) {
445 if (numberof(list) < numberof(r)) {
446 rt= r;
447 zt= z;
448 i= rt(mnx); /* the interp operation seems to be safe even if many
449 of the rt are clustered at this minimum value */
450 if (min(list)<=i && max(rt(:i))>=plray_rx) {
451 r= [plray_rx];
452 z= interp(zt(:i), rt(:i), r);
453 } else {
454 r= [];
455 z= [];
456 }
457 grow, r, rt(list);
458 grow, z, zt(list);
459 if (max(list)>=i && max(rt(i:))>=plray_rx) {
460 grow, r, [plray_rx];
461 grow, z, interp(zt(i:), rt(i:), [plray_rx]);
462 }
463 }
464
465 /* finally, plot the result */
466 /* Note: Gist substitutes the marker for \001 if that is the
467 first character of the legend. */
468 legend=
469 swrite(format="\001: theta= %.4g, phi= %.4g, thru (%.4g, %.4g, %.4g)",
470 theta, phi, x0, y0, z0)(1);
471 plg, r, z, rays=1, legend=legend;
472
473 } else {
474 write, "WARNING (plray) ray not within limits (use plray_lims)";
475 }
476 }
477
plray_lims(zmin,zmax,rmax,rx)478 func plray_lims(zmin, zmax, rmax, rx)
479 /* DOCUMENT plray_lims, zmin, zmax, rmax
480 sets the (z,r) coordinate limits for the plray command. Subsequent
481 rays will be clipped to the region from z=ZMIN to z=ZMAX, and from
482 r=0 to r=RMAX.
483 SEE ALSO: plray
484 */
485 {
486 /**/ extern plray_zx, plray_zn, plray_rx;
487 if (!is_void(rx)) rmax= rx; /* allow for limits syntax as well */
488 if (zmin!=zmax && rmax>0.0) {
489 plray_zx= double(max(zmin, zmax));
490 plray_zn= double(min(zmin, zmax));
491 plray_rx= double(rmax);
492 limits, plray_zn, plray_zx, 0.0, plray_rx;
493 } else {
494 error, "zmin==zmax or rmax<=0.0 are not legal ray limits";
495 }
496 }
497
498 plray_zx= 1.0;
499 plray_zn= -1.0;
500 plray_rx= 1.0;
501 plray_n= 31; /* This should always be an odd number to ensure that
502 the turning point of the ray is always one of the
503 points plotted. */
504