Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-04-09 07:49:37

0001 #include <limits>
0002 
0003 #include "scuda.h"
0004 #include "squad.h"
0005 #include "sqat4.h"
0006 #include "stran.h"
0007 
0008 #ifdef WITH_OLD_FRAME
0009 #include "sframe.h"
0010 #else
0011 #include "sfr.h"
0012 #endif
0013 
0014 #include "sc4u.h"
0015 #include "ssincos.h"
0016 #include "ssys.h"
0017 
0018 #include "SLOG.hh"
0019 #include "SRng.hh"
0020 #include "SGenstep.h"
0021 #include "OpticksGenstep.h"
0022 #include "SFrameGenstep.hh"
0023 #include "NP.hh"
0024 
0025 
0026 const plog::Severity SFrameGenstep::LEVEL = SLOG::EnvLevel("SFrameGenstep", "DEBUG" );
0027 
0028 /**
0029 SFrameGenstep::CE_OFFSET
0030 -------------------------
0031 
0032 Interprets the CE_OFFSET envvar converting the values from the string (default "0,0,0")
0033 into the ce_offset vector with one or more float3.
0034 
0035 Typically CE_OFFSET "0.,0.,0." corresponding to local frame origin.
0036 The string "CE" or "ce" in the envvar is special cased, causing
0037 the offset to be set to the ce provided by the argument.
0038 
0039 **/
0040 
0041 void SFrameGenstep::CE_OFFSET(std::vector<float3>& ce_offset, const float4& ce ) // static
0042 {
0043     const char* ekey = "CE_OFFSET" ;
0044     const char* val = ssys::getenvvar(ekey);
0045 
0046     bool is_CE_string = val == nullptr ? false : ( strcmp(val, "CE")== 0 || strcmp(val, "ce")== 0 )  ;
0047     float3 offset = make_float3(0.f, 0.f, 0.f );
0048 
0049 
0050     if(is_CE_string)   // this is not typically used anymore
0051     {
0052         offset.x = ce.x ;
0053         offset.y = ce.y ;
0054         offset.z = ce.z ;
0055         ce_offset.push_back(offset);
0056     }
0057     else
0058     {
0059         std::vector<float>* fvec = ssys::getenvfloatvec(ekey, "0,0,0");
0060         unsigned num_values = fvec->size() ;
0061         bool multiple_of_3 = num_values % 3 == 0 ;
0062         unsigned num_offset = num_values/3 ;
0063         assert(fvec);
0064 
0065         LOG(LEVEL)
0066            << " ekey " << ekey
0067            << " num_values " << num_values
0068            << " multiple_of_3 " << multiple_of_3
0069            << " num_offset " << num_offset
0070            ;
0071 
0072         LOG_IF(fatal, !multiple_of_3) << " not multiple_of_3 num_values " << num_values << " ekey " << ekey  ;
0073         assert(multiple_of_3);
0074 
0075         for(unsigned i=0 ; i < num_offset ; i++)
0076         {
0077             offset.x = (*fvec)[i*3+0] ;
0078             offset.y = (*fvec)[i*3+1] ;
0079             offset.z = (*fvec)[i*3+2] ;
0080             ce_offset.push_back(offset);
0081         }
0082     }
0083 
0084     bool with_offset = ce_offset.size() > 0 ;
0085     LOG(LEVEL)
0086          << "ekey " << ekey
0087          << " val " << val
0088          << " is_CE_string " << is_CE_string
0089          << " ce_offset.size " << ce_offset.size()
0090          << " ce " << ce
0091          << " with_offset " << ( with_offset ? "YES" : "NO " )
0092          ;
0093 
0094     LOG(LEVEL) << Desc(ce_offset) ;
0095 
0096     LOG_IF(fatal, !with_offset) << " at least one ce_offset is required " ;
0097     assert(with_offset);
0098 
0099 
0100 }
0101 
0102 std::string SFrameGenstep::Desc(const std::vector<float3>& ce_offset )
0103 {
0104     std::stringstream ss ;
0105     ss << "SFrameGenstep::Desc ce_offset.size " << ce_offset.size() << std::endl ;
0106     for(unsigned i=0 ; i < ce_offset.size() ; i++)
0107     {
0108         const float3& offset = ce_offset[i] ;
0109         ss << std::setw(4) << i << " : " << offset << std::endl ;
0110     }
0111 
0112     std::string s = ss.str();
0113     return s ;
0114 }
0115 
0116 std::string SFrameGenstep::Desc(const std::vector<int>& cegs )
0117 {
0118     std::stringstream ss ;
0119     ss << " size " << cegs.size() << "[" ;
0120     for(unsigned i=0 ; i < cegs.size() ; i++) ss << cegs[i] << " " ;
0121     ss << "]" ;
0122     std::string s = ss.str();
0123     return s ;
0124 }
0125 
0126 
0127 
0128 void SFrameGenstep::GetGridConfig_(std::vector<int>& cegs,  const char* ekey, char delim, const char* fallback )
0129 {
0130     ssys::getenvintvec(ekey, cegs, delim, fallback );
0131     LOG(LEVEL)
0132         << " ekey " << ( ekey ? ekey : "-" )
0133         << " fallback " << ( fallback ? fallback : "-" )
0134         << " delim " << delim
0135         << " cegs.size " << cegs.size()
0136         ;
0137 
0138 
0139     StandardizeCEGS(cegs);
0140     assert( cegs.size() == 0 || cegs.size() == 8 );
0141     LOG(LEVEL) << " ekey " << ekey << " Desc " << Desc(cegs)  ;
0142 }
0143 
0144 std::string SFrameGenstep::GetGridConfig(std::vector<int>& cegs)
0145 {
0146     const char* CEGS_fallback = CEGS_XZ ;
0147     GetGridConfig_(cegs, CEGS, CEGS_delim, CEGS_fallback );
0148 
0149     std::stringstream ss ;
0150     assert( cegs.size() == 8 );
0151     for(int i=0 ; i < 8 ; i++ ) ss << cegs[i] << ( i < 8 - 1 ? "," : "" ) ;
0152     std::string str = ss.str();
0153     return str ;
0154 }
0155 
0156 
0157 
0158 
0159 
0160 /**
0161 SFrameGenstep::MakeCenterExtentGensteps_FromFrame
0162 --------------------------------------------------
0163 
0164 Canonical usage is from SEvt::createInputGenstep_simtrace
0165 
0166 Default is to apply extent scaling.
0167 To switch that off (eg whilse debugging) use::
0168 
0169     export CE_SCALE_OFF=1
0170 
0171 Default GRIDSCALE is 0.1 to conform with default CEGS of 16:0:9:500
0172 
0173 The sframe argument is used for its *ce* and the details
0174 of the grid config are saved into the sframe
0175 which makes the config available from python.
0176 
0177 **/
0178 
0179 const char* SFrameGenstep::CEGS_XY = "16:9:0:1000" ;
0180 const char* SFrameGenstep::CEGS_XZ = "16:0:9:1000" ;  // default
0181 
0182 
0183 bool SFrameGenstep::HasConfigEnv()
0184 {
0185     return ssys::hasenv_("CEGS");
0186 }
0187 
0188 
0189 
0190 /**
0191 SFrameGenstep::MakeCenterExtentGenstep_FromFrame
0192 --------------------------------------------------
0193 
0194 This uses config obtained from envvars such as CEGS to
0195 create an array of gensteps. sframe::set_grid is
0196 invoked to pass config details into the frame for
0197 persisting and use from python such as simtrace
0198 plotter cxt_min.py. Those details are also
0199 written into metadata of the gensteps.
0200 
0201 HMM: could avoid changing the frame by instead
0202 relying on the genstep metadata which might
0203 simplify FRAME_TRANSITION
0204 
0205 Q: What does simtrace cxt_min.py actually need from the frame ?
0206 
0207 
0208 simtrace stack::
0209 
0210     SFrameGenstep::MakeCenterExtentGenstep_FromFrame
0211     SEvt::createInputGenstep_simtrace
0212     SEvt::createInputGenstep_configured
0213     SEvt::addInputGenstep
0214     SEvt::beginOfEvent
0215     QSim::simtrace
0216     CSGOptiX::simtrace
0217     CSGOptiX::SimtraceMain
0218     main
0219 
0220 
0221 **/
0222 #ifdef WITH_OLD_FRAME
0223 NP* SFrameGenstep::MakeCenterExtentGenstep_FromFrame(sframe& fr)  // static
0224 {
0225     const float4& ce = fr.ce ;
0226     Tran<double>* geotran = fr.getTransform();
0227     char* _GRIDSCALE = getenv("GRIDSCALE") ;
0228     float gridscale = ssys::getenvfloat("GRIDSCALE", 0.1f ) ;
0229     int prim = -1 ;
0230 
0231     LOG(LEVEL)
0232        << " _GRIDSCALE [" << ( _GRIDSCALE ? _GRIDSCALE : "-" ) << "]"
0233        << " GRIDSCALE " << gridscale
0234        ;
0235 
0236     std::vector<int> cegs ;
0237     std::string str_cegs = GetGridConfig(cegs);
0238     fr.set_grid(cegs, gridscale);
0239 
0240     NP* gs = MakeCenterExtentGenstep_From_CE_geotran( ce, cegs, gridscale, geotran, prim );
0241 
0242     gs->set_meta<int>("midx", fr.midx() );
0243     gs->set_meta<int>("mord", fr.mord() );
0244     gs->set_meta<int>("gord", fr.gord() );
0245     gs->set_meta<float>("gridscale", fr.gridscale() );
0246     gs->set_meta<std::string>("cegs", str_cegs );
0247 
0248     return gs ;
0249 }
0250 #else
0251 NP* SFrameGenstep::MakeCenterExtentGenstep_FromFrame(sfr& fr)  // static
0252 {
0253     float4 ce = make_float4( fr.ce.x, fr.ce.y, fr.ce.z, fr.ce.w );
0254     Tran<double>* geotran = fr.getTransform();
0255     char* _GRIDSCALE = getenv("GRIDSCALE") ;
0256     float gridscale = ssys::getenvfloat("GRIDSCALE", 0.1f ) ;
0257     int prim = fr.get_prim() ;
0258 
0259     LOG(LEVEL)
0260        << " _GRIDSCALE [" << ( _GRIDSCALE ? _GRIDSCALE : "-" ) << "]"
0261        << " GRIDSCALE " << gridscale
0262        ;
0263 
0264 
0265     std::vector<int> cegs ;
0266     std::string str_cegs = GetGridConfig(cegs);
0267     fr.set_gridscale(gridscale);
0268 
0269     LOG(LEVEL)
0270        << " cegs.size " << cegs.size()
0271        ;
0272 
0273 
0274     NP* gs = MakeCenterExtentGenstep_From_CE_geotran( ce, cegs, gridscale, geotran, prim );
0275 
0276     gs->set_meta<float>("gridscale", fr.get_gridscale() );
0277     gs->set_meta<std::string>("cegs", str_cegs );
0278 
0279     return gs ;
0280 }
0281 #endif
0282 
0283 NP* SFrameGenstep::MakeCenterExtentGenstep_From_CE_geotran(const float4& ce, const std::vector<int>& cegs, float gridscale, const Tran<double>* geotran, int prim)  // static
0284 {
0285     std::vector<float>* cegs_radial_range = ssys::getenv_vec<float>(CEGS_RADIAL_RANGE, nullptr, CEGS_delim );
0286 
0287     std::vector<float3> ce_offset ;
0288     CE_OFFSET(ce_offset, ce);
0289 
0290     LOG(LEVEL)
0291         << " ce " << ce
0292         << " ce_offset.size " << ce_offset.size()
0293         ;
0294 
0295     bool with_offset = ce_offset.size() > 0 ;
0296     LOG_IF(fatal, !with_offset) << "ce_offset vector of float3 needs at least one entry " ;
0297     assert( with_offset );
0298 
0299     assert(!getenv("CE_SCALE") && "CE_SCALE ENVVAR IS NO LONGER USED :CHANGE YOUR SCRIPT" );
0300 
0301     bool ce_scale_off = ssys::getenvbool("CE_SCALE_OFF") ;
0302 
0303 
0304 
0305     std::vector<NP*> gsl ;
0306     NP* gs_base = MakeCenterExtentGenstep(ce, cegs, gridscale, geotran, ce_offset, !ce_scale_off, cegs_radial_range );
0307     gsl.push_back(gs_base) ;
0308 
0309     std::vector<std::string> keys = {{
0310         "CEHIGH_0",
0311         "CEHIGH_1",
0312         "CEHIGH_2",
0313         "CEHIGH_3",
0314         "CEHIGH_4",
0315         "CEHIGH_5",
0316         "CEHIGH_6",
0317         "CEHIGH_7",
0318         "CEHIGH_8",
0319         "CEHIGH_9"
0320       }} ;
0321 
0322     const char* CEHIGH_fallback = "" ;
0323     // need the envvar CEHIGH_0,... for the below to add gensteps
0324 
0325     for(unsigned i=0 ; i < keys.size() ; i++)
0326     {
0327         const char* key = keys[i].c_str() ;
0328         std::vector<int> cehigh ;
0329         GetGridConfig_(cehigh, key, CEGS_delim, CEHIGH_fallback );
0330         LOG(LEVEL) << " key " << key << " cehigh.size " << cehigh.size() ;
0331         if(cehigh.size() == 8)
0332         {
0333             NP* gs_cehigh = MakeCenterExtentGenstep(ce, cehigh, gridscale, geotran, ce_offset, !ce_scale_off, cegs_radial_range );
0334             gsl.push_back(gs_cehigh) ;
0335         }
0336     }
0337 
0338 
0339     NP* CEGS_NPY = ssys::hasenv_("CEGS_NPY") ? NP::Load("$CEGS_NPY") : nullptr ;
0340     NP* gs_CEGS_NPY = CEGS_NPY ? Make_CEGS_NPY_Genstep(CEGS_NPY, geotran) : nullptr ;
0341     if(gs_CEGS_NPY) gsl.push_back(gs_CEGS_NPY);
0342 
0343 
0344     Maybe_Add_PRIOR_SIMTRACE_Genstep( gsl, prim );
0345 
0346     LOG(LEVEL) << " gsl.size " << gsl.size() ;
0347 
0348     LOG(LEVEL) << "[ NP::Concatenate " ;
0349     NP* gs = NP::Concatenate(gsl) ;
0350     LOG(LEVEL) << "] NP::Concatenate " ;
0351 
0352     gs->set_meta<int>("ce_scale", int(!ce_scale_off) );
0353 
0354     return gs ;
0355 }
0356 
0357 
0358 /**
0359 SFrameGenstep::Maybe_Add_PRIOR_SIMTRACE_Genstep
0360 ------------------------------------------------
0361 
0362 PRIOR_SIMTRACE
0363     envvar which when defined causes loading of the array
0364     the intersect positions from which are used to define
0365     genstep positions
0366 
0367 
0368 PRIOR_SIMTRACE_PRIM:-2
0369     default, use target prim from frame to select the
0370     intersects with which to base the new gensteps
0371 
0372 PRIOR_SIMTRACE_PRIM:-1
0373     do not select from the prior simtrace, use them all
0374     without checking the prim. This is appropriate when
0375     PRIOR_SIMTRACE is $MFOLD/simtrace_overlap.npy
0376 
0377 PRIOR_SIMTRACE_PRIM:0,1,2,3,
0378     use the provided prim index (not typically used)
0379 
0380 
0381 Examples::
0382 
0383     export PRIOR_SIMTRACE=$MFOLD/simtrace.npy
0384     export PRIOR_SIMTRACE_PRIM=-2  # default
0385 
0386     export PRIOR_SIMTRACE=$MFOLD/simtrace_overlap.npy
0387     export PRIOR_SIMTRACE_PRIM=-1  # use all without selection
0388 
0389 **/
0390 
0391 
0392 void SFrameGenstep::Maybe_Add_PRIOR_SIMTRACE_Genstep( std::vector<NP*>& gsl, int prim )
0393 {
0394     bool with_PRIOR_SIMTRACE = ssys::hasenv_("PRIOR_SIMTRACE") ;
0395     NP* PRIOR_SIMTRACE = with_PRIOR_SIMTRACE ? NP::Load("$PRIOR_SIMTRACE") : nullptr ;
0396 
0397     int PRIOR_SIMTRACE_PRIM = ssys::getenvint("PRIOR_SIMTRACE_PRIM", -2 );
0398     int uprim = prim ;
0399     if( PRIOR_SIMTRACE_PRIM == -2 )
0400     {
0401         uprim = prim ;  // select prior intersects onto target prim
0402     }
0403     else if( PRIOR_SIMTRACE_PRIM == -1 )
0404     {
0405         uprim = -1 ;    // use all prior intersects without prim selection
0406     }
0407     else if( PRIOR_SIMTRACE_PRIM > -1 )
0408     {
0409         uprim = PRIOR_SIMTRACE_PRIM ;
0410     }
0411 
0412     NP* gs_PRIOR_SIMTRACE = PRIOR_SIMTRACE  ?  Make_PRIOR_SIMTRACE_Genstep( PRIOR_SIMTRACE, uprim ) : nullptr ;
0413     if(gs_PRIOR_SIMTRACE) gsl.push_back(gs_PRIOR_SIMTRACE);
0414     LOG(info)
0415          << "with_PRIOR_SIMTRACE " << ( with_PRIOR_SIMTRACE ? "YES" : "NO " )
0416          << " prim " << prim
0417          << " uprim " << uprim
0418          << " PRIOR_SIMTRACE_PRIM " << PRIOR_SIMTRACE_PRIM
0419          << " PRIOR_SIMTRACE " << ( PRIOR_SIMTRACE ? PRIOR_SIMTRACE->sstr() : "-" )
0420          << " gs_PRIOR_SIMTRACE " << ( gs_PRIOR_SIMTRACE ? gs_PRIOR_SIMTRACE->sstr() : "-" )
0421          ;
0422 }
0423 
0424 /**
0425 SFrameGenstep::Make_PRIOR_SIMTRACE_Genstep
0426 -------------------------------------------
0427 
0428 **/
0429 
0430 
0431 NP* SFrameGenstep::Make_PRIOR_SIMTRACE_Genstep( const NP* _simtrace, int prim ) // static
0432 {
0433     const NP* _ontoprim = prim > -1 ? quad4::select_prim( _simtrace, prim ) : _simtrace ;
0434     int ni = _ontoprim->shape[0] ;
0435     const quad4* ontoprim = (const quad4*)(_ontoprim->bytes());
0436 
0437     quad6 gs ;
0438     gs.zero();
0439 
0440     glm::tmat4x4<double> identity(1.);
0441     qat4* qc = Tran<double>::ConvertFrom( identity ) ;
0442 
0443     std::vector<quad6> gensteps ;
0444 
0445     float dx = 0.f;
0446     float dy = 0.f;
0447     float dz = 0.f;
0448 
0449     float s = 1.f ;
0450 
0451     for(int i=0 ; i < ni ; i++)
0452     {
0453         const quad4& q = ontoprim[i];
0454         float x = q.q1.f.x ;
0455         float y = q.q1.f.y ;
0456         float z = q.q1.f.z ;
0457 
0458         for(int j=0 ; j < 6 ; j++)
0459         {
0460             switch(j)
0461             {
0462                case 0: dx = -s  ; dy = 0.f ; dz = 0.f ; break ;
0463                case 1: dx = +s  ; dy = 0.f ; dz = 0.f ; break ;
0464                case 2: dx = 0.f ; dy = -s  ; dz = 0.f ; break ;
0465                case 3: dx = 0.f ; dy = +s  ; dz = 0.f ; break ;
0466                case 4: dx = 0.f ; dy = 0.f ; dz = -s  ; break ;
0467                case 5: dx = 0.f ; dy = 0.f ; dz =  s  ; break ;
0468             }
0469 
0470             gs.q1.f.x = x + dx ;
0471             gs.q1.f.y = y + dy ;
0472             gs.q1.f.z = z + dz ;
0473             gs.q1.f.w = 1.f ;
0474 
0475             qc->write(gs);  // copy qc transform into gs.q2,q3,q4,q5
0476 
0477             int gridaxes = 0 ;
0478             int gsid = 0 ;
0479             int photons_per_genstep = 100 ;
0480             SGenstep::ConfigureGenstep(gs, OpticksGenstep_FRAME, gridaxes, gsid, photons_per_genstep );
0481 
0482             gensteps.push_back(gs);
0483         }
0484     }
0485 
0486     return SGenstep::MakeArray(gensteps);
0487 }
0488 
0489 
0490 /**
0491 SFrameGenstep::Make_CEGS_NPY_Genstep
0492 -------------------------------------
0493 
0494 Construct gensteps from an array of positions of interest of shape (-1,3)
0495 
0496 **/
0497 
0498 
0499 NP* SFrameGenstep::Make_CEGS_NPY_Genstep( const NP* t, const Tran<double>* geotran ) // static
0500 {
0501     assert( t && t->has_shape(-1,3) && t->uifc == 'f' && t->ebyte == 4 );
0502     int ni = t->shape[0] ;
0503     int nj = t->shape[1] ;
0504     assert( nj == 3 );
0505 
0506     const float* tt = t->cvalues<float>();
0507 
0508     qat4* qc = Tran<double>::ConvertFrom( geotran->t ) ;
0509 
0510     quad6 gs ;
0511     gs.zero();
0512 
0513     std::vector<quad6> gensteps ;
0514 
0515     for(int i=0 ; i < ni ; i++)
0516     {
0517         gs.q1.f.x = tt[nj*i+0] ;
0518         gs.q1.f.y = tt[nj*i+1] ;
0519         gs.q1.f.z = tt[nj*i+2] ;
0520         gs.q1.f.w = 1.f ;
0521 
0522         qc->write(gs);  // copy qc transform into gs.q2,q3,q4,q5
0523 
0524         int gridaxes = 0 ;
0525         int gsid = 0 ;
0526         int photons_per_genstep = 100 ;
0527         SGenstep::ConfigureGenstep(gs, OpticksGenstep_FRAME, gridaxes, gsid, photons_per_genstep );
0528 
0529         gensteps.push_back(gs);
0530     }
0531     return SGenstep::MakeArray(gensteps);
0532 }
0533 
0534 
0535 
0536 
0537 
0538 /**
0539 SFrameGenstep::MakeCenterExtentGenstep
0540 -----------------------------------------
0541 
0542 Creates grid of gensteps centered at ce.xyz with the grid specified
0543 by integer ranges that are used to scale the extent parameter to yield
0544 offsets from the center.
0545 
0546 ce(float4)
0547    cx:cy:cz:extent
0548 
0549 cegs(uint4)
0550    nx:ny:nz:photons_per_genstep
0551    specifies a grid of integers -nx:nx -ny:ny -nz:nz inclusive used to scale the extent
0552 
0553    The number of gensteps becomes: (2*nx+1)*(2*ny+1)*(2*nz+1)
0554 
0555 gridscale
0556    float multiplier applied to the grid integers, values less than 1. (eg 0.2)
0557    increase the concentration of the genstep grid on the target geometry giving a
0558    better intersect rendering of a smaller region
0559 
0560    To expand the area when using a finer grid increase the nx:ny:nz, however
0561    that will lead to a slower render.
0562 
0563 ce_offset:vector of float3
0564    typically local frame origin (0.,0.,0.)
0565    the grid is repeated for each offset float3
0566    This is used for example to create intersects in multiple planes.
0567 
0568 ce_scale:true
0569    grid translate offsets *local_scale* set to ce.w*gridscale
0570 
0571    SEEMS LIKE THIS SHOULD ALWAYS BE USED ?
0572    ACTUALLTY NOT : APPARENTLY SOME TYPES OF TRANSFORM INCORPORATE THE EXTENT
0573    SCALE ALREADY
0574 
0575 
0576 ce_scale:false
0577    grid translate offsets *local_scale* set to gridscale
0578 
0579 
0580 
0581 The gensteps are consumed by qsim::generate_photon_torch
0582 Which needs to use the gensteps data in order to transform the axis
0583 aligned local frame grid of positions and directions
0584 into global frame equivalents.
0585 
0586 Instance transforms are best regarded as first doing rotate
0587 about a local origin and then translate into global position.
0588 When wish to create multiple transforms with small local frame offsets
0589 to create a grid or plane between them need to first pre-multiply by the
0590 small local translation followed by the rotation and large global translation
0591 into position.
0592 
0593 For example when using reverse=true get individual tilts out of the plane
0594 ie the single local XZ plane becomes lots of planes in global frame as the local_translate is done last.
0595 When using reverse=false get all the tilts the same so local XZ single plane stays one plane in the global
0596 frame as are doing the local_translate first.
0597 
0598 The *high = cegs[7]* parameter is normally 1.  For high resolution regions it is 2/3/4
0599 which multiplies the grid points and divides the gridscale.
0600 This allows high resolution regions without changes to indexing.
0601 
0602 **/
0603 
0604 NP* SFrameGenstep::MakeCenterExtentGenstep(
0605     const float4& ce,
0606     const std::vector<int>& cegs,
0607     float gridscale,
0608     const Tran<double>* geotran,
0609     const std::vector<float3>& ce_offset,
0610     bool ce_scale,
0611     std::vector<float>* cegs_radial_range) // static
0612 {
0613 
0614     quad6 gs ;
0615     gs.zero();
0616 
0617     assert( cegs.size() == 8 );
0618 
0619     int high = cegs[7] ;  // > 1 scales grid resolution
0620     assert( high >= 1 && high <= 8);
0621     double scale = double(gridscale)/double(high) ;
0622 
0623     int ix0 = cegs[0]*high ;
0624     int ix1 = cegs[1]*high ;
0625     int iy0 = cegs[2]*high ;
0626     int iy1 = cegs[3]*high ;
0627     int iz0 = cegs[4]*high ;
0628     int iz1 = cegs[5]*high ;
0629 
0630     int photons_per_genstep = cegs[6] ;
0631 
0632     int nx = (ix1 - ix0)/2 ;
0633     int ny = (iy1 - iy0)/2 ;
0634     int nz = (iz1 - iz0)/2 ;
0635 
0636     int gridaxes = SGenstep::GridAxes(nx, ny, nz);   // { XYZ, YZ, XZ, XY }
0637     int num_offset = int(ce_offset.size()) ;
0638 
0639     LOG(LEVEL)
0640         << " num_offset " << num_offset
0641         << " ce_scale " << ce_scale
0642         << " nx " << nx
0643         << " ny " << ny
0644         << " nz " << nz
0645         << " GridAxes " << gridaxes
0646         << " GridAxesName " << SGenstep::GridAxesName(gridaxes)
0647         << " high " << high
0648         << " gridscale " << gridscale
0649         << " scale " << scale
0650         << " photons_per_genstep " << photons_per_genstep
0651         ;
0652 
0653     double local_scale = ce_scale ? scale*ce.w : scale ; // ce_scale:true is almost always expected
0654     // hmm: when using SCenterExtentFrame model2world transform the
0655     // extent is already handled within the transform so must not apply extent scaling
0656     // THIS IS CONFUSING : TODO FIND WAY TO AVOID THE CONFUSION BY MAKING THE DIFFERENT TYPES OF TRANSFORM MORE CONSISTENT
0657 
0658     [[maybe_unused]] unsigned photon_offset = 0 ;
0659     std::vector<quad6> gensteps ;
0660 
0661     float rmin = 0. ;
0662     float rmax = std::numeric_limits<float>::max();
0663     if(cegs_radial_range)
0664     {
0665         size_t crr = cegs_radial_range->size();
0666         if(crr > 0) rmin = (*cegs_radial_range)[0] ;
0667         if(crr > 1) rmax = (*cegs_radial_range)[1] ;
0668 
0669         LOG(info)
0670             << " [" << CEGS_RADIAL_RANGE << "] "
0671             << " crr " << crr
0672             << " rmin " << std::fixed << std::setw(10) << std::setprecision(2) << rmin
0673             << " rmax " << std::fixed << std::setw(10) << std::setprecision(2) << rmax
0674             << " ce.w " << std::fixed << std::setw(10) << std::setprecision(2) << ce.w
0675             << " scale " << std::fixed << std::setw(10) << std::setprecision(2) << scale
0676             << " local_scale " << std::fixed << std::setw(10) << std::setprecision(2) << local_scale
0677             << "\n"
0678             ;
0679     }
0680 
0681 
0682 
0683     for(int ip=0 ; ip < num_offset ; ip++)   // planes
0684     {
0685         const float3& offset = ce_offset[ip] ;
0686 
0687         gs.q1.f.x = offset.x ;
0688         gs.q1.f.y = offset.y ;
0689         gs.q1.f.z = offset.z ;
0690         gs.q1.f.w = 1.f ;
0691 
0692         for(int ix=ix0 ; ix < ix1+1 ; ix++ )
0693         for(int iy=iy0 ; iy < iy1+1 ; iy++ )
0694         for(int iz=iz0 ; iz < iz1+1 ; iz++ )
0695         {
0696             double tx = double(ix)*local_scale ;
0697             double ty = double(iy)*local_scale ;
0698             double tz = double(iz)*local_scale ;
0699 
0700             double radius = std::sqrt(tx*tx + ty*ty + tz*tz);
0701             bool in_radial_range = cegs_radial_range ?  ( radius >= rmin && radius <= rmax ) : true ;
0702             if(!in_radial_range) continue ;
0703 
0704 
0705 
0706             const Tran<double>* grid_shift = Tran<double>::make_translate( tx, ty, tz );
0707 
0708             bool reverse = false ;
0709             const Tran<double>* transform = Tran<double>::product( geotran, grid_shift, reverse );
0710 
0711             qat4* qc = Tran<double>::ConvertFrom( transform->t ) ;
0712 
0713             unsigned gsid = SGenstep::GenstepID(ix,iy,iz,ip) ;
0714 
0715             SGenstep::ConfigureGenstep(gs, OpticksGenstep_FRAME, gridaxes, gsid, photons_per_genstep );
0716 
0717             qc->write(gs);  // copy qc into gs.q2,q3,q4,q5
0718 
0719             gensteps.push_back(gs);
0720             photon_offset += std::abs(photons_per_genstep) ;
0721         }
0722     }
0723 
0724     LOG(LEVEL)
0725          << " num_offset " << num_offset
0726          << " gensteps.size " << gensteps.size()
0727          ;
0728     return SGenstep::MakeArray(gensteps);
0729 }
0730 
0731 
0732 
0733 /**
0734 SFrameGenstep::StandardizeCEGS
0735 --------------------------------
0736 
0737 The cegs vector configures a grid.
0738 Symmetric and offset grid input configs are supported using
0739 vectors of length 4 and 7.
0740 
0741 This method standardizes the specification
0742 into an absolute index form which is used by
0743 SFrameGenstep::MakeCenterExtentGensteps
0744 
0745 nx:ny:nz:num_photons
0746      symmetric grid -nx:nx, -ny:ny, -nz:nz
0747      with total of (2*nx+1)*(2*ny+1)*(2*nz+1)*num_photons
0748      eg with 16:0:9:2000 thats
0749      (2*16+1)*(2*9+1)*2000 = 1,254,000
0750 
0751 nx:ny:nz:dx:dy:dz:num_photons
0752      offset grid -nx+dx:nx+dx, -ny+dy:ny+dy, -nz+dz:nz+dz
0753      with total still (2*nx+1)*(2*ny+1)*(2*nz+1)*num_photons
0754 
0755      To land a grid above in Z to 16:0:9:2000 use 16:0:9:0:0:10:2000
0756 
0757 ix0:ix1:iy0:iy1:iz0:iz1:num_photons:high
0758      standardized absolute form of grid specification
0759      that SFrameGenstep::StandardizeCEGS creates.
0760 
0761      This eight element form is also now used as an
0762      input form when wishing to increase resolution with
0763      the "high" 8th element cegs[7].
0764 
0765 Examples of CEHIGH envvars using the third standardized form::
0766 
0767     export CEHIGH_0=-11:-9:0:0:-2:0:1000:4
0768     export CEHIGH_1=9:11:0:0:-2:0:1000:4
0769     export CEHIGH_2=-1:1:0:0:-2:0:1000:4
0770 
0771 
0772 
0773 For example a portrait base CEGS of 9:0:16:1000 means,
0774 becomes -9:9:0:0:-16:16:1000:1
0775 
0776 
0777                                16
0778                                15
0779                                14
0780                                13
0781                                12
0782                                11
0783                                10
0784                                 9
0785                                 8
0786                                 7
0787                                 6
0788                                 5
0789                                 4
0790                                 3
0791                                 2
0792                                 1
0793      9  8  7  6  5  4  3  2  1  0  1  2  3  4  5  6  7  8  9
0794                                 1
0795                                 2
0796                                 3
0797                     +           4           +
0798                                 5
0799                                 6
0800                                 7
0801                                 8
0802                                 9
0803                                10
0804                                11
0805                     +          12           +
0806                                13
0807                                14
0808                                15
0809                                16
0810 
0811 
0812 
0813 To quadruple grid resolution in the above square::
0814 
0815       -4:4:0:0:-12:-4:1000:4
0816 
0817 
0818 **/
0819 
0820 void SFrameGenstep::StandardizeCEGS( std::vector<int>& cegs ) // static
0821 {
0822     if( cegs.size() == 0 ) return ;
0823     if( cegs.size() == 8 ) return ;  // when start with 8, just exit assuming already standardized
0824 
0825     bool expect = cegs.size() == 4 || cegs.size() == 7 ;
0826     LOG_IF(error, !expect) << " unexpected cegs.size " << cegs.size() ;
0827     assert( expect );
0828 
0829     int ix0 = 0 ;
0830     int ix1 = 0 ;
0831     int iy0 = 0 ;
0832     int iy1 = 0 ;
0833     int iz0 = 0 ;
0834     int iz1 = 0 ;
0835     int photons_per_genstep = 0 ;
0836 
0837     if( cegs.size() == 4 )
0838     {
0839         ix0 = -cegs[0] ; ix1 = cegs[0] ;
0840         iy0 = -cegs[1] ; iy1 = cegs[1] ;
0841         iz0 = -cegs[2] ; iz1 = cegs[2] ;
0842         photons_per_genstep = cegs[3] ;
0843     }
0844     else if( cegs.size() == 7 )
0845     {
0846         int nx = std::abs(cegs[0]) ;
0847         int ny = std::abs(cegs[1]) ;
0848         int nz = std::abs(cegs[2]) ;
0849         int dx = cegs[3] ;
0850         int dy = cegs[4] ;
0851         int dz = cegs[5] ;
0852         photons_per_genstep = cegs[6] ;
0853 
0854         ix0 = -nx + dx ;
0855         iy0 = -ny + dy ;
0856         iz0 = -nz + dz ;
0857         ix1 =  nx + dx ;
0858         iy1 =  ny + dy ;
0859         iz1 =  nz + dz ;
0860     }
0861 
0862     cegs.resize(8) ;
0863     cegs[0] = ix0 ;
0864     cegs[1] = ix1 ;
0865     cegs[2] = iy0 ;
0866     cegs[3] = iy1 ;
0867     cegs[4] = iz0 ;
0868     cegs[5] = iz1 ;
0869     cegs[6] = photons_per_genstep ;
0870     cegs[7] = 1 ;   // high
0871 
0872     //  +---+---+---+---+
0873     // -2  -1   0   1   2
0874 
0875     unsigned grid_points = (ix1-ix0+1)*(iy1-iy0+1)*(iz1-iz0+1) ;
0876     unsigned tot_photons = grid_points*std::abs(photons_per_genstep) ;
0877 
0878     LOG(LEVEL)
0879         << " CEGS "
0880         << " ix0 ix1 " << ix0 << " " << ix1
0881         << " iy0 iy1 " << iy0 << " " << iy1
0882         << " iz0 iz1 " << iz0 << " " << iz1
0883         << " photons_per_genstep " << photons_per_genstep
0884         << " grid_points (ix1-ix0+1)*(iy1-iy0+1)*(iz1-iz0+1) " << grid_points
0885         << " tot_photons (grid_points*photons_per_genstep) " << tot_photons
0886         ;
0887 }
0888 
0889 
0890 /**
0891 SFrameGenstep::GetBoundingBox
0892 -------------------------------
0893 
0894 Uses CE center-extent and gridscale together with the cegs grid parameters to provide
0895 the float3 mn and mx bounds of the CE grid.  NB no use of any transforms here.
0896 
0897 **/
0898 
0899 void SFrameGenstep::GetBoundingBox( float3& mn, float3& mx, const float4& ce, const std::vector<int>& standardized_cegs, float gridscale, const float3& ce_offset ) // static
0900 {
0901     assert( standardized_cegs.size() == 7 ) ;
0902 
0903     int ix0 = standardized_cegs[0] ;
0904     int ix1 = standardized_cegs[1] ;
0905     int iy0 = standardized_cegs[2] ;
0906     int iy1 = standardized_cegs[3] ;
0907     int iz0 = standardized_cegs[4] ;
0908     int iz1 = standardized_cegs[5] ;
0909     int photons_per_genstep = standardized_cegs[6] ;
0910 
0911 
0912     float x0 = float(ix0)*gridscale*ce.w + ce_offset.x ;   // ce_offset is usually local fram origin (0.f,0.f,0.f)
0913     float x1 = float(ix1)*gridscale*ce.w + ce_offset.x ;
0914 
0915     float y0 = float(iy0)*gridscale*ce.w + ce_offset.y ;
0916     float y1 = float(iy1)*gridscale*ce.w + ce_offset.y ;
0917 
0918     float z0 = float(iz0)*gridscale*ce.w + ce_offset.z ;
0919     float z1 = float(iz1)*gridscale*ce.w + ce_offset.z ;
0920 
0921     mn.x = x0 ;
0922     mx.x = x1 ;
0923 
0924     mn.y = y0 ;
0925     mx.y = y1 ;
0926 
0927     mn.z = z0 ;
0928     mx.z = z1 ;
0929 
0930     LOG(LEVEL)
0931         << " ce_offset " << ce_offset
0932         << " x0 " << std::setw(10) << std::fixed << std::setprecision(3) << x0
0933         << " x1 " << std::setw(10) << std::fixed << std::setprecision(3) << x1
0934         << " y0 " << std::setw(10) << std::fixed << std::setprecision(3) << y0
0935         << " y1 " << std::setw(10) << std::fixed << std::setprecision(3) << y1
0936         << " z0 " << std::setw(10) << std::fixed << std::setprecision(3) << z0
0937         << " z1 " << std::setw(10) << std::fixed << std::setprecision(3) << z1
0938         << " photons_per_genstep " << photons_per_genstep
0939         << " gridscale " << std::setw(10) << std::fixed << std::setprecision(3) << gridscale
0940         << " ce.w(extent) " << std::setw(10) << std::fixed << std::setprecision(3) << ce.w
0941         ;
0942 }
0943 
0944 
0945 
0946 
0947 
0948 
0949 
0950 
0951 
0952 
0953 
0954 
0955 
0956 /**
0957 SFrameGenstep::GenerateCenterExtentGenstepPhotons
0958 -----------------------------------------------------
0959 
0960 **/
0961 
0962 NP* SFrameGenstep::GenerateCenterExtentGenstepPhotons_( const NP* gsa, float gridscale )
0963 {
0964     std::vector<quad4> pp ;
0965     GenerateCenterExtentGenstepPhotons( pp, gsa, gridscale );
0966 
0967     NP* ppa = NP::Make<float>( pp.size(), 4, 4 );
0968     memcpy( ppa->bytes(),  (float*)pp.data(), ppa->arr_bytes() );
0969     return ppa ;
0970 }
0971 
0972 
0973 
0974 
0975 /**
0976 SFrameGenstep::GenerateCenterExtentGenstepPhotons
0977 ---------------------------------------------------
0978 
0979 Contrast this CPU implementation of CEGS generation with qudarap/qsim.h qsim<T>::generate_photon_torch
0980 
0981 The inner loop should be similar to qudarap/qsim.h/generate_photon_simtrace
0982 
0983 TODO: arrange a header such that can actually use the same code via some curand_uniform
0984 macro refinition trickery
0985 
0986 -ve num_photons uses regular, not random azimuthal spray of directions
0987 
0988 TODO: consolidate the below two methods, paradir is the hangup
0989 
0990 
0991 **/
0992 
0993 void SFrameGenstep::GenerateCenterExtentGenstepPhotons( std::vector<quad4>& pp, const NP* gsa, float gridscale )
0994 {
0995     LOG(LEVEL) << " gsa " << gsa->sstr() ;
0996 
0997     assert( gsa->shape.size() == 3 && gsa->shape[1] == 6 && gsa->shape[2] == 4 );
0998     assert( gsa->has_shape(-1,6,4) );
0999 
1000     std::vector<quad6> gsv(gsa->shape[0]) ;
1001     memcpy( gsv.data(), gsa->bytes(), gsa->arr_bytes() );
1002 
1003     quad4 p ;
1004     p.zero();
1005 
1006     unsigned seed = 0 ;
1007     SRng<float> rng(seed) ;
1008 
1009 
1010     float3 paradir ;
1011     qvals(paradir, "PARADIR", "0,0,0" );
1012     bool with_paradir = dot(paradir,paradir) > 0.f ;
1013     if(with_paradir)
1014     {
1015         paradir = normalize(paradir);
1016         LOG(LEVEL) << " PARADIR enabled " << paradir ;
1017     }
1018     else
1019     {
1020         LOG(LEVEL) << " PARADIR NOT-enabled " ;
1021     }
1022 
1023 
1024     for(unsigned i=0 ; i < gsv.size() ; i++)
1025     {
1026         const quad6& gs = gsv[i];
1027         qat4 qt(gs) ;  // copy 4x4 transform from last 4 quads of genstep
1028 
1029         C4U gsid ;   // genstep integer grid coordinate IXYZ and IW photon index up to 255
1030 
1031         int gencode           = gs.q0.i.x ;
1032         int gridaxes          = gs.q0.i.y ;
1033         gsid.u                = gs.q0.u.z ;     // formerly gs.q1.u.w
1034         int      num_photons_ = gs.q0.i.w ;     // -ve num_photons_ doesnt use random u0 so creates phi wheels
1035 
1036         unsigned num_photons  = std::abs(num_photons_);
1037 
1038         bool expect = gencode == OpticksGenstep_TORCH || gencode == OpticksGenstep_FRAME ;
1039 
1040         LOG_IF(error, !expect)
1041             << "unexpected gencod " << gencode
1042             << " OpticksGenstep_::Name " << OpticksGenstep_::Name(gencode) ;
1043 
1044         assert(expect);
1045 
1046         //std::cout << " i " << i << " num_photons " << num_photons << std::endl ;
1047 
1048         double u0, u1 ;
1049         double phi, sinPhi,   cosPhi ;
1050         double sinTheta, cosTheta ;
1051 
1052 
1053 
1054         for(unsigned j=0 ; j < num_photons ; j++)
1055         {
1056             u0 = num_photons_ < 0 ? double(j)/double(num_photons-1) : rng() ;
1057 
1058             phi = 2.*M_PIf*u0 ;     // azimuthal 0->2pi
1059             ssincos(phi,sinPhi,cosPhi);
1060 
1061 
1062             // cosTheta sinTheta are only used for 3D (not 2D planar gensteps)
1063             u1 = rng();
1064             cosTheta = u1 ;
1065             sinTheta = sqrtf(1.0-u1*u1);
1066 
1067             if( with_paradir )
1068             {
1069                 // what scaling is needed to span the grid ?
1070                 p.q0.f.x = 0.f ;
1071                 p.q0.f.y = u0*float(num_photons-1)*gridscale ;
1072                 p.q0.f.z = 0.f ;
1073                 p.q0.f.w = 1.f ;
1074 
1075                 p.q1.f.x = paradir.x ;
1076                 p.q1.f.y = paradir.y ;
1077                 p.q1.f.z = paradir.z ;
1078                 p.q1.f.w = 0.f ;
1079             }
1080             else
1081             {
1082                 // copy position from genstep into the photon, historically has been origin
1083                 p.q0.f.x = gs.q1.f.x ;
1084                 p.q0.f.y = gs.q1.f.y ;
1085                 p.q0.f.z = gs.q1.f.z ;
1086                 p.q0.f.w = 1.f ;       // <-- dont copy the "float" gsid
1087 
1088                 SetGridPlaneDirection( p.q1.f, gridaxes, cosPhi, sinPhi, cosTheta, sinTheta );
1089             }
1090 
1091             // tranform photon position and direction into the desired frame
1092             qt.right_multiply_inplace( p.q0.f, 1.f );  // position
1093             qt.right_multiply_inplace( p.q1.f, 0.f );  // direction
1094 
1095             unsigned char ucj = (j < 255 ? j : 255 ) ;  // photon index local to the genstep
1096             gsid.c4.w = ucj ;     // setting C4U union element to change gsid.u
1097             p.q3.u.w = gsid.u ;   // include photon index IW with the genstep coordinate into photon (3,3)
1098 
1099             pp.push_back(p) ;
1100         }
1101     }
1102 
1103     LOG(LEVEL) << " pp.size " << pp.size() ;
1104 }
1105 
1106 
1107 /**
1108 SFrameGenstep::GenerateSimtracePhotons
1109 -----------------------------------------
1110 
1111 Canonical invokation is from SEvt::setFrame_HostsideSimtrace
1112 using the genstep created by SFrameGenstep::MakeCenterExtentGensteps
1113 
1114 **/
1115 
1116 void SFrameGenstep::GenerateSimtracePhotons( std::vector<quad4>& simtrace, const std::vector<quad6>& genstep )
1117 {
1118     LOG(LEVEL)
1119         << " genstep.size " << genstep.size()
1120         << " simtrace.size " << simtrace.size()
1121         ;
1122 
1123     unsigned seed = 0 ;
1124     SRng<float> rng(seed) ;
1125 
1126     bool simtrace_layout = true ;
1127     unsigned count = 0 ;
1128 
1129     for(unsigned i=0 ; i < genstep.size() ; i++)
1130     {
1131         const quad6& gs = genstep[i];
1132         C4U gsid ;   // genstep integer grid coordinate IXYZ and IW photon index up to 255
1133         int gencode           = gs.q0.i.x ;
1134         int gridaxes          = gs.q0.i.y ;
1135         gsid.u                = gs.q0.u.z ;     // formerly gs.q1.u.w
1136         int      num_photons_ = gs.q0.i.w ;     // -ve num_photons_ doesnt use random u0 so creates phi wheels
1137 
1138         qat4 qt(gs) ;  // copy 4x4 transform from last 4 quads of genstep
1139 
1140         unsigned num_photons  = std::abs(num_photons_);
1141         bool expect = gencode == OpticksGenstep_TORCH || gencode == OpticksGenstep_FRAME ;
1142 
1143         LOG_IF(error, !expect)
1144             << "unexpected gencode " << gencode
1145             << " OpticksGenstep_::Name " << OpticksGenstep_::Name(gencode) ;
1146 
1147         assert(expect);
1148         //std::cout << " i " << i << " num_photons " << num_photons << std::endl ;
1149 
1150         double u0, u1 ;
1151         double phi, sinPhi,   cosPhi ;
1152         double sinTheta, cosTheta ;
1153 
1154         for(unsigned j=0 ; j < num_photons ; j++)
1155         {
1156             u0 = num_photons_ < 0 ? double(j)/double(num_photons-1) : rng() ;
1157 
1158             phi = 2.*M_PIf*u0 ;     // azimuthal 0->2pi
1159             ssincos(phi,sinPhi,cosPhi);
1160 
1161             // cosTheta sinTheta are only used for 3D (not 2D planar gensteps)
1162             u1 = rng();
1163             cosTheta = u1 ;
1164             sinTheta = sqrtf(1.0-u1*u1);
1165 
1166             float4 ori ; // copy position from genstep, historically has been origin
1167             ori.x = gs.q1.f.x ;
1168             ori.y = gs.q1.f.y ;
1169             ori.z = gs.q1.f.z ;
1170             ori.w = 1.f ; // <-- dont copy the "float" gsid
1171 
1172             float4 dir ;
1173             SetGridPlaneDirection(dir, gridaxes, cosPhi, sinPhi, cosTheta, sinTheta );
1174 
1175             // transform photon position and direction into the genstep frame
1176             qt.right_multiply_inplace( ori, 1.f );  // position
1177             qt.right_multiply_inplace( dir, 0.f );  // direction
1178 
1179             unsigned char ucj = (j < 255 ? j : 255 ) ;  // photon index local to the genstep
1180             gsid.c4.w = ucj ;     // setting C4U union element to change gsid.u
1181 
1182             unsigned identity = gsid.u ;   // include photon index IW with the genstep coordinate
1183 
1184             quad4& p = simtrace[count] ;
1185             p.zero();
1186 
1187             if(simtrace_layout)  // follow layout of sevent::add_simtrace
1188             {
1189                 p.q0.f = {0.f,0.f,0.f,0.f} ;   // intersect normal and distance
1190                 p.q1.f = {0.f,0.f,0.f,0.f} ;   // intersect position
1191                 p.q2.f = ori ;                 // initial position
1192                 p.q3.f = dir ;                 // initial direction
1193                 p.q3.u.w = identity ;
1194             }
1195             else
1196             {
1197                 p.q0.f = ori ;
1198                 p.q1.f = dir ;
1199                 p.q2.f = {0.f,0.f,0.f,0.f} ;
1200                 p.q3.f = {0.f,0.f,0.f,0.f} ;
1201                 p.q3.u.w = identity ;
1202             }
1203             count += 1 ;
1204         }
1205     }
1206     LOG(LEVEL) << " simtrace.size " << simtrace.size() ;
1207 }
1208 
1209 
1210 
1211 /**
1212 SFrameGenstep::SetGridPlaneDirection
1213 ----------------------------------
1214 
1215 The cosTheta and sinTheta arguments are only used for the 3D gridaxes == XYZ
1216 
1217 **/
1218 
1219 void SFrameGenstep::SetGridPlaneDirection( float4& dir, int gridaxes, double cosPhi, double sinPhi, double cosTheta, double sinTheta )
1220 {
1221     if( gridaxes == YZ )
1222     {
1223         dir.x = 0.f ;
1224         dir.y = float(cosPhi) ;
1225         dir.z = float(sinPhi) ;
1226         dir.w = 0.f    ;
1227     }
1228     else if( gridaxes == XZ )
1229     {
1230         dir.x = float(cosPhi) ;
1231         dir.y = 0.f    ;
1232         dir.z = float(sinPhi) ;
1233         dir.w = 0.f    ;
1234     }
1235     else if( gridaxes == XY )
1236     {
1237         dir.x = float(cosPhi) ;
1238         dir.y = float(sinPhi) ;
1239         dir.z = 0.f     ;
1240         dir.w = 0.f    ;
1241     }
1242     else if( gridaxes == XYZ )    // formerly adhoc used XZ here
1243     {
1244         dir.x = float(sinTheta * cosPhi)  ;
1245         dir.y = float(sinTheta * sinPhi)  ;
1246         dir.z = float(cosTheta) ;
1247         dir.w =  0.f   ;
1248     }
1249     else
1250     {
1251         LOG(fatal) << " invalid gridaxes value " << gridaxes ;
1252         assert(0);
1253     }
1254 }
1255