raypath_prop_jd

 This function calculates a ray trajectory in a torus assuming straight propagation along unit vector g. 

 INPUTS

       - Rp: Plasma major radius (m) [1,1]
       - yx_L: launching toroidal major radius (m) [0...infty] [1,nb]
       - yy_L: launching vertical position (m) [-infty...infty] [1,nb]
       - yz_L: launching axial position (m) [-pi...pi] [1,nb]
       - yalpha_L: launching horizontal angle (with respect to x) [-pi...pi] [1,nb]
       - ybeta_L: launching vertical angle (with respect to y) [0...pi] [1,nb]
       - s: curvilinear positions along ray path (m) [1,ns]

 OUTPUTS
       
       - sysh: horizontal distance along trajectory (sh > 0) [ns,nb]
       - sysy: vertical distance along trajectory [ns,nb] 
       - syx: toroidal major radius [ns,nb]
       - syy: vertical position [ns,nb] 
       - syz: axial position [ns,nb]  
       - syalpha: horizontal angle (with respect to R) [-pi...pi] [ns,nb] 
       - sybeta: vertical angle (with respect to Z) [0...pi] [ns,nb] 
       - sygx: component of g along x [ns,nb] 
       - sygy: component of g along y [ns,nb] 
       - sygz: component of g along z [ns,nb] 

 by J. Decker (joan.decker@cea.fr) and Y. Peysson (yves.peysson@cea.fr)

0001 function [sysh,sysy,syx,syy,syz,syalpha,sybeta,sygx,sygy,sygz] = raypath_prop_jd(Rp,yx_L,yy_L,yz_L,yalpha_L,ybeta_L,s)
0002 %
0003 % This function calculates a ray trajectory in a torus assuming straight propagation along unit vector g.
0004 %
0005 % INPUTS
0006 %
0007 %       - Rp: Plasma major radius (m) [1,1]
0008 %       - yx_L: launching toroidal major radius (m) [0...infty] [1,nb]
0009 %       - yy_L: launching vertical position (m) [-infty...infty] [1,nb]
0010 %       - yz_L: launching axial position (m) [-pi...pi] [1,nb]
0011 %       - yalpha_L: launching horizontal angle (with respect to x) [-pi...pi] [1,nb]
0012 %       - ybeta_L: launching vertical angle (with respect to y) [0...pi] [1,nb]
0013 %       - s: curvilinear positions along ray path (m) [1,ns]
0014 %
0015 % OUTPUTS
0016 %
0017 %       - sysh: horizontal distance along trajectory (sh > 0) [ns,nb]
0018 %       - sysy: vertical distance along trajectory [ns,nb]
0019 %       - syx: toroidal major radius [ns,nb]
0020 %       - syy: vertical position [ns,nb]
0021 %       - syz: axial position [ns,nb]
0022 %       - syalpha: horizontal angle (with respect to R) [-pi...pi] [ns,nb]
0023 %       - sybeta: vertical angle (with respect to Z) [0...pi] [ns,nb]
0024 %       - sygx: component of g along x [ns,nb]
0025 %       - sygy: component of g along y [ns,nb]
0026 %       - sygz: component of g along z [ns,nb]
0027 %
0028 % by J. Decker (joan.decker@cea.fr) and Y. Peysson (yves.peysson@cea.fr)
0029 %
0030 if nargin < 7
0031     error('not enough arguments')
0032 end
0033 if min(yx_L) < -Rp
0034     error('R_L must be positive')
0035 end
0036 if min(yalpha_L) <= -pi || max(yalpha_L) > pi
0037     error('need -pi < alpha_L <= pi')
0038 end
0039 if min(ybeta_L) < 0 || max(ybeta_L) > pi
0040     error('need 0 <= beta_L <= pi')
0041 end    
0042 %
0043 ns = length(s);
0044 nb = max([length(yx_L),length(yy_L),length(yz_L),length(yalpha_L),length(ybeta_L)]);
0045 %
0046 if ~any(length(yx_L) == [1,nb]) || ~any(length(yy_L) == [1,nb]) || ~any(length(yz_L) == [1,nb]) || ~any(length(yalpha_L) == [1,nb]) || ~any(length(ybeta_L) == [1,nb])
0047     error('Launching inputs must have the same dimension')
0048 end
0049 %
0050 yx_L = repmat(yx_L,[1,nb/length(yx_L)]);
0051 yy_L = repmat(yy_L,[1,nb/length(yy_L)]);
0052 yz_L = repmat(yz_L,[1,nb/length(yz_L)]);
0053 yalpha_L = repmat(yalpha_L,[1,nb/length(yalpha_L)]);
0054 ybeta_L = repmat(ybeta_L,[1,nb/length(ybeta_L)]);
0055 %
0056 sysy = s'*cos(ybeta_L);
0057 sysh = s'*sin(ybeta_L);
0058 %
0059 syx_L = ones(ns,1)*yx_L;
0060 syy_L = ones(ns,1)*yy_L;
0061 syz_L = ones(ns,1)*yz_L;
0062 syalpha_L = ones(ns,1)*yalpha_L;
0063 sybeta = ones(ns,1)*ybeta_L;
0064 %
0065 syy = syy_L + sysy;
0066 %
0067 if isinf(Rp),
0068     syx = syx_L + sysh.*cos(syalpha_L);
0069     syz = syz_L + sysh.*sin(syalpha_L);
0070     syalpha = syalpha_L; 
0071 else
0072     syR_L = syx_L + Rp;
0073     syR = sqrt(syR_L.^2 + sysh.^2 + 2*syR_L.*sysh.*cos(syalpha_L));
0074     syx = syR - Rp;
0075     %
0076     syphi = acos((syR.^2 + syR_L.^2 - sysh.^2)./(2*syR.*syR_L)).*sign(syalpha_L);
0077     syphi = real(syphi);
0078     %
0079     syz = Rp*syphi + syz_L;
0080     %
0081     syalpha = syalpha_L - syphi;
0082 end
0083 %
0084 sygx = sin(sybeta).*cos(syalpha);
0085 sygy = cos(sybeta);
0086 sygz = sin(sybeta).*sin(syalpha);