filling a matrix with information regarding angles using atan-function

Dear reader,

I am currently working on a line of sight effect model, using a circular plasma cross-section. The plasma is considered to be symmetric. As I would like to add effects just outside of the plasma, I use cartesian coordinates and generate a NxN matrix containing radial information:

clc
clear all
close all

N=500;      %number of cells used (matrix will be NxN)
R=4;        %length of x and y axis [cm]
r=3;        %plasma radius [cm]
R_LOS=-1;    %LOS position [cm]

%generating x and y locations:
x_axis=-R:2*R/(N-1):R;
y_axis=-R:2*R/(N-1):R;

for i=1:1:N;
    position_matrix(i,:)=sqrt(y_axis(i)^2+x_axis.^2);
end

this part just calculates the r-values with respect to the centre of the plasma. I also need to calculate the angles of each point with respect to the observer (which will just be the x-component of a velocity vector). I do know the velocity as function of radius, so i can generate a NxN matrix containing the information, using the ''position_matrix''. So i only need the angles. I thought the atan(y/x) the best way to deal with this, however this only helps me to obtain the angles in 1 quarter of the plasma. So i am currently doing the following to obtain the other 3/4 of my matrix:

        for i=1:1:N/2;
            for j=1:1:N/2;
                theta(i,j)=abs(atan(x_axis(i)/y_axis(j)));
            end
        end

            theta=fliplr(theta);
            figure(3)
            subplot(2,2,1)
            imagesc(theta);colorbar;
            subplot(2,2,2); imagesc(rot90(theta)+pi/2);colorbar;
                theta_2=[rot90(theta)+pi/2 theta];
            subplot(2,2,3); imagesc(theta_2); colorbar;
                theta_3=[theta_2  ; rot90(rot90(theta_2))+pi];
            subplot(2,2,4); imagesc(theta_3); colorbar; 

            % i use different theta_# to keep track of the step and its
            % affect. 

As you can see i am currently using 2 for loops, which isnt ideal (generally loops take long to go through) and i am forced to work with N-values that are a multiple of 2, or I wont have identical matrix sizes (theta would be an N-1xN-1 matrix for uneven values of N). So this made me wonder if there is a better (and faster) way to deal with this kind of a problem?

For those wondering what the other part of the model looks like (warning I have only just started, hence this is just a really simplistic model) here are the other lines of code with a simple test profile (temperature) to check if it works:

           %lets create a matrix with temperature information:
T_matrix=3.*exp(-position_matrix.^2/(r/2));


%replacing all r values greater than 3 (outside our plasma) by zero)
plasma_matrix=position_matrix;
plasma_matrix(plasma_matrix>r)=0;
            %please note that we can no longer use this matrix for profile
            %calculations. this is just for the plot.
 
 
figure(1)
subplot(1,2,1)
imagesc(x_axis,y_axis,plasma_matrix)
colorbar
subplot(1,2,2)
imagesc(x_axis,y_axis,T_matrix)
colorbar

%now we can check what happens if we put a LOS through it:
%we need to find the location of the y_axis values that is closest to
%R_LOS:
    [minValue, closestIndex] = min(abs(y_axis - R_LOS.'));
    closestValue = y_axis(closestIndex);
        %lets see how our LOS looks like:
        figure(1)
        subplot(1,2,1); hold on;
        plot(x_axis,ones(N,1).*closestValue,'-r'); hold off ; title('plasma cross-section');ylabel('radius [cm]');xlabel('radius [cm]');
        subplot(1,2,2); hold on;
        plot(x_axis,ones(N,1).*closestValue,'-r'); hold off; title('T_e [eV]');ylabel('radius [cm]');xlabel('radius [cm]');

        %now we can check our T_e profile along this line:
        figure(2)
        plot(x_axis,T_matrix(:, closestIndex),'-k'); grid on;
            hold on; plot(x_axis,T_matrix(:, closestIndex),'*r'); hold off ;xlabel('radius [cm]'); ylabel('T_e [eV]'); title('T_e along LOS'); 

thanks in advance!

btw if you have any remarks regarding those types of models in general (tips and tracks) I am more than happy to read through them!

Answer 1

[Matt J]

[X,Y]=ndgrid( -R:2*R/(N-1):R );

[theta, position_matrix]=cart2pol(X,Y);

[Jesse Verstappen]

Dear Matt J,

thanks for your reply! That seems to asnwer both my questions in only 2 lines of codes. Thanks for taking your time to help me out.

yours,

Jesse