% line of sight - NEW 15/09/20

% instead of pre-generating plane of view, incorporate into loop through
% all points in the 3D-foreshock emissions
% video = VideoWriter('test3.avi');
% video.FrameRate = 1;
% open(video);
RR = size(jF,1);
xR = size(jF,2);
dR = 2*Rrange/(RR-1);
dx = xrange/xR;

% RR = size(jF_nlt,1);
% xR = size(jF_nlt,2);
% dR = 2*Rrange/(RR-1);
% dx = xrange/xR;
figure()
tiledlayout(3,3)
% 
for omega = 0:-10:-80
  nexttile
% omega = 0;
% omega = 360-theta_deg;
RE = 6.3780e6;
dist = 200*RE;
% RE = 1;
% O = [100000*Rp*sind(omega),-100000*Rp*cosd(omega),0];
O = [dist*cosd(omega),dist*sind(omega),0];

% O = [1e6,0,0];
% O = [200*RE,0,0];
% O = [20*RE,20*RE,0];
% O = [Xt - stheta*200*Rp, Yt - ctheta*200*Rp, 0];
% O = [200*Rp,200*Rp,0];
% O = [0,-200*Rp,0];
% O = [1,0,200*Rp];

% for i = 2:102
%     jF(:,:,i) = jF(:,:,1);
%     jH(:,:,i) = jH(:,:,1);
%     jF_lmc(:,:,i) = jF_lmc(:,:,1);
% end
% mirror emissions to Z- space
    Znmax = size(jF,3);
%     Znmax = size(jF_nlt,3);
%     Znmax = 102;


if exist('jF','var') == 1        
    jF_3D = flip(jF(:,:,1:Znmax),3);
    jF_3D(:,:,Znmax+1:2*Znmax-1) = jF(:,:,2:Znmax);
    
    jH_3D = flip(jH(:,:,1:Znmax),3);
    jH_3D(:,:,Znmax+1:2*Znmax-1) = jH(:,:,2:Znmax);
    
    jLMC_3D = flip(jF_lmc(:,:,1:Znmax),3);
    jLMC_3D(:,:,Znmax+1:2*Znmax-1) = jF_lmc(:,:,2:Znmax);
else
    

    jF_3D = flip(jF_nlt(:,:,1:Znmax),3);
    jF_3D(:,:,Znmax+1:2*Znmax-1) = jF_nlt(:,:,2:Znmax);
    jH_3D = flip(jH(:,:,1:Znmax),3);
    jH_3D(:,:,Znmax+1:2*Znmax-1) = jH(:,:,2:Znmax);    
    jLMC_3D = flip(jF_lmc(:,:,1:Znmax),3);
    jLMC_3D(:,:,Znmax+1:2*Znmax-1) = jF_lmc(:,:,2:Znmax);
end

% for i = 1:size(jF_3D,3)
%     jF_3D(:,:,i) = jF_3D(:,:,Znmax);
%     jH_3D(:,:,i) = jH_3D(:,:,Znmax);
%     jLMC_3D(:,:,i) = jLMC_3D(:,:,Znmax);
% end

Zrange = (Znmax-1)*dZ;
Znmax = size(jF_3D,3);

% distance observer to object
O2 = sqrt(O(1)^2+O(2)^2+O(3)^2);
F_S_c = zeros(RR,xR,Znmax,3);

% Euler angle sequence for coordinate system
% transformation
yaw = atand(O(2)/O(1));
pitch = atand(O(3)/sqrt(O(1)^2+O(2)^2));
roll = 0;
% Transformation Matrix
C = [cosd(pitch) 0 -sind(pitch); 0 1 0; sind(pitch) 0 cosd(pitch)]*[cosd(yaw) sind(yaw) 0;-sind(yaw) cosd(yaw) 0;0 0 1];

%counter
ZRx = 0;
S_YZ = zeros(Znmax*RR*xR,5);
% figure()
% hold on
F_tot = 0;
H_tot = 0;
LMC_tot = 0;
V = xrange/xR*(2*Rrange/RR)*dZ;
for Zn = 1:Znmax    
% for Zn = (Znmax+1)/2+2
    Z = -Zrange + (Zn-1)*dZ;
%     Z = -Zrange + 99*dZ;    
    
    for Rn = 1:RR
%     for Rn = 125:875
        R = -Rrange + (Rn-1)*dR;
        
        for xn = 1:xR
%         for xn = 1:750
            ZRx = ZRx+1;
            x = xn*dx;            
            Xk = Xt + x*stheta + R*ctheta - b*Z^2;
            Y = Yt + x*ctheta - R*stheta;
%             Xk = Xt - x*stheta - R*ctheta -b*Z^2;
%             Y = Yt + x*ctheta - R*stheta;
%             if Y == Yt
%                 z = z;
%             end
            % check where the line of sight crosses the shock
            D = sqrt((Xk-O(1))^2+(Y-O(2))^2+(Z-O(3))^2);
            los = [Xk-O(1),Y-O(2),Z-O(3)]/D;
            
            % input O+t*los in parabola equation =>
            % Ox+t*losx - as + b*(Oy+t*losy)^2 + b*(Oy+t*losz)^2 = 0
            % c = Ox - as + b*(Oy^2+Z^2)
            % b = losx + b*(2*Oy*losy+2*Oz*losz)
            % a = b*(losy^2 + losz^2)
            
            ta = b*(los(2)^2 + los(3)^2);
            tb = b*(2*O(2)*los(2) + 2*O(3)*los(3)) + los(1);
            tc = O(1) - as + b*(O(2)^2 + O(3)^2);
            t = (-tb + sqrt(tb^2 - 4*ta*tc))/(2*ta);
            
            % calculate observed flux
            F_Rxz = abs(jF_3D(Rn,xn,Zn))/D^2*coneF*V;
            F_Rnxn(Rn,xn)=F_Rxz;
            F_tot = F_tot + F_Rxz;
            H_Rxz = abs(jH_3D(Rn,xn,Zn))/D^2*coneH*V;
            H_tot = H_tot + H_Rxz;
            LMC_Rxz = abs(jLMC_3D(Rn,xn,Zn))/D^2*coneF*V;
            LMC_tot = LMC_tot + LMC_Rxz;
            % find location of emission on 2D plane            
            % plane for which ObsOrg is the normal vector
            % Obs(1)*x + Obs(2)*y + Obs(3)*z = 0
            % intersection of ObsEms with plane
            % S = O(1)*(O(1)+s*losx) + Obs(2)*(Obs(2)+s*losy) +
            % Obs(3)*(Obs(3)+s*losz) = 0
            s = -O*O'/(O*los');
            S = O + s*los;
            % new location of S on Y'Z' plane
            S_c = C*S';            
            
            if imag(t) == 0
%                 if S_c(1) == 0
                if t > D
%                     if mod(xn,10) == 0 && mod(Rn,10) == 0 && mod(Zn,10) == 0
%                         scatter(S_c(2),S_c(3),'o')
%                     end
                    S_YZ(ZRx,:) = [S_c(2),S_c(3),F_Rxz,H_Rxz,LMC_Rxz];
%                     else
%                         scatter(S_c(2),S_c(3),'o','filled')
                end
%                 end                                                      
            else
%                 if mod(xn,10) == 0 && mod(Rn,10) == 0 && mod(Zn,10) == 0
%                     scatter(S_c(2),S_c(3),'o')
%                 end
%                 S_YZ(ZRx,:) = [S_c(2),S_c(3),F_Rxz,H_Rxz,LMC_Rxz];
            end
%             F_S_c(Rn,xn,Zn,:) = [S_c(2),S_c(3),F_Rxz];                                     
        end        
    end
end
F_tot
H_tot
LMC_tot
dyz = Rp;
dy = 0.7*Rp;
dz = 1.1*dZ;
% dy = 0.1*Rp;
% dz = 0.25*dZ;
Y_min = min(S_YZ(:,1));
Y_max = max(S_YZ(:,1));
Z_min = min(S_YZ(:,2));
Z_max = max(S_YZ(:,2));
Y_dash = Y_min:dy:Y_max;
Z_dash = Z_min:dz:Z_max;

F_YZ = zeros(length(Y_dash)+1,length(Z_dash)+1);
H_YZ = zeros(length(Y_dash)+1,length(Z_dash)+1);
LMC_YZ = zeros(length(Y_dash)+1,length(Z_dash)+1);

for ZRx = 1:size(S_YZ,1)
    Y_n = round((S_YZ(ZRx,1)-Y_min)/dy)+1;
    Z_n = round((S_YZ(ZRx,2)-Z_min)/dz)+1;
    if F_YZ(Y_n,Z_n) ~= 0
        F_YZ(Y_n,Z_n) = F_YZ(Y_n,Z_n) + S_YZ(ZRx,3);        
    else
        F_YZ(Y_n,Z_n) = S_YZ(ZRx,3);
    end
    if H_YZ(Y_n,Z_n) ~= 0
        H_YZ(Y_n,Z_n) = H_YZ(Y_n,Z_n) + S_YZ(ZRx,4);        
    else
        H_YZ(Y_n,Z_n) = S_YZ(ZRx,4);
    end
    if LMC_YZ(Y_n,Z_n) ~= 0
        LMC_YZ(Y_n,Z_n) = LMC_YZ(Y_n,Z_n) + S_YZ(ZRx,5);        
    else
        LMC_YZ(Y_n,Z_n) = S_YZ(ZRx,5);
    end
end

%% Filling out missing planes
% for i = 1:2
% tempFYZ = zeros(size(F_YZ));
% for y = 2:size(F_YZ,1)-1
%   for z = 2:size(F_YZ,2)-1
%     if F_YZ(y,z) == 0
%       Ycomp = (F_YZ(y-1,z)+F_YZ(y+1,z));
%       Zcomp = abs(F_YZ(y,z-1))+abs(F_YZ(y,z+1));
%       if F_YZ(y-1,z) ~= 0 && F_YZ(y+1,z) ~= 0
%         tempFYZ(y,z) = Ycomp/2;
%       elseif F_YZ(y-1,z) == 0
%         tempFYZ(y,z) = F_YZ(y+1,z);
%       elseif F_YZ(y+1,z) == 0
%         tempFYZ(y,z) = F_YZ(y-1,z);  
%       end
%     else
%       tempFYZ(y,z) = F_YZ(y,z);
%     end    
%   end
% end
% F_YZ = tempFYZ;  
% tempFYZ = zeros(size(F_YZ));
% for y = 2:size(F_YZ,1)-1
%   for z = 2:size(F_YZ,2)-1
%     if F_YZ(y,z) == 0
%       Ycomp = (F_YZ(y-1,z)+F_YZ(y+1,z));
%       Zcomp = abs(F_YZ(y,z-1))+abs(F_YZ(y,z+1));
%       if F_YZ(y,z-1) ~= 0 && F_YZ(y,z+1) ~= 0
% %       if F_YZ(y+1,z) == 0
%         tempFYZ(y,z) = Zcomp/2;
%       elseif F_YZ(y,z-1) == 0
%         tempFYZ(y,z) = F_YZ(y,z+1);
%       elseif F_YZ(y,z+1) == 0
%         tempFYZ(y,z) = F_YZ(y,z-1);  
%       end
% %       elseif F_YZ(y,z+1) == 0
% %         F_YZ(y,z) = Ycomp/2;
% %       end
% %       F_YZ(y,z) = (Ycomp+Zcomp)/2;
%     else
%       tempFYZ(y,z) = F_YZ(y,z);
%     end    
%   end
% end
% F_YZ = tempFYZ;
% end

% %% Smoothing out
% fact = 0.6;
% % for i=1:2
%   F_YZ_n = zeros(size(F_YZ));
%   for y = 2:size(F_YZ,1)-1
%     for z = 2:size(F_YZ,2)-1
%       if F_YZ(y,z) ~= 0
%         Ycomp = (F_YZ(y-1,z)+F_YZ(y+1,z));
%         Zcomp = (F_YZ(y,z-1)+F_YZ(y,z+1));
%         if F_YZ(y,z) > fact*Ycomp && F_YZ(y,z) > fact*Zcomp
%           F_YZ_n(y,z) = (Ycomp+Zcomp)/4;
%         elseif F_YZ(y,z) > fact*Ycomp
%           F_YZ_n(y,z) = Ycomp/2;
%         elseif F_YZ(y,z) > fact*Zcomp
%           F_YZ_n(y,z) = Zcomp/2;
%         else
%           F_YZ_n(y,z) = F_YZ(y,z);
%         end
%       end
%     end
%   end
%   F_YZ = F_YZ_n;
% % end

% figure()
% imagesc(Y_dash/Rp,Z_dash/Rp,log10(F_YZ'),[-23,-20])
% colors = [jet];
% colormap(colors)
% set(gca,'YDir','normal','FontSize',18,'FontWeight','bold');
% xlabel('Y [R_M]')
% ylabel('Z [R_M]')
% h = colorbar;
% ylabel(h,'Flux [W m^{-2} Hz^{-1}]')
% axis equal
% set(gca,'Color',[0 0 0.5])
% % 
% figure()
% imagesc(Y_dash/Rp,Z_dash/Rp,log10(H_YZ'),[-23,-20])
% colors = [jet];
% colormap(colors)
% set(gca,'YDir','normal','FontSize',18,'FontWeight','bold');
% xlabel('Y [R_M]')
% ylabel('Z [R_M]')
% h = colorbar;
% ylabel(h,'Flux [W m^{-2} Hz^{-1}]')
% axis equal
% set(gca,'Color',[0 0 0.5])
% % 
% figure()
% imagesc(Y_dash/Rp,Z_dash/Rp,log10(LMC_YZ'),[-23,-20])
% colors = [jet];
% colormap(colors)
% set(gca,'YDir','normal','FontSize',18,'FontWeight','bold');
% xlabel('Y [R_M]')
% ylabel('Z [R_M]')
% h = colorbar;
% ylabel(h,'Flux [W m^{-2} Hz^{-1}]')
% axis equal
% set(gca,'Color',[0 0 0.5])
% % writeVideo(video,getframe(gcf));
% % 
imagesc(Y_dash/Rp,Z_dash/Rp,log10(H_YZ'),[-23,-20])
colors = [jet];
colormap(colors)
% set(gca,'YDir','normal','FontSize',18,'FontWeight','bold');

end
% % close(video);
