Initial commit of raycon.

Author: arichar6

adjust_disp_m.m

@@ -0,0 +1,49 @@
+function y_ant = adjust_disp_m(y,m)
+    global plasma %cnst
+    
+    % ----- Adjusted position & wave vector in cylindrical coordinates -----------
+    r=y(1); z=y(2); 
+    kr=y(3); kz=y(4); % starting wave vector
+    rho=sqrt((r-plasma.r0).^2+z.^2);
+    theta=atan2(plasma.r0-r,-z)+pi/2;
+
+    krho_old=-sqrt(kr^2 + kz^2 -(m/rho)^2); % assume wave starts out going in
+    %krho=+sqrt(kr^2 + kz^2 -(m/rho)^2); 
+
+    % Want to adjust kr and kz so that, for given m, k_phi, the initial
+    % condtions are on the dispersion surface.
+    % m = function argument = rho*k_theta
+    
+    %kf = plasma.kant(2);
+
+    %y0=[r,z,kr,(kr*sin(theta)+(m/rho))/cos(theta)];
+    
+    % kz = (kr*sin(theta)+(m/rho))/cos(theta)
+%     f = @(krho) dispersion([r,z,...
+%         krho*cos(theta)-(m/rho)*sin(theta),...
+%         -(krho*sin(theta)+(m/rho)*cos(theta))],'Dsp');
+%     krho = fzero(f,krho);
+    
+    f = @(krho) disp_eig([r,z,...
+        krho*cos(theta)-(m/rho)*sin(theta),...
+        -(krho*sin(theta)+(m/rho)*cos(theta)), reshape(eye(4),1,16)].','Dsp');
+%     f = @(krho) disp_eig([r,z,...
+%         krho*cos(theta)-(m/rho)*sin(theta),...
+%         -(krho*sin(theta)+(m/rho)*cos(theta)), zeros(1,16)].','Dsp');
+    krho_new = fzero(f,krho_old);
+    
+    %kz = (kr*sin(theta)+(m/rho))/cos(theta);
+    
+    y_ant = [r,z,krho_new*cos(theta)-(m/rho)*sin(theta),...
+             -(krho_new*sin(theta)+(m/rho)*cos(theta))].';
+    
+    % ----- Wave vector and refraction index -------------------------------------
+%     kn=kr.*ener +kf.*enef +kz.*enez;
+%     kb=kr.*eber +kf.*ebef +kz.*ebez;
+%     kp=kr.*eper +kf.*epef +kz.*epez;
+%     T =[ener enef enez;...
+%         eber ebef ebez;...
+%         eper epef epez];
+%     k_cart = T^(-1) * [kn;kb;kp];
+end
+

calcFlux.m

@@ -0,0 +1,16 @@
+function plasma = calcFlux(plasma)
+
+    %% Calculate the flux surface coordinates ('mesh' and 'boundary')
+    s=0.0001:(.999/(plasma.NS-1)):1.;                 %  radial mesh
+    theta=0:2*pi/plasma.NT:2*pi;                      %  poloidal mesh
+    plasma.r = zeros(plasma.NS,plasma.NT+1);
+    plasma.z = zeros(plasma.NS,plasma.NT+1);
+    for kk=1:plasma.NS
+        si=s(kk).*ones(size(theta));             %  flux surfaces
+        [rho,plasma.r(kk,:),plasma.z(kk,:)]=mapFlux(si,theta); %    coordinates
+    end;
+
+    plasma.Rho = sqrt((plasma.r-plasma.r0).^2 +plasma.z.^2);
+    plasma.Theta = atan2(-plasma.z,(plasma.r-plasma.r0));
+
+    

cgamma.m

@@ -0,0 +1,63 @@
+function out = cgamma(z)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%  Complex gamma function for -10<Re(z)<10, 10<Im(z)<10
+%
+%  Tested 12 digits accuracy for
+%   cgamma([0 .5+.5i -.5+.5i -.5-.5i .5-.5i 1 1+i i -1+i -1 -1-i -i 1-i]')
+%
+%  A. JAUN, Numerical Analysis, KTH, 100 44 Stockholm, Sweden
+%  A.N. KAUFMAN, Lawrence Berkeley Laboratory, Berkeley, CA 94720, USA
+%  E.R. TRACY, College of William & Mary, Williamsburg, VA 23187-8795, USA
+%
+%  Documented under "http://www.nada.kth.se/~jaun"
+%
+%  (C) Version 7.0,  14-Aug-2006. All Rights Reserved.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+x=real(z); y=imag(z); n=prod(size(x));
+%
+% ----- Map to domain of asymptotic expansion
+%
+for k=1:n
+  if (y(k)<0) 
+    z(k)=conj(z(k));
+  end
+end
+%
+lnCgamma=zeros(size(x));
+for k=1:n
+  for m=1:floor(9-x)
+   lnCgamma(k)=lnCgamma(k)-log(z(k)); z(k)=z(k)+1;
+  end
+end
+%
+% ----- First compute log(cgamma) from expansion for 9<Re(z)<10, 0<Im(z)<10
+%
+zm =1./z;   zm1 =zm; zm2=zm1.*zm1;        % Avoid power function
+zm=zm.*zm2; zm3 =zm;
+zm=zm.*zm2; zm5 =zm;
+zm=zm.*zm2; zm7 =zm;
+zm=zm.*zm2; zm9 =zm;
+zm=zm.*zm2; zm11=zm;
+zm=zm.*zm2; zm13=zm;
+zm=zm.*zm2; zm15=zm;
+%
+lnCgamma=lnCgamma -zm15*(3617./122400);   % Add series in reverse order
+lnCgamma=lnCgamma +zm13/156;
+lnCgamma=lnCgamma -zm11*(691./360360);
+lnCgamma=lnCgamma +zm9 /1188;
+lnCgamma=lnCgamma -zm7 /1680;
+lnCgamma=lnCgamma +zm5 /1260;
+lnCgamma=lnCgamma -zm3 /360;
+lnCgamma=lnCgamma +zm1 /12 +(z-0.5).*log(z) -z +0.5*log(2*pi);
+%
+for k=1:n
+%  if (y(k)>0)                             % Don't understand why not <0
+  if (y(k)<0)                             % modified by Steve Richardson
+    lnCgamma(k)=conj(lnCgamma(k));
+  end
+end
+%
+% ----- Take exponential
+%
+out = exp(lnCgamma);