clear all
close all 
clc

tic % start timer


fullFileName     = 'C:\Users\Dell\Desktop\ayushman\LOLA_NAC\mosaic_files\NAC_ROI_VIKRAM__LOC_P708S0237_5M.IMG';

fid              = fopen(fullFileName);
Data             = fread(fid,'*float32');

samples_per_line = 9596;
lines            = 8637;

arr              = Data(9597:end); % In label it is mentioned that the data starts from 2nd record 
                                   % and first record(till 9596 samples) is label content. 
reshaped_array         = reshape(arr,samples_per_line,lines); 
final = reshaped_array;  % get 9596 x 8637 array

%% Preallocating arrays with 0 value to improve time taken for code to run
Lon = zeros(size(reshaped_array,1),size(reshaped_array,2));
Lat = zeros(size(reshaped_array,1),size(reshaped_array,2));
C   = zeros(size(reshaped_array,1),size(reshaped_array,2));
x   = zeros(size(reshaped_array,1),size(reshaped_array,2));
y   = zeros(size(reshaped_array,1),size(reshaped_array,2));

scale_in_km = 5/1000;
R_in_km               = 1737.4;

LonP = 23.7;
LatP = -70.8; % OR -90 

S0 = 4805.5;
L0 = 118749.5;

lines            = size(reshaped_array,2);
samples_per_line = size(reshaped_array,1);

%% using polar stereographic projection

for sample = 1:samples_per_line
    for line = 1:lines

        x(sample,line) = (sample-S0-1)*scale_in_km;   
        y(sample,line) = (1-L0-line)*scale_in_km; 
        ysquare = (y(sample,line))^2;   % temporary variable
        xsquare = (x(sample,line))^2;   % temporary variable
        P = sqrt(xsquare + ysquare);        
        t = P/(2*R_in_km);              % t is a temporary variable
        C(sample,line) = 2*atand(t);
        
%%   uncomment below line if using equation 20-17 for longitude calculation

   %     Lon(sample,line) = LonP + atand((x(sample,line))/(y(sample,line))); 
%% Uncomment below line if using Equation 20-15 for latitude calculation       

        a1 = ((x(sample,line))*sind(C(sample,line)));
        a2 = (P)*(cosd(LatP))*cosd(C(sample,line));
        a3 = (y(sample,line))*(sind(LatP))*(sind(C(sample,line)));
        
        Lon(sample,line) = LonP + atand(a1/(a2-a3)) ; % storing longitude value

%% For latitude calculation
        v1 = (cosd(C(sample,line)))*sind(LatP)  ; % temporary variables v1,v2,v3,v4
        v2 = y(sample,line)/P    ;
        v3 = (sind(C(sample,line)))*(cosd(LatP));
        v4 = v2*v3;
        
        Lat(sample,line) = asind(v1+v4); % storing latitude value
    end
end

max_Lon = (max(max(Lon)));
max_Lat = (max(max(Lat)));

min_Lat = (min(min(Lat)));
min_Lon = (min(min(Lon)));

%% 
for i = 1: samples_per_line
    for j = 1:lines
        if (reshaped_array(i,j) < -3.4028227e+31)
             final(i,j)=0.010000000000000; % assigning a fixed value so that imshow can work fine
                                           % without this step will just
                                           % see a black and white image
        end
    end
end

figure;
imshow(final',[]);  % displaying transpose of array
title('Vikram Landing site (November Mosaic made using LRO NAC images)')

toc % end timer