Problem with plotting, i am getting nothing on the plot region

Hi everyone, here i attached my code in which i am not able to see plot between theta and Fe, and psi and Fe. I don't see any errors and any help is appreciated. Thanking you in advance.

[Ameer Hamza]

The main problem is that you are not storing the value of Fe inside the for-loop. It is overwritten after each iteration. However, the correct solution depends on your equations. Can you show them?

[satish kumar]

Thanks a lot Ameer Hamza, how can i save the my Fe array. i already wrote my formulas in there. i am attaching a figure regarding my formulas.

[satish kumar]

can i plot for those values of theta =(0:5*pi/180:pi/2), i only needed for those particular 19 values.

[Ameer Hamza]

The value of Fe depends on theta, phi, and psi. What is the variable in the x-axis? Why are you using a double for-loop? It is usually used to make a surface plot. What does the loop

for x = (1:6)

calculates?

[satish kumar]

 for x = (1:6) 

goes through each value of lambdat at a time and then goes through if statements.i want to plot theta values on x-axis and corresponding Fe values on y-axis. The same for psi also. Thanks for the response Ameer.

[Ameer Hamza]

But as you can see, for each value of theta in for-loop, there will be many values of phi and Fe. How do you suggest to plot one value against several values? Do you have an example plot you want to get at the end?

[satish kumar]

Yes, I do. i am attaching it here. y- axis contains Fe values in the png. Thanks again man.

[Ameer Hamza]

Yes, now it is a bit clear. But one more thing needs clarification. You see that for every value of theta and phi, there are 6 values of 'x', what about those 6 values? How to plot those 6 values?

[satish kumar]

Sorry, i forgot to put this in code, now i realized that. i want the sum all of these six values and make this as Fe1 or something. Then plot for theta vs Fe1 and psi vs Fe1. Thanks a lot Ameer for your efforts.

[satish kumar]

I just updated my code regarding this and i am attaching my code down here.

Answer 1

[Ameer Hamza]

Try this

clc
close all

Re=6378;         % Earth radius [km]
H=400;             % altitude
r=(Re+H);          % orbital radius
counter = 0;

phi=(0*pi/180);   % euler angles for a rotation
thetaV = (0:5*pi/180:pi/2);
psiV = (0:5*pi/180:pi/2);
Fe1 = zeros(numel(thetaV), numel(psiV));
for i = 1:numel(thetaV) % euler angles for a rotation
    theta = thetaV(i);
    
    for j = 1:numel(psiV)  % euler angles for a rotation
        psi = psiV(j);
        
        R=[cos(phi).*cos(psi)-sin(phi).*cos(theta).*sin(psi)   sin(phi).*cos(psi)+cos(phi).*cos(theta).*sin(psi)   sin(theta).*sin(psi)
            -cos(phi).*sin(psi)-sin(phi).*cos(theta).*cos(psi) -sin(phi).*sin(psi)+cos(phi).*cos(theta).*cos(psi) sin(theta).*cos(psi)
            sin(phi).*sin(theta)  -cos(phi).*sin(theta)  cos(theta)];
        
        v=[0 0 1]'; % local vertical in stabilized frame
        
        V_xyz=R*v;  % versor v components in body frame
        
        nxp=[1 0 0]; %vector normal to body frame with +x-face
        nxm=[-1 0 0]; %vector normal to body frame with -x-face
        nyp=[0 1 0] ;%vector normal to body frame with +y-face
        nym=[0 -1 0]; %vector normal to body frame with -y-face
        nzp=[0 0 1] ;%vector normal to body frame with +z-face
        nzm=[0 0 -1]; %vector normal to body frame with -z-face
        
        
        lambdatnxp=acos(nxp*V_xyz); % angle b/w normal nxp and local vertical in body frame
        lambdatnxm=acos(nxm*V_xyz); % angle b/w normal nxm and local vertical in body frame
        lambdatnyp=acos(nyp*V_xyz);% angle b/w normal nyp and local vertical in body frame
        lambdatnym=acos(nym*V_xyz); % angle b/w normal nym and local vertical in body frame
        lambdatnzp=acos(nzp*V_xyz); % angle b/w normal nzp and local vertical in body frame
        lambdatnzm=acos(nzm*V_xyz); % angle b/w normal nzm and local vertical in body frame
        
        lambdat=[lambdatnxp,lambdatnxm,lambdatnyp,lambdatnym,lambdatnzp,lambdatnzm];
        lambdam=asin(Re/(r));
        counter = counter +1;
        
        Fe_set = zeros(1,6);
        for x = (1:6)
            
            if lambdat(1,x) <= pi/2-lambdam
                Fe=(Re/r).^2.*cos(lambdat(1,x)); % first domain
                %      elseif pi/2 - lambdam < x && x < pi/2 + lambdam
            elseif lambdat(1,x)>pi/2+lambdam
                Fe=0;
                %Fe=(2/pi)*(0.25*pi-0.5*asin(sqrt((r/Re)^2-1)/((r/Re)*sin(lambdat)))+0.5*(Re/r)^2*(cos(lambdat)*acos(-cot(lambdat)*sqrt((r/Re)^2-1))+(-sqrt((1-(r/Re)^2*cos(lambdat)^2)*((r/Re)^2-1)))));   % second domain
            else
                Fe =(2/pi)*((pi/4)-0.5.*asin((sqrt(r/Re)^2-1)/((r/Re).*sin(lambdat(1,x))))+0.5*(Re/r)^2.*((cos(lambdat(1,x))).*acos(-cot(lambdat(1,x)).*(sqrt(r/Re).^2-1))-sqrt((1-(r/Re)^2.*cos(lambdat(1,x)).^2).*((r/Re).^2-1))));  % third domain
            end
            Fe_set(x+1) = Fe;
        end
        Fe1(i, j) = sum(Fe_set);
    end
    
end

theta=(0:5*pi/180:pi/2);
psi=(0:5*pi/180:pi/2);

% figure(1)
%  plot(lambdat*180/pi, Fe),'r';xlabel('Attitude angle \lambda_t (deg)'), ylabel('Earth view factor')
figure(2)
plot(theta*180/pi, Fe1),'b--*';xlabel('Attitude angle \theta (deg)'), ylabel('Earth view factor')
figure(3)
plot(psi*180/pi, Fe1),'o-.+';xlabel('Attitude angle \psi (deg)'), ylabel('Earth view factor')

According to the equations you have written, the graph comes out to be

[satish kumar]

Thanks a lot Ameer, this means a lot for me. But I am getting both the plots as similar, like the first graph you presented here. But here for the plot" psi vs Fe1" you got a slightly different graph compared to "theta vs Fe1".Let me know if there is any real difference between those two plots for you. Apart from that, i am very much satisfied with the solution. Once again thanks a lot Ameer.

[Ameer Hamza]

I am glad to be of help. The similarity is not caused by the plot() function. I suspect there is some minor mistake in your equations, which might be giving wrong outputs. You may carefully analyze your equations to see if there is an unexpected output at any step.