% file: v3d_shapes_demo.sl
%
% 12/19/06, 10/8/08 - dd
%

v3d_setup(100.0, 37);
variable shape, shape_name;

% The v3d_view images are 256x256 - save them as png.s
% The v3d_project images are nominally 2*37+1 on a side:
%   manually adjust them to 256x256 and then save as png.


% Sphere
shape = v3d_sphere(20.0, 40.0, [20.,40.,20.]) + 
  v3d_sphere(0.0, 20.0, [0.,0.,-50.0]);
shape_name = "v3d_sphere";
%
v3d_view(shape);
v3d_project(shape);


% Spheroid
shape = v3d_spheroid(40.0, 50.0, 1.8, ,[-20.0,70.]);
shape_name = "v3d_spheroid";
%
v3d_view(shape);
v3d_project(shape);


% Cylinder (disk)
shape = v3d_cylinder(10.0,20.0, 0.0,80.0, , [290.0,0.0]) + 
v3d_cylinder(40.0, 90.0, -10.0, 10.0, ,[-55.0, -30.0]);
shape_name = "v3d_cylinder";
%
v3d_view(shape);
v3d_project(shape);


% Cone  (jet)
shape = v3d_cone(20.0, 30.0, 30.0, 80.0, , [-40.,-10.]) + 
2.0* v3d_cone(0.0, 3.0, -200.0,200.0, , [50.0,60.0] );
shape_name = "v3d_cone";
%
v3d_view(shape);
v3d_project(shape);


% Torus
shape = ( v3d_torus(30.0,70.0, [20.0,0.,0.],[-40.,10.0]) -
 v3d_torus(45.0,55.0, [20.,0.,0.], [-40.,10.0]) )  + 
  v3d_torus(30.0,50.0, [-10.,0.,40.], [140.0,80.0]);
shape_name = "v3d_torus";
%
v3d_view(shape);
v3d_project(shape);


% Cube-slice
shape = v3d_cube(100.0,[-70.,0.,-100.]) * v3d_sphere(50.0,80.0);
shape_name = "v3d_cube-slice";
%
v3d_view(shape);
v3d_project(shape);


% Roche lobe
shape = v3d_roche(90., 0.333, 1.0, [-75.0,0.,.0],[0.0,0.0]) + 
  ( v3d_roche(90., 3.0, 1.0, [15.0,0.,0.0],[180.0,0.0]) - 
    v3d_sphere(0.0, 20.0, [15.0,0.,0.0])  );
shape_name = "roche";
%
v3d_view(shape);
v3d_project(shape);


% Ring of spheres
shape = shape = v3d_sphere_ring(70.0, 5.0, 17, , [260.0, -44.0]);


% Sphere with radial lookup tables
shape = v3d_sphere_rlup(1.0,90.0, [0.0,3.0,10.0,30.0,100.0], [10.0,5.0,1.0,3.0,0.5]);


% Cylinder with Azimuthal limits/lookup, v3d_cyl_azlup
azlups = 9.0*[-10:20:1];
vlups = 0.3 + cos(5.0*azlups*PI/180.0)^2;
shape = v3d_cyl_azlup (70.0, 90.0, -40.0, 40.0, , [-70.0,30.0], azlups, vlups);


% Make a 2d array and turn it to 3D...
% simple test...
%  array is i_lattitude, i_radius
variable arr2d = Float_Type [91,100];
% from 0 to 10 degrees fill it to 1/5 radius:
arr2d[[0:9],[0:20]] = 1.0;
% from 10 to 45 degrees fill it to 1/2 radius:
arr2d[[10:45],[0:49]] = 1.0;
% from 46 to 90 degrees fill it to full radius:
arr2d[[46:90],*] = 1.0;
test3d = v3d_2dto3d(arr2d, 20.0, 90.0, , [-20.0, 30.0]);
%
v3d_view(test3d);


% 2D hydro simulation to 3D volume
variable arr2d = h5_read("wrbub_1026.h5");
arr2d = arr2d  / 1.660e-24 ;   % Density in amu/cm^3
% flip the lattitude values
variable arrRT = 0.0*arr2d;
_for irad (0,199,1) arrRT[*,irad] = reverse(arr2d[*,irad]);
%
shape = v3d_2dto3d(arrRT, 22.0, 90.0, , [90.0,80.0]);
%
v3d_view(shape*v3d_cube(100.0,[0.,0.,-100.]) - 6.0);
v3d_project(shape);


% look at projected intensity
v3d_project(rt3d);
% look at 3D volume - subtract a low level
v3d_view(rt3d - 5.0);  
% same thing but zoom in
v3d_setup(90.0, 37);
rt3d = v3d_2dto3d(arrRT, 40.0, 200.0, [-90.,-90.,-90.], [90.0,80.0]);
% look at projected intensity
v3d_project(rt3d);
% look at 3D volume - subtract a low level
v3d_view(rt3d - 5.0);  


% Function of 3-d radius
% A sum of two Gaussians ~ e-r^2  with 10:1 flux ratio
shape = 10.0*(1./30.0^3)*exp(-1.0*v3d_r3dsq([0.,0.,0.])/(30.0^2));
shape = shape + (1./10.^3)*exp(-1.0*v3d_r3dsq([60.,50.,30.])/(10.0^2));
%
% look at it with v3d_view and v3d_project
v3d_view(shape);
v3d_project(shape);


% Generic function of variables, v3d_f_generic
shape = v3d_f_generic("custom function", , [-25.0,15.0]);


% Demonstrate histxy and reproject

% Read in some events using Event2D routines (could use fits readcol instead)
% E0102:
% X, Y events - E0102 from obsid 3828, centered on the central bright spot
% Load some events
e2di=15;
e2d_load_data("../Data/obs_3828/evt2", "X", "Y", [4094.1,4166.1], [0,0], [0.5,5.0]);
e2d_view_data (15, [ -30.0 , 30.0 , 1.5 ] , [ -30.0, 30.0 , 1.5 ] );

% The events are in arc seconds, set the v3d scale for this demo appropriately:
v3d_setup(25.0, 27);

% v3d_histxy:  bin x, y values into the central (z=0) x-y plane of the cube:
% 
e0102_flat = v3d_histxy (e2d_data[15].dxs, e2d_data[15].dys);

% set low-level values to 0:
sel = where(e0102_flat < 0.2*max(e0102_flat));
e0102_flat[sel]=0.0;

% save these to png files:
%   v3d_histxy_project.png
v3d_project(e0102_flat);
%   v3d_histxy_view.png , Y slider at -46 deg.:
v3d_view(e0102_flat);

% v3d_reproject: 
%
e0102_sphere = v3d_reproject(e0102_flat, v3d_sphere(20.0,25.0));

%   v3d_reproject_project.png
v3d_project(e0102_sphere);
%   v3d_histxy_view.png , Y slider at -46 deg.:
v3d_view(e0102_sphere);