% file: e0102_model.sl        10/28/12 - dd
%
% Define a general model useful for E0102, consisting of a
% Ring plus a Spherical shell and a Blob too.
% See the parameter definitions below for more specifics.

% - - - - -
% Setup general info about the source model...

s3d_model_name="Ring-Sphere-Blob";
% The distance, in kpc.
s3d_distance = 60.0;
% Date of model (in decimal year)
s3d_date = 2000.0;   

% The v3d coordinate system size and resolution
s3d_radius = 30.0;  % arc seconds
s3d_resolution = 57; % ~ 1/2 arc sec resolution

% The foreground absorption
s3d_NH22 = 0.08;
% This value is used in a nominal absorption model.
% Set to NULL to ignore.
% The model can be user customized if desired, see
% s3d_update_nh() in source3d.sl for details.
s3d_SHOW_NHARF_PLOT = 0;
% Nominal NH-arf method and max-slope:
% (A slope of 3 is 55-75% efficient generally for E0102, NH22=0.08.
%  Changing the slope may improve efficieny depending on the instrument
%  and the NH, e.g., for NH22 > 1.0 a slope
%  of 4 is more efficient.)
s3d_NHarf_method=1, s3d_NHarf_slope=3.0;

% Nominal spectral binning parameters
s3d_lam_short=Const_hc/8.0, s3d_lam_long=Const_hc/0.3, s3d_lam_dlog=0.0003;


% - - - - -
% Define the Custom Parameters 
% The structure "s3d_ps" is hard-defined, but
% the user can define the tags as desired:
%
% For the generalized ring+sphere+blob model
% the parameters here are:
% - xc,yc : center of ring and sphere
%    Ring:
% - Rtype:
%     1-cylinder(Ra=axial 'radius'),
%     2-spherical-belt(Ra=half-opening angle)
% - Rthphi: ring axis orientation array, [theta, phi].
% - Ri,Ro, Ra: inner, outer radii; axial extent.
% - az0, azamp: azimuthal density variation parameters,
%     az0 is the  highest-density azimuth (-180 to 180 degrees)
%     azamp is fractional density amplitude, azamp=NULL for none.
%     (see also the specific equation in code below.)
%    Spherical-shell:
% - Si,So: inner outer radii of spherical shell
%    Blob (sphere):
% - Bxyz: location of blob, [x, y, z].
% - Br: radius of the blob.
% 
% - norm1,2,3 are included and directly set the 
%    s3d_comp[1,2,3].norm values.
% 
s3d_ps = struct {xc,yc, Rtype,Rthphi,Ri,Ro,Ra,az0,azamp,
   Si,So,Bxyz,Br,norm1,norm2,norm3};


% - - - - -
% Define and read in the spectra that can/will be used with the components.

% Generate an isis par file
xspec_xset ("NEIVERS", "1.1");
() = xspec_abund("angr");
% Create and use a .par file,
% a single vnpshock model (remember, the N_H is already included)
fit_fun("vnpshock(1)");
set_par("vnpshock(1).norm", 0.0169);
set_par("vnpshock(1).kT_a", 0.636);
set_par("vnpshock(1).kT_b", 0.211);
set_par("vnpshock(1).Tau_l", 3.27e9);
set_par("vnpshock(1).Tau_u", 4.59e11);
% abundances C--Ni to 0:
set_par([6:16],0.0);
% some non-zero abundances:
set_par("vnpshock(1).O", 0.7);
set_par("vnpshock(1).Ne", 1.3);
set_par("vnpshock(1).Mg", 0.7);
set_par("vnpshock(1).Si", 0.15);
set_par("vnpshock(1).Fe", 0.0746);
save_par("temp_e0102_vnpshock.par");

% load and evaluate the par file to set the spectrum:
s3d_load_spectrum(1,"temp_e0102_vnpshock.par");


% - - - - -
% Define my 3D geometric components
%
%   Note that the s3d_update_comp routine will be called
% as part of model updating and its chief purpose is
% to update the array values.
%   Component values that are set once-for-all or that
% may be changed by the user should be defined here outside
% of the 'update routine.
variable iarr;

% Component index 1 - one-time/manual parameters:
iarr=1;
s3d_comp[iarr].name = "Ring";
s3d_comp[iarr].type = "xray";
s3d_comp[iarr].ignore = 0;
s3d_comp[iarr].norm = 1.0;
s3d_comp[iarr].ispec= 1;
s3d_comp[iarr].rndint= 1.0;
s3d_comp[iarr].rndclr= [1., 0., 0.];
s3d_comp[iarr].vsxyz= [0.,0.,0.];
s3d_comp[iarr].vtype= -1;  % no velocity field
s3d_comp[iarr].velval=0.0;
s3d_comp[iarr].voxyz= NULL;   
s3d_comp[iarr].vrthphi= NULL;   

% Component index 2 - one-time/manual parameters:
iarr=2;
s3d_comp[iarr].name = "Sphere-shell";
s3d_comp[iarr].type = "xray";
s3d_comp[iarr].ignore = 0;
s3d_comp[iarr].norm = 1.0;
s3d_comp[iarr].ispec= 1;
s3d_comp[iarr].rndint= 3.0;
s3d_comp[iarr].rndclr= [0., 0., 1.];
s3d_comp[iarr].vsxyz= [0.,0.,0.];
s3d_comp[iarr].vtype= -1;  % no velocity field
s3d_comp[iarr].velval=0.0;
s3d_comp[iarr].voxyz= NULL;   
s3d_comp[iarr].vrthphi= NULL;   

% Component index 23- one-time/manual parameters:
iarr=3;
s3d_comp[iarr].name = "Blob";
s3d_comp[iarr].type = "xray";
s3d_comp[iarr].ignore = 0;
s3d_comp[iarr].norm = 1.0;
s3d_comp[iarr].ispec= 1;
s3d_comp[iarr].rndint= 0.01;
s3d_comp[iarr].rndclr= [0., 1., 0.];
s3d_comp[iarr].vsxyz= [0.,0.,0.];
s3d_comp[iarr].vtype= -1;  % no velocity field
s3d_comp[iarr].velval=0.0;
s3d_comp[iarr].voxyz= NULL;   
s3d_comp[iarr].vrthphi= NULL;   


% - - - - -
% Define the calculation/update of the component geometry arrays:
public define s3d_update_comp()
{
   % This is an example routine that the user will replace with their
   % own version, generally by modifying the lines between the - - - -'s   
   variable iarr, this_arr;
   % set the v3d parameters
   v3d_setup(s3d_radius, s3d_resolution);
   % Note that some values can be hard wired 
   % and others come from the (user defined) s3d_ps structure.


   % - - - - - - - - -
   % Set the component norms from the s3d_ps.normN values:
   % Save the comp[].norm values in to s3d_ps:
   s3d_comp[1].norm = s3d_ps.norm1;
   s3d_comp[2].norm = s3d_ps.norm2;
   s3d_comp[3].norm = s3d_ps.norm3;


   % - - - - - - - - -
   % Component index 1
   iarr=1;

   % Using an azimuthal variation?
   if (s3d_ps.azamp != NULL)
     {
	% calculate the azimuthal lookup values
	variable azlups=179.99*[-30:30:1]/30.0;
	variable intlups=(1.0+s3d_ps.azamp*cos((azlups-s3d_ps.az0)*PI/180.0))^2;
     }
   
   % which type of ring?
   if (s3d_ps.Rtype == 1)
     {	
	% the Ring is a simple cylinder, Ra is the axial half-length.
	% Include an azimuthal variation of intensity?
	if (s3d_ps.azamp == NULL)
	  {
	     % just a simple cylinder
	     this_arr = v3d_cylinder(s3d_ps.Ri, s3d_ps.Ro,
				     -1.0*s3d_ps.Ra,s3d_ps.Ra, 
				     [s3d_ps.xc, s3d_ps.yc, 0.0], s3d_ps.Rthphi);	     
	  }
	else
	  {
	     % cylinder with azimuthal lookup  
	     this_arr = v3d_cyl_azlup(s3d_ps.Ri, s3d_ps.Ro, -1.0*s3d_ps.Ra,s3d_ps.Ra, 
				      [s3d_ps.xc, s3d_ps.yc, 0.0], s3d_ps.Rthphi,
				      azlups, intlups);
	  }
     }
   
   % Only one other Ring option, so use "else" instead of
   %          "if (s3d_ps.Rtype == 2) then"
   else
     {
	% the Ring is a truncated spherical-shell, Ra is half-opening angle.
	% start with a spherical shell with opening angle	
	this_arr = v3d_sphere(s3d_ps.Ri, s3d_ps.Ro, [s3d_ps.xc, s3d_ps.yc, 0.0]);
	variable anticone = 1.0 - v3d_cone(0.0,90.0-s3d_ps.Ra,-1.0*s3d_ps.Ro,s3d_ps.Ro,
					   [s3d_ps.xc, s3d_ps.yc, 0.0], s3d_ps.Rthphi);
	anticone[where(anticone < 0.8)]=0.0;
	this_arr = this_arr * anticone;

	% Include an azimuthal variation of intensity?
	if (s3d_ps.azamp != NULL)
	  {
	     % multiply by a cylinder with the azimuthal variation
	     this_arr = this_arr * v3d_cyl_azlup(s3d_ps.Ri/10.0, 1.1*s3d_ps.Ro, -1.0*s3d_ps.Ro, s3d_ps.Ro, 
				      [s3d_ps.xc, s3d_ps.yc, 0.0], s3d_ps.Rthphi, azlups, intlups);
	  }
     }
   
   % Normalize the array:
   s3d_comp[iarr].arr = this_arr/sum(this_arr);

   % Ancillary component information that is is updated with parameter changes:
   %%s3d_comp[iarr].velval=2000.0/s3d_radius;

   % - - - - - - - - -      
   % Component index 2
   iarr=2;
   this_arr = v3d_sphere(s3d_ps.Si, s3d_ps.So, [s3d_ps.xc, s3d_ps.yc, 0.0]);
   % Normalize the array:
   s3d_comp[iarr].arr = this_arr/sum(this_arr);

   % Ancillary component information that is is updated with parameter changes:
   %%s3d_comp[iarr].velval=2000.0/s3d_radius;

   % - - - - - - - - -      
   % Component index 3
   iarr=3;
   this_arr = v3d_sphere(0.0, s3d_ps.Br, s3d_ps.Bxyz);
   % Normalize the array:
   s3d_comp[iarr].arr = this_arr/sum(this_arr);

   % Ancillary component information that is is updated with parameter changes:
   %%s3d_comp[iarr].velval=2000.0/s3d_radius;
   
}


% - - - - -
% Set the initial parameter values:

% center offsets
s3d_ps.xc = 0.0;
s3d_ps.yc = 3.0;

% Ring type and orientation:
% 
%   E0102 tilted cylinder:
%%s3d_ps.Rtype=1;
%%s3d_ps.Rthphi=[270.0-18.0, 5.3];
% cylinder half-length:
%%s3d_ps.Ra = 8.0;
%
%   in-plane Hughes'94-like ring:
s3d_ps.Rtype=2;
s3d_ps.Rthphi=[270.001, 0.001];
% cylinder half-opening angle in degrees:
s3d_ps.Ra = 56.0;

% ring inner/outer radii:
s3d_ps.Ri = 12.0;
s3d_ps.Ro = 17.0;

% ring azimuthal variation:
s3d_ps.az0 = -135.0;
s3d_ps.azamp = 0.10;

% Spherical shell inner, outer radii:
s3d_ps.Si = 17.0;  % have it start at Ro
s3d_ps.So = 22.0;

% Blob location and radius:
s3d_ps.Bxyz = [0.0,0.0,10.0];
s3d_ps.Br = 1.5;

% set the norm(s)
s3d_ps.norm1 = 0.62;
s3d_ps.norm2 = 0.52;
s3d_ps.norm3 = 0.020;


% - - - - -
% Evaluate the model
s3d_update;

% Print the model info
s3d_list_model;


% Can now look at the model:
%
% a 3-D rotatable view:
%%s3d_view_comps;
% "sliced" to show interior:
%%s3d_view_comps(1);
%
% a projected-on-the-sky view:
s3d_proj_comps;

% Put the image above into a file by doing this before the
% 'proj_comps line above:
% save_image_file="e0102_model.png";

% Can change parameter(s) and update/re-view, e.g.,
%   hydra> s3d_ps.Ri=5.0;
%   hydra> s3d_ps.Ro=25.0;
%   hydra> s3d_update;
%   hydra> s3d_view_comps(1);
%   hydra> s3d_proj_comps;
%   
% etc.
%