 %Webb3DModelInputLT
%11/13/07 by David Galvan

%Import Phillip Webb's model density information, to be used as "a priori"
%data for 3D tomography.  This version outputs reformat5 in Local Time
%instead of longitude, which allows us to imagine the prior density model
%is fixed in local time over a longer period of universal time.  Basically
%allows us to take a longer swath of TEC data.

%VI = interp3(X,Y,Z,V,XI,YI,ZI);

%fid = fopen('/Users/dgalvan/Documents/Research/3DTomography/webb_model_reformatted_test.txt'

%[Elev_Angle(:,1), STEC(:,1), SecJ2000(:,1), Xleo(:,1), Yleo(:,1), Zleo(:, 1), Xgps(:,1), Ygps(:,1), Zgps(:, 1)] = textread('/Users/dgalvan/Documents/Research/SpringAGU2007/tec9files/jaso04285.tec.9.txt', '%*f %*f %*f %*f %*f %*f %f %*f %f %*f %*f %*f %*f %*f %*f %f %*f %*f %*f %*f %f %f %f %f %f %f %*f %*f %*f %*f %*f %*f');

%It's a bit of a pain to automatically import the data from Webb's output,
%since it is in blocks, not columns.  So just import the file:
%/Users/dgalvan/Documents/Research/3DTomography/webb_model2.txt (for the 2004/346 event on the
%mac) or
%/Users/dgalvan/Documents/Research/3DTomography/GPID200215204_eMOD.txt (for
%the 2002/152 event on the mac)
%OR
%C:\My Documents\Research\3DTomography\webb_model2.txt (on the PC)

%using "File
%--> Import Data".  Once it is in the workspace:

%We're going to produce a new matrix called "reformat1" in which each row
%has a lon, a lat, and all the radial density values for that lat/lon.
%webb_model2 = data;
%webb_model2 = GPID200434601_eMOD;
webb_model2 = GPID200215204_eMOD;

reformat1 = []; %initialize this new matrix

for i = 1:(size(webb_model2, 1)./15); %For each block of data in webb_model2 (each single lat/lon point with a bunch of radial densities)
    %get the lon value from the 1st row of that block, first column
    lon = webb_model2(((i-1)*15 + 1),1);
        
    %get the lat value from the 2nd row of that block, first column
    lat = webb_model2(((i-1)*15 + 2),1);
    
    full_density_row = [];
    for j = 3:15; %now, within a given block of 15 lines, j represents a specific single line.  We're collecting the density values from each block, in order to put them in the same row.
        subrow = webb_model2(((i-1)*15 + j), 1:8); %collect one of the 13 rows of density values
        full_density_row = [full_density_row, subrow]; %and append it to build the "full_density_row" of all density values for that lat/lon.
    end
    
    reformat1 = [reformat1; lon, lat, full_density_row]; %take the lat, lon, and all density values for that block, and append it to build the reformat1 matrix.
end

%Now, since the last line of density values had only 4 columns of useful
%data, instead of the 8 columns the previous lines had, shave off the last
%4 columns of reformat1.

reformat2 = reformat1(:,1:102);

%and change units from electrons / m^3 to electrons / cm^3 by dividing by
%10^6

reformat2(:,3:102) = (10^(-6)) * reformat2(:,3:102);

%Now, one more reformat, such that we get a matrix with the following
%columns in each row: lon, lat, r, and density.

%first get the r_values for the 100 density columns 

R_earth = 6371; %radius of the Earth in km

r_shells_webb = R_earth + [500   547   595   645   697   750   805   862   921   981 1044  1109  1176  1245  1316  1390  1466  1545  1626  1710 1797  1887  1980  2075  2174  2277  2382  2491  2604  2720 2841  2965  3093  3226  3363  3504  3650  3801  3958  4119 4285  4457  4635  4819  5008  5204  5407  5616  5832  6056 6286  6524  6771  7025  7288  7559  7840  8129  8429  8738 9057  9387  9728 10080 10444 10820 11208 11610 12024 12452 12895 13352 13824 14312 14815 15336 15874 16429 17003 17596 18209 18842 19496 20171 20869 21590 22335 23104 23899 24721 25569 26446 27351 28287 29253 30251 31283 32349 33449 36000]';



reformat3=[];

%Warning!  The following for-loop takes ~30-40 minutes to run!
for k = 1:size(reformat2,1); %for each lon/lat value in reformat2
    
    for m = 3:size(reformat2,2);
        %create the row of reformat3, which has: lon, lat, r(kim), and
        %density (electrons / cm^2)
        reformat3_row = [reformat2(k,1), reformat2(k,2), r_shells_webb(m-2), reformat2(k,m)];
        reformat3 = [reformat3; reformat3_row];
    end
end

%Now, Tomography3D.m, which will be using the output of this code, uses
%longitudes from 0 to 360 degrees, instead of -180 to 180.  You can change
%this later if you want, but you'll have to spend some careful time with
%Tomography3D.m.  For now, the easiest way to ensure we get a prior without
%NaN's is to convert all negative longitude values to positive values.  IE:
%-90 deg longitude --> 270 longitude.  So just add the negative value to
%360.
reformat3_poslon = [];
for q = 1:size(reformat3,1);
    if reformat3(q,1) < 0;
        reformat3_poslon(q,1) = 360 + reformat3(q,1);
        reformat3_poslon(q,2:4) = reformat3(q,2:4);
    elseif reformat3(q,1) >= 0;
        reformat3_poslon(q,1:4) = reformat3(q,1:4);
    end
    
end
        
%Finally, convert the lon,lat,r values to XYZ:

[cart_points(:,1), cart_points(:,2), cart_points(:,3)] = sph2cart((pi/180).*reformat3(:,1) , (pi/180).*reformat3(:,2), reformat3(:,3));


reformat4 = [cart_points , reformat3(:,4)];

% %make the webb intersection points
% [THETA,PHI,R] = meshgrid(0:10:360, -90:5:90, r_shells_webb);

%LT = longitude * (24 hrs / 360 degrees) + UT

UT =0;  %The Universal time at which Webb's model is set.
%NOTE: because this particular webb model happens to be set at UT = 0, we
%can simply say that LT_lon = geographic lon.  We can't do this if the webb
%model is not set at 00 UT.  In any other case, we'd have to convert using:

%LT = lon * (24/360) + UT
%LT_lon = LT * (360/24)
% or, simplified:
%LT_lon = lon + UT*(360/24)
%
%or lon = LT_lon - UT*(360/24)

%which will be used below when finding the density values stored in
%reformat3's geographic longitude bin

%make the webb intersection points using local time longitude, LT_lon.  We still need to
%remember that webb's model outputs every 10 degrees of longitude, which
%corresponds to 10*(24/360) = 0.667 hours of local time.
%[LT_lon,PHI,R] = meshgrid(0:10:360, -90:5:90, r_shells_webb);
lon_res = 10;
lat_res = 5;
[lon,PHI,R] = meshgrid(0:lon_res:360, -90:lat_res:90, r_shells_webb);
%LT_lon = lon;

LT_lon = lon + UT*(360/24) ;

% 
% %Doing the above addition to create LT_lon can push LT_lon values past 360.
% %Also, since we'll be using interp3 in Tomography3DLT, I'd like to keep the
% %longitudes monotonic.  I don't really know if interp3 requires this but I
% %have a feeling it would cause problems, so I'm just going to do it this
% %way.
% 
% if LT_lon_temp(1,end,1) > 360;
%     LT_lon =  LT_lon_temp - ones*lon_res;
% elseif LT_lon_temp(1,end,1) <= 360;
%     LT_lon = LT_lon_temp;


%Now you must create reformat5, which is a 3D array that has the densities
%only at 1st dimension: LT_lon, 2nd: PHI/lat, 3rd: r_shells

%first, create an array of the correct dimensions, but filled with ones.
reformat5 = ones(size(LT_lon, 1), size(LT_lon, 2), size(LT_lon, 3));

for l = 1:length(LT_lon(1,:,1));  %for each LT_lon value
    
    for m = 1:length(PHI(:,1,1)); %for each lat value
       
        for n = 1:length(R(1,1,:)); %for each R value
            
            if (LT_lon(1,l,1) - UT*(360/24))==360;
                %find the index of reformat 3 that has those LT_lon, lat, R
                %values. . .  If it happens that LT_lon is 360 degrees, use
                %0 degrees longitude instead, since that's how it's stored
                %in reformat3_poslon.
                index = find(reformat3_poslon(:,1) == 0 & reformat3_poslon(:,2) == PHI(m,1,1) & reformat3_poslon(:,3) == R(1,1,n));
            else
                %find the index of reformat 3 that has those LT_lon, lat, R
                %values. . .
                index = find(reformat3_poslon(:,1) == (LT_lon(1,l,1) - UT*(360/24)) & reformat3_poslon(:,2) == PHI(m,1,1) & reformat3_poslon(:,3) == R(1,1,n));
            end
                %. . . and put that density value into the appropriate row (m,
            %lat, phi), column (l, THETA, lon), and slice (n,R)
            
            reformat5(m, l, n) = reformat3_poslon(index, 4);
        end
    end
end

% %To make it a closed loop, we must include a column for longitude 360
% %degrees at the end, which is identical to the 1st column (longitude 0
% %degrees).  This is just to make sure the interpolation covers the entire
% %range of longitudes.
% 
% reformat5 = [reformat5, reformat5(:,1,:)];

% LT_lon = [LT_lon, LT_lon(:,1,:)];
% PHI = [PHI, PHI(:,1,:)];
% R = [R, R(:,1,:)];

% save the workspace

save webbDONE2002152.mat;

%set slice positions for lon, lat, and r.
LT_lonslice = [180,0]; PHIslice = [0, 30]; Rslice = [8000];

figure;   
slice(LT_lon,PHI,R,reformat5,LT_lonslice, PHIslice,Rslice);
colorbar;
shading flat;

caxis([1 8000]);  %set minimum and maximum density values for colormap
view([0,0]);  %set the azimuth and elevation for the viewpoint
       


        
    