%
% file: v3dex_mc_vels.sl
% 
% Examples of V3D routines 
%
% - - - - - Monte Carlo and Velocity Effects - - - - - 
%

% Set the scale and resolution (default):
v3d_setup(100.0,37);


% Make a thin disk (cylinder) oriented in space:
variable daxis = [290.0,60.0];
variable mydisk = v3d_cylinder(10.0,90.0, -5.0,5.0, , daxis);

% look at it in 3D and in projection:
#ifexists volview
v3d_view(mydisk);
#endif
v3d_project(mydisk);


% Generate MC points from the disk
variable xs,ys,zs;
(xs,ys,zs) = v3d_mc_points(mydisk, 15000L);

% Look at these points on the sky:
v3d_evtimg(xs,ys);


% Assign line-of-sight velocity energy factors to the points
% and look at the events with color coding velocity
% for different velocity cases
%
variable enfs;

% Expansion at constant velocity along r_hat:
%  v_const = 1000 km/s
enfs = v3d_doppler_points(xs,ys,zs, , 1000.0, 0);
%
% print the max/min los velocities of these:
message(" LOS velocity range: "+string((min(enfs)-1.0)*3.e5)+" to " + string((max(enfs)-1.0)*3.e5) );
% And view these "events" color coded by los velocity: 
v3d_evtimg(xs,ys,enfs, 1.0+[-3000.0,-500.0]/3.e5, 1.0+[-500.0,500.0]/3.e5, 1.0+[500.0,3000.0]/3.e5);


% Expansion proportional to r  w/ velocity along r_hat:
%   |v| = 1000 km/s * (r/70.0)
enfs = v3d_doppler_points(xs,ys,zs, , 1000.0/70.0, 1);
%
% print the max/min los velocities of these:
message(" LOS velocity range: "+string((min(enfs)-1.0)*3.e5)+" to " + string((max(enfs)-1.0)*3.e5) );
% And view these "events" color coded by los velocity: 
v3d_evtimg(xs,ys,enfs, 1.0+[-3000.0,-500.0]/3.e5, 1.0+[-500.0,500.0]/3.e5, 1.0+[500.0,3000.0]/3.e5);


% Rigid-body rotation about the disk axis
%   |v| = 1000 km/s * (r/70.0)
enfs = v3d_doppler_points(xs,ys,zs, , 1000.0/70.0, 2, , daxis);
%
% print the max/min los velocities of these:
message(" LOS velocity range: "+string((min(enfs)-1.0)*3.e5)+" to " + string((max(enfs)-1.0)*3.e5) );
% And view these "events" color coded by los velocity: 
v3d_evtimg(xs,ys,enfs, 1.0+[-3000.0,-500.0]/3.e5, 1.0+[-500.0,500.0]/3.e5, 1.0+[500.0,3000.0]/3.e5);


% Keplerian rotation about the disk axis
%   |v| = 500 km/s / sqrt(r/100.0)
enfs = v3d_doppler_points(xs,ys,zs, , 500.0*sqrt(100.0), 3, , daxis);
%
% print the max/min los velocities of these:
message(" LOS velocity range: "+string((min(enfs)-1.0)*3.e5)+" to " + string((max(enfs)-1.0)*3.e5) );
% And view these "events" color coded by los velocity: 
v3d_evtimg(xs,ys,enfs, 1.0+[-3000.0,-500.0]/3.e5, 1.0+[-500.0,500.0]/3.e5, 1.0+[500.0,3000.0]/3.e5);


% If we have acess to usual 1D plotting functions, e.g. via ISIS, then
% We can bin and plot the los values:
#ifexists hplot
variable bins =  [0.990:1.010:0.0001];
variable dbins, blos, bhis;
dbins = bins[1] - bins[0];
blos = bins - 0.5*dbins;
bhis = bins + 0.5*dbins;
variable histenfs = hist1d(enfs, bins);
hplot( blos, bhis, histenfs);
#endif