function [Z,d,Ro] = gyrotronRadiusRTM(R,theta,m)
%[z,d,Ro] = gyrotronRadiusRTM(R,theta,m)
% Based on the Radius and taper angle, the distance 
% between changes in angle is calculated.  This is 
% sort of the inverse of gyrotronRadiusTDM, which 
% calculates radius based on distance and angle and 
% then interpolates to create a mesh of points for 
% each section.
% The inputs are 'Radius', 'Theta' (angle) and 'Mesh',
% hence the 'RTM' in the function name.
%   Note: for the case of indeterminate distance
% (theta = 0), INF is returned as the distance.  Also,
% If the radius is increasing, theta should be positive.
%
% Inputs:
%   'R': [distance] the corresponding radius vector 
%   defining the points where theta changes.
%   'theta': [radians] a vector of theta angles 
%   that describe the increase in radius with increasing 
%   axial distance.
%   'm': [] the mesh size, that is, the number of points
%   of interpolation for each section.
%
% Output:
%   'z': [distance] the cumulative axial distance vector.
%   'd': [distance] the axial length vector containing 
%   the lengths of each section.
%   'Ro': [distance] radius at each interpolated mesh point.
%
% Example:
%   R     = [1 2 1.5 1];
%   theta = [2 1 -1 -2]*pi/180;
%   m     = [1 2 2   1]*10;
%   [z,d,Ro] = gyrotronRadiusRTM(R,theta,m);
%   plot(z,Ro,'b.',z,Ro,'g');
%
% Last edit by Eunmi Choi and Colin Joye 7/15/03
Z = [];
Ro= [];

% 'g' is the radius offset for each new section to make it continuous
g = [0 R];%.*sign(theta)];
g = g(1:(end-1));

% check to make sure if radius is increasing that theta>0 and vv.  Also
% make sure theta ~= 0.
A = R-[0 R(1:(end-1))];
B = tan(theta);
if( any( B./A <= 0 ) ),
    disp(['gyrotronRadiusRTM.m: Radius and Theta values are inconsistent' ...
        'or theta is zero.']);
    d=[]; 
    break;
end
d = abs( A ./ B );
d = [0 d];

% loop to generate outputs:
%   'inc' is the step size based on the input mesh number
%   'z' is the local values to generate the incremental radius
%   'Ro' is the new interpolated radius vector.
for k = 1:length(g),
  inc = d(k+1)/m(k);
  z   = [0 : inc : (d(k+1)-inc)];
  Ro  = [Ro z*tan(theta(k))+g(k)];
  Z   = [Z z+sum(d(1:k))];
end

% ------------ End gyrotronRadiusRT.m -------------