clear all
close all

load euv2001091

re=6371;
lons=0;
lats=deg2rad(90);
diss=8*re;

%Calculate the He+ total content along the ray paths
[thic,rpel,rpaz]=euvgb1(lons,lats,diss,heden,lonmat,latmat,disv,timmat,dno,year);

%For comparson conduct a simailar calculation for H+
[thyic,rpel,rpaz]=euvgb1(lons,lats,diss,hyden,lonmat,latmat,disv,timmat,dno,year);

%Rotate the image so Sun is downwards and interpolate to
%produce finer image.
[thp,rdp]=cart2pol(rpaz,rpel);               %convert to polar co-ordinates
thp2=thp-(timmat-6)*pi/12;                   %rotate points by required angle
[rpaz2,rpel2]=pol2cart(thp2,rdp);            %cartesian coorrdinates of rotated points

%Use interpolation to produce a finer profile and restrict to -45 to 45 deg.
[rpazf,rpelf]=meshgrid(linspace(-pi/4,pi/4,100));      %finer grid points
thicf=griddata(rpaz2,rpel2,log10(thic),rpazf,rpelf);   %interpolate finer grid points
thicf=10.^thicf;

%Use interpolation to produce a finer profile and restrict to -45 to 45 deg.
thyicf=griddata(rpaz2,rpel2,log10(thyic),rpazf,rpelf);   %interpolate finer grid points
thyicf=10.^thyicf;

%Ratio to plot results in Earth radii
eang=atan(re/diss);

figure(1)
pcolor(rpel2/eang,rpaz2/eang,log10(thic));
shading flat
axis equal
colorbar

figure(2)
pcolor(rpelf/eang,rpazf/eang,log10(thicf));
shading flat
axis equal
colorbar
xlabel('Earth radii','fontweight','bold','fontsize',12)
ylabel('Earth radii','fontweight','bold','fontsize',12)
set(gca,'fontweight','bold','fontsize',12)
tits=['GPID He^+ EUV   UT=',num2str(timmat2),'   ',int2str(day),'/',int2str(month),'/',int2str(year)];
title(tits,'fontweight','bold','fontsize',14)

figure(3)
pcolor(rpelf/eang,rpazf/eang,log10(thyicf));
shading flat
axis equal
colorbar
xlabel('Earth radii','fontweight','bold','fontsize',12)
ylabel('Earth radii','fontweight','bold','fontsize',12)
set(gca,'fontweight','bold','fontsize',12)
tits=['GPID H^+ EUV   UT=',num2str(timmat2),'   ',int2str(day),'/',int2str(month),'/',int2str(year)];
title(tits,'fontweight','bold','fontsize',14)

figure(4)
plot(rpaz2,rpel2,'x',rpazf,rpelf,'o')
