Up for discussion: Cut_throughs with curves, spectrum calculations, static streamlines vol centering#147
Conversation
… vector outputs (I dont remember anymore why I did that.)
…ll at least plot with out tex.
…ent usage and began work on curve support for cut-through
…ill, tested against a regular cut-through, seems to work. Tightened the epsilon in regular cut-through.
…rt for static_field_tracer, with a) automatic detection of vol in variable name and b) a keyword parameter to force centering. Added a spherical inner boundary trace stop parameter to static_field_tracer.
…ra v-cell coordinate sampling, commented out some debug lines. Added possibility of selecting spectra weights.
…e.g. flux function values at null points to be plotted with flux function contours
…its (1/m^3 is hard to automagicize)
…, get_scaled_var comments added
|
so def get_scaled_var() replicates what was somewhat available already in Line 187 in bf074f0 Could you then also rework the current plotting tools to use the new method so we don't have overlapping approaches and fragmentation of code? |
|
So, the get_scaled_var is a bit suboptimal approach as it copies the whole array. Should we implement slicing operators for variable_info, so that for some variable_info var: var[index] would return var.data[index]*vscale? |
…mpling 3) variableinfo specialist units reworking
…aleunits behaviour, and pointed occurrences of that to the new scaler using the new metadata getter.
|
So the above |
…izes, fixed variable scaler unknown keyword exception
|
So this PR is getting a bit crowded, since I have built whatever I have used in eVlasiator survey on this branch - up for merge? |
markusbattarbee
left a comment
markusbattarbee
left a comment
There was a problem hiding this comment.
Not a lot of changes but e.g. some suggestions about default behaviour
| print("Error in rotation: testvector ",count,testvect," not a unit vector") | ||
| if abs(1.0-np.amax(testvect))>1.e-3: | ||
| print("Error in rotation: testvector ",count,testvect," largest component is not unity") | ||
| else: |
There was a problem hiding this comment.
So what does this resampler actually do that is different from the old version?
There was a problem hiding this comment.
It handles the re-meshing the original v-mesh to a cubic grid aligned with the output bins, so there is no alias effect from the histogramming. I have to say though, that that is not implemented for each case, or slicethickness. Blimey, that needs still some work.
There was a problem hiding this comment.
Ok, looking forward to seeing that in action in a complete form!
There was a problem hiding this comment.
Hmm. If the resampleReducer does not need the VX, VY, Voutofslice values etc, shouldn't it be called under a separate if at around line 330? And then it can return before doing all these vector calculations and rotations.
Also, because the resampleReducer doesn't support flux weighting (not that anyone has used it here at all), there should probably be a check to see if that combo was requested? or if flux weighting is requested, re-route to the old histogramming.
|
A bunch of updates, incl. more controls for vdf plotting, and better (yet incomplete) support for the vdf resampler. |
|
So the resampler is now a bit more finished, but I noticed also I haven't made the switch of defaults yet. |
markusbattarbee
left a comment
markusbattarbee
left a comment
There was a problem hiding this comment.
A bit of head-scratching about flux weighting, among other stuff.
For plot_vdf we could even consider dropping flux weighting as an option altogether as it doesn't really make much sense, but on the other hand might as well have it tag along. But I think (?) that there are issues with the integrations inside spectra.py
| Ekin[Ekin > max(EkinBinEdges)] = max(EkinBinEdges) | ||
| # normalization | ||
| if(weight == 'flux'): | ||
| fw = f*np.sqrt(V2) # use particle flux as weighting |
There was a problem hiding this comment.
For differential flux density, we should have:
fw = f * (dV**3) * np.sqrt(V2) / (Ekin * 4*np.pi) for particles per steradian per eV per square metre per second;
units = "1/(s m^2 sr eV)"
latexunits = '$\mathrm{s}^{-1}\,\mathrm{m}^{-2}\,\mathrm{sr}^{-1}\,\mathrm{eV}^{-1}$'
latex='$\mathcal{F}(\vec{r})$'
| latexunits = '$\mathrm{s}\,\mathrm{m}^{-2}\,(4\pi\,\mathrm{eV})^{-1}$' | ||
| latex='$f(\vec{r})v\,\DeltaE\,sr^-1$' | ||
| else: | ||
| fw = f |
There was a problem hiding this comment.
and for number density we should have
fw = f * (dV**3) / Ekin
units = "1/(m^3 eV)"
latexunits = '$\mathrm{m}^{-3}\,\mathrm{eV}^{-1}$'
latex='$f(\vec{r})$'
| weight = 'particles' | ||
| dV = np.prod(vlsvReader.get_velocity_mesh_dv(population)) | ||
| # compute histogram | ||
| (nhist,edges) = np.histogram(Ekin,bins=EkinBinEdges,weights=fw*dV,normed=0) |
There was a problem hiding this comment.
And then we should get rid of the dV here and on lines 85, 88, methinks. Thoughts?
Similar corrections needed to the other two routines as well.
There was a problem hiding this comment.
So the dV(**3) goes in on line 83 (and ftotal wasn't even used from 85, got rid of that), so that was accounted there. But yes, maybe better move that to a more explicit position.
... but blergh. The energy bin width is accounted for, but yeah, that Ekin is missing all the way from the old create_time_energy_spectrogram.py. It does make sense in that that then integration over the E bins (taking the edge values) gives the total, again. With log bins, \Delta E even goes \propto E, so it's... almost the same, cough.
…flux weighting for alongaxis-velocity; feel free to reimplement if those are needed, _in the specific manner required_.
…ral warning about using flux weighting.
…ersion to eV in spectra.py instead of a mess-able q kword, removed unused frame kword (anyway better for the user to do whatever conversions explicitly?)
2 participants
Cut_through for arbitrary curves, which precipitated some fixes to cut_through indexing etc.,
and proper volume-centering options for static_streamline. Lifted vdf energy spectrum calculations from old scripts to pyCalculations for ease of use, need to at least finalize comments/help texts.