How to plot this equation in MATLAB? Keep getting matrix dimension errors

Question

The code I have is:

T=[0:0.1:24];
omega=((12-T)/24)*360;
alpha = [0:1:90];
latitude=35;
delta=[-23.45:5:23.45];
sind(alpha)=sind(delta).*sind(latitude)+cosd(delta).*cosd(latitude).*cosd(omega)
cosd(phi)=(sind(alpha).*sind(latitude)-cosd(delta))./(cosd(alpha).*cosd(latitude))

alpha represents my y-axis and is an angle between 0 and 90 degrees. phi represents my x-axis and is an angle between say -120 to +120. The overall result should look something like a half sine-wave.

Whenever I try to input that last line I get the error stating inner matrix dimensions must agree. So I tried to use reshape on my matrices for those variables I defined so that they work. But then I get '??? ??? Subscript indices must either be real positive integers or logicals.'

It seems very tedious to have to reshape my matrix every time I define a new set of variables in order to use them with an equation. Those variables are used for defining my axis range, is there a better way I can lay them out or an automatic command that will make sure they work every time?

I want to plot alpha and phi as a graph using something like

plot(alpha,phi) 

but can't get past those errors? can't I just use a command that says something like define x-axis [0:90], define y-axis [-120:120] or something? I have spent far too much time on this problem and can't find a solution. I just want to plot the graph. Somebody please help! thanks.

Thanks

Answer 1

Here is a start for you:

% Latitude
L=35;

% Hour Angle
h = [-12:5/60:12];
w = 15*h;

% Initialize and clear plot window
figure(1); clf;

% Plot one day per month for half of the year
for N = 1:30:365/2
    % Declination
    d = 23.45*sind(360*(284+N)/365);

    % Sun Height
    alpha = asind(sind(L)*sind(d) + cosd(L)*cosd(d)*cosd(w));

    % Solar Azimuth 
    x = ( sind(alpha)*sind(L)-sind(d) )./( cosd(alpha)*cosd(L) );
    y = cosd(d)*sind(w)./cosd(alpha);
    phi = real(atan2d(y,x));

    % Plot
    plot(phi,alpha); hold on;
end

hold off;
grid on;
axis([-180, 180, 0, 90]);
xlabel('Solar Azimuth')
ylabel('Solar Elevation')

The function asind is inherently limited to return values in the range of -90 to 90. That means that you will not get a plot that spans over 240 degrees like the one you linked to. In order for the plot to not have the inflection at +/- 90 degrees, you will need to find a way to infer the quadrant. Update: I added the cosine term to get an angle using atan2d.

Hopefully this will be enough to give you an idea of how to use Matlab to get the kind of result that you are after.

Answer 2

Element-wise multiplication (.*) and division ./, not the matrix versions:

sind(alpha)=sind(delta).*sind(latitude)+cosd(delta).*cosd(latitude).*cosd(omega)
cosd(phi)=(sind(alpha).*sind(latitude)-cosd(delta))./(cosd(alpha).*cosd(latitude))

But a bigger problem is that you can't index with decimal or non-positive values. Go back to the code before you got the subscript indices error. The indices "must either be real positive integers or logicals".

Answer 3

I think you're a bit confused about MATLAB sintax (as nispio mentioned). Maybe you want to do something like this?

T=[0:0.1:24];
omega=((12-T)/24)*360;
alpha = [0:1:90];
latitude=35;
delta=[-23.45:5:23.45];
[o1,d1]=meshgrid(omega,delta);
[a2,d2]=meshgrid(alpha,delta);
var1=sind(d1(:)).*sind(latitude)+cosd(d1(:)).*cosd(latitude).*cosd(o1(:));
var2=(sind(a2(:)).*sind(latitude)-cosd(d2(:)))./(cosd(a2(:)).*cosd(latitude));
plot(var1);hold on;plot(var2);

If not, you should post the algorithm or pseudocode