%% euvgb1.m calculates the total helium ion content (THIC) along multiple
%% raypaths corresponding to that observed by the EUV Imager on IMAGE. 
%% The He+ input arrays should be generated with plasigb.m
%% 

function [thic,rpelmat,rpazmat]=euvgb1(lons,lats,diss,hedeno,lonmato,latmato,disv,timmat,dno,year)

% lons = geographic longitude location of IMAGE (rads) (scalar)
% lats = geographic latitude location of IMAGE (rads) (scalar)
% diss = radial distance of IMAGE (km) (scalar)
% hedeno = He+ density array (m^-3)
% lonmato = geographic longitude of He+ from plasigb.m (rads)
% latmato = geographic latitude of He+ from plasigb.m (rads)
% disv = radial distance of He+ from plasigb.m (km)
% timmat = UT (fractional hours) (scalar)
% dno = day of interest (scalar)
% year = year of interest (scalar)
% 
% thic = total helium ion content 
% rpelmat = ray path elevations thic calaculated along
% rpazmat = ray path azimuth thic calaculated along

%Generate required 3-D arrays
[lonmat,latmat,dismat]=meshgrid([lonmato(1,:) pi],latmato(:,1),disv);

%For interpolation need to make copy of -pi longitude He+ values and place at pi
sz=size(lonmat);
heden=zeros(sz);
heden(:,1:(sz(2)-1),:)=hedeno;
heden(:,sz(2),:)=hedeno(:,1,:);

%Specify the angular seperation and length of raypaths to be used.
re=6371;            %Earth radii (km)
angsep=1;           %angular seperation (deg)
rpdis=diss+8*re;    %length of ray paths (km)
nop=100;            %number of points to break each raypath into

%Determine the co-ordinates of raypaths
[rplon,rplat,rpdis,rpel,rpaz]=euvrayp(lons,lats,diss,angsep,rpdis,nop);

%Reset He+ densities of zero (points less 350km latitude) to 1 for interpolation.
heden(heden==0)=1;

%Use 3-D linear interpolation to determine the He+ density at each point on each ray path.
%So as not to overload the RAM interpolate 1000 rayspaths at a time
nopfl=size(rplon,1);                         %number of field lines
nopl=1000;                                   %number of points per loop
nol=ceil(nopfl/nopl);                        %number of loops required
herp=zeros(nopfl,nop);                       %He+ results matrix
for i=1:nol                                  %loop through each group of raypaths
 li=(i-1)*nopl+1;                            %lower index of i'th loop
 if i~=nol
  ui=i*nopl;                                 %lower index of i'th loop
 else
  ui=nopfl;
 end
 herpi=interp3(lonmat,latmat,dismat,log10(heden),rplon(li:ui,:),rplat(li:ui,:),rpdis(li:ui,:));
 herp(li:ui,:)=10.^herpi;
end

%Determine the required elevations and azimuth arrays
rpelmat=reshape(rpel,sqrt(size(rplon,1)),sqrt(size(rplon,1)));
rpazmat=reshape(rpaz,sqrt(size(rplon,1)),sqrt(size(rplon,1)));

%Determine which raypath points are "inside" the Earth, which is taken 
%to be any points with radial distance less than re+350 (since GPID doesn't
%determine He+ densities below this height). Reset these points to 0.
poie=rpdis<(re+350);          %points inside the Earth
herp(poie)=0;                 %reset poie points to 0

%If any points along a given raypath where inside the Earth, then reset
%all raypoint points on the farside of the Earth to 0 as well.
poie2=sum(poie,2)>0;                         %raypaths that passed through the Earth  
ind2=find(poie2);                            %index of poie2
if ~isempty(poie2)                           %check if any raypaths passed through the Earth
 for i=1:length(ind2)                        %loop through each such ray path
  ind2x=min(find(poie(ind2(i),:)==1));       %index of first point on raypath inside the Earth
  herp(ind2(i),ind2x:nop)=0;                 %reset first point and all subsquent points along ray path to zero
 end
end

%Determine if any points along raypaths are in the Earth-Sun shadow. Reset this
%points to zero.
[day,month]=daynum2(dno,year);                    %day and month number
timl=timut2l(timmat,rplon);                       %local time of each point on raypaths
pois=eshadow(day,month,year,timl,rplat,rpdis);    %determine if points in Earths shadow
herp(pois==1)=0;                                  %reset points in Earths shadow to 0

%Reset any NaN to the polar He+ density. These will be points that were outside the interpolation range
herp0=herp;
poinan=isnan(herp);                                    %determine if there are any NaN 
if sum(sum(poinan))>0 
 lonpol=rplon(poinan);                                 %polar geographic longitudes
 latpol=rplat(poinan);                                 %polar geographic latitudes
 dispol=rpdis(poinan);                                 %polar radial distances 
 [lonm,latm]=geo2mag3(lonpol,latpol);                  %geomagnetic co-ordinates of NaN points
 svp=ones(size(lonm));
 lsh=lshpar(-svp,latm,dispol);                         %L-shell of NaN points
 dflm=distfl2(latm,lsh,(re+350)*svp);                  %distance along field lines of NaN points
 kp=solarval(timmat,dno,year);                         %Kp index
 [hypol,hepol]=polarmod(lsh,latm,dflm,kp);             %polar densities of NaN points
 herp(poinan)=hepol;                                   %reset NaN points to polar He+ densities
end

% Calculate the actual distance between points along the ray paths.
[xc,yc,zc]=sph2cart(rplon,rplat,rpdis);                %cartesian co-ords
dsdf=sqrt((xc(:,2:nop)-xc(:,1:(nop-1))).^2+(yc(:,2:nop)-yc(:,1:(nop-1))).^2+(zc(:,2:nop)-zc(:,1:(nop-1))).^2);                                              %distance between points along ray paths

%Determine the THIC for each raypath
hemn=0.5*(herp(:,2:nop)+herp(:,1:(nop-1)));       %mean He+ density matrix
thicm=[zeros(size(rplon,1),1) hemn.*1e3.*dsdf];   %THIC difference matrix
thicmcs=cumsum(thicm,2);                          %cumulative sum of THIC matrix
thic=thicmcs(:,nop);                              %total THIC for each ray path

%Reshape THIC to required output size
thic=reshape(thic,sqrt(size(rplon,1)),sqrt(size(rplon,1)));
