clear

set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
    'centimeters', 'PaperPosition', [0.2 0.2 20.5 29]); 
%define colormap
map = [    0.2510         0    0.2510
    0.2331         0    0.3045
    0.2151         0    0.3580
    0.1972         0    0.4115
    0.1793         0    0.4650
    0.1613         0    0.5185
    0.1434         0    0.5720
    0.1255         0    0.6255
    0.1076         0    0.6790
    0.0896         0    0.7325
    0.0717         0    0.7860
    0.0538         0    0.8395
    0.0359         0    0.8930
    0.0179         0    0.9465
         0         0    1.0000
    0.0091    0.0909    1.0000
    0.0182    0.1818    1.0000
    0.0273    0.2727    1.0000
    0.0364    0.3636    1.0000
    0.0455    0.4545    1.0000
    0.0545    0.5455    1.0000
    0.0636    0.6364    1.0000
    0.0727    0.7273    1.0000
    0.0818    0.8182    1.0000
    0.0909    0.9091    1.0000
    0.1000    1.0000    1.0000
    0.2500    1.0000    1.0000
    0.4000    1.0000    1.0000
    0.5500    1.0000    1.0000
    0.7000    1.0000    1.0000
    0.8500    1.0000    1.0000
    1.0000    1.0000    1.0000
    1.0000    1.0000    0.8333
    1.0000    1.0000    0.6667
    1.0000    1.0000    0.5000
    1.0000    1.0000    0.3333
    1.0000    1.0000    0.1667
    1.0000    1.0000         0
    1.0000    0.9333         0
    1.0000    0.8667         0
    1.0000    0.8000         0
    1.0000    0.7333         0
    1.0000    0.6667         0
    1.0000    0.6000         0
    1.0000    0.5333         0
    1.0000    0.4667         0
    1.0000    0.4000         0
    1.0000    0.3333         0
    1.0000    0.2667         0
    1.0000    0.2000         0
    1.0000    0.1333         0
    1.0000    0.0667         0
    1.0000         0         0
    1.0000         0    0.0909
    1.0000         0    0.1818
    1.0000         0    0.2727
    1.0000         0    0.3636
    1.0000         0    0.4545
    1.0000         0    0.5455
    1.0000         0    0.6364
    1.0000         0    0.7273
    1.0000         0    0.8182
    1.0000         0    0.9091
    1.0000         0    1.0000];
colormap(map);


fid = fopen('3D_Sq_UT_Br_400_20.txt');
%fid = fopen('test.txt');
a = fscanf(fid,'%e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e',[360 inf]);
Br_3D = a';
fclose(fid)

fid = fopen('1D_Sq_UT_Br_400_20.txt');
%fid = fopen('test.txt');
a = fscanf(fid,'%e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e',[360 inf]);
Br_1D = a';
fclose(fid)

phi = [0.5:1:359.5]';
theta = [.5:1:179.5]';

load coast
subplot(3,1,1)
imagesc(phi, theta, Br_3D, [-30, 30])
hold on;
plot([long; 360+long], 90-[lat; lat], '-k')
colorbar('h');

title 'UT: 20 h; 400 km; Br(3D), Br(1D), Br(3D)-Br(1D)'

subplot(3,1,2)
imagesc(phi, theta, Br_1D, [-30, 30])
hold on;
plot([long; 360+long], 90-[lat; lat], '-k')
colorbar('h');

subplot(3,1,3)
imagesc(phi, theta, Br_3D - Br_1D, [-10, 10])
hold on;
plot([long; 360+long], 90-[lat; lat], '-k')
colorbar('h');

print -depsc -tiff  'Sq_Br_400_UT_20.eps'

return

tt = load('fort.77')
plot(tt(:,1), tt(:,2), 'b')
hold on
plot(tt(:,1), tt(:,3), 'r')
hold on
plot(tt(:,1), tt(:,4), 'g')
title '3D; LT; Bp - b, Bt - r, Br - g; theta = 40; phi = 90; all coeff'

return 
