% file: vhx3d_proj_fluxes.sl
%
% Demonstration script file to convert the shl_fluxes
% to 3D projected color-coded image.

#iffalse
% - - - Summary of usage steps:

% [linux] hydra

% Setup the desired simulation, e.g., for SN96crE_out :
.source dia_96crE
;

% Timestep to do:
tid=65;
% Can change the energy bands if desired:
vhx_fluxehis = [8.0, 1.3, 8.0];
vhx_fluxelos = [0.5, 0.9, 5.0];

% this fills the shl_fluxes values
vhx_calc_shellfluxes(tid);

% then do all the stuff in here (below) by doing:
.source vhx3d_proj_fluxes
;

% Save an image to a file:
% encode the energy ranges, flux-exponent value,
% and the color assignments in the name (somehow/if-desired):
save_image_file="96crE_step65_BGR058009135080flexp25_project.png";
% ALL=blue LO=green HI=red 
v3d_project(arrEB2/sum(arrEB2),arrEB1/sum(arrEB1),0.5*arrEB0/sum(arrEB0));

% - - -
#endif


% Prototype routine to make 3D projected image from shell fluxes
%  vhx3d_proj_fluxes <-- script ?  or a routine(flexp) ?  or just as far as:
%    (arrflx0, arrflx1, arrflx2) = vhx3d_arrs_from_fluxes (flexp)

% Set the non-linear contrast reduction exponent:
flexp = 0.25;  %  <-- ***


% useful number, indices:
szshells=length(shl_rad);
shlinds = [1:szshells-1:1];

% max, min radii that have flux (in pc)
sel = where( (shl_fluxes[0,*] + shl_fluxes[1,*] + shl_fluxes[2,*]) > 0.0);
shlradmax = shl_rad[max(sel)]/Const_pc;
shlradmin = shl_rad[min(sel)-1]/Const_pc;

% Volume of each shell, use to convert to the flux per volume (4pir^2 dr, in pc^3):
shlvolpc3 = 0.0*shl_rad;
shlvolpc3[0]=(4.0*PI/3.0)*(shl_rad[0]/Const_pc)^3;
shlvolpc3[shlinds] = (shl_rad[shlinds] -
   shl_rad[shlinds-1])/Const_pc;
shlvolpc3[shlinds] *= 4.0*PI*(0.5*(shl_rad[shlinds] +
   shl_rad[shlinds-1])/Const_pc)^2;

% make arrays with double the number of points
% so that lookup will be histogram-ish
doubrad = 0.95*shl_rad[shlinds-1] + 0.05*shl_rad[shlinds];
doubrad = [doubrad, 0.05*shl_rad[shlinds-1] + 0.95*shl_rad[shlinds] ];
doubrad /= Const_pc;
% put into sorted order (i.e., shuffled or interleaved):
sord = array_sort(doubrad);
doubrad = doubrad[sord];
% make the doubled fluxes arrays and order/interleave them too:
doubEB0 = [shl_fluxes[0,shlinds]/shlvolpc3[shlinds],
 shl_fluxes[0,shlinds]/shlvolpc3[shlinds]];
doubEB0 = doubEB0[sord];
doubEB1 = [shl_fluxes[1,shlinds]/shlvolpc3[shlinds],
 shl_fluxes[1,shlinds]/shlvolpc3[shlinds]];
doubEB1 = doubEB1[sord];
doubEB2 = [shl_fluxes[2,shlinds]/shlvolpc3[shlinds],
 shl_fluxes[2,shlinds]/shlvolpc3[shlinds]];
doubEB2 = doubEB2[sord];

% clip any very low (or 0) values:
%%doubEB0[where(doubEB0 < 1.e-4*max(doubEB0))] = 1.e-4*max(doubEB0);
%%doubEB1[where(doubEB1 < 1.e-4*max(doubEB1))] = 1.e-4*max(doubEB1);
%%doubEB2[where(doubEB2 < 1.e-4*max(doubEB2))] = 1.e-4*max(doubEB2);

% plot the flux-per-volume vs radius lookup curves
xlin; xrange(0.97*shlradmin, 1.02*shlradmax);
sel = [doubEB0, doubEB1, doubEB2][where([doubEB0, doubEB1, doubEB2] > 0.0)];
yrange(min(sel), 1.2*max(sel)); ylog;
label("Radius (pc)","Flux-from-shell per volume (erg/s/cm^2 per pc^3)",
  "Based on shells par-file: "+vhx_outdir+"/model_"+string(1000+tid)+".par");
plot_bin_integral;  % incase doing a Histogram plot

plot(doubrad, doubEB0, 1);
oplot(doubrad, doubEB1, 2);
oplot(doubrad, doubEB2, 4);


% Make sure we have v3d for this...
#ifexists v3d_setup

% setup the v3d range:
v3d_setup(1.05*shlradmax,67);

% Make 3D volume arrays from the 3 energy bands
%   note non-linear compressin with flexp:
arrEB0 = v3d_sphere_rlup (shlradmin, shlradmax, doubrad, (doubEB0)^flexp);
arrEB1 = v3d_sphere_rlup (shlradmin, shlradmax, doubrad, (doubEB1)^flexp);
arrEB2 = v3d_sphere_rlup (shlradmin, shlradmax, doubrad, (doubEB2)^flexp);


% and look at the projection of these:
% make it big and save in a file:
v3d_image_npix=512;  % only for the on-screen display
% show to screen:
save_image_file=NULL;

% pick an image:
v3d_project(arrEB0);
v3d_project(arrEB1);
v3d_project(arrEB2);
% LO=red HI=green
v3d_project(arrEB1/sum(arrEB1),arrEB2/sum(arrEB2));
% LO=red HI=green ALL=blue
v3d_project(arrEB1/sum(arrEB1),arrEB2/sum(arrEB2),0.5*arrEB0/sum(arrEB0));
% LO=green HI=red
v3d_project(arrEB2/sum(arrEB2),arrEB1/sum(arrEB1));
% LO=green HI=red ALL=blue
v3d_project(arrEB2/sum(arrEB2),arrEB1/sum(arrEB1),0.5*arrEB0/sum(arrEB0));

% can put the next image in a file by pre-specifying:
%%  save_image_file="96crE_step65_project.png";

#endif