updates: igor 7, map projections, hemi cuts, longitudinal section

- pearl procedures compile under igor 7.
  some features may not work, e.g. 3D graphics with gizmo.
- add orthographic map projection to angle scans.
- add functions for azimuthal and polar cuts through hemispherical scan.
- add function to load longitudinal section from pshell data file.
This commit is contained in:
2016-10-14 16:56:20 +02:00
parent 600061f684
commit c8a69460bc
4 changed files with 459 additions and 88 deletions

View File

@ -808,9 +808,9 @@ function convert_angles_ttpa2polar(theta, tilt, phi, analyser, polar, azi)
// this is simply a polar-cartesian mapping, independent of the manipulator
// phi=0 is in the polar rotation plane
make /n=(3,na) /d /free w_orig_polar, w_orig_cart, w_rot_cart, w_rot_polar
w_orig_polar[0] = radius
w_orig_polar[1] = analyser[q]
w_orig_polar[2] = 0
w_orig_polar[0][] = radius
w_orig_polar[1][] = analyser[q]
w_orig_polar[2][] = 0
polar2cart_wave(w_orig_polar, w_orig_cart)
// if the angle-dispersive axis was horizontal, we'd need to rotate the detector
//rotate_z_wave(w_orig_cart, 90)
@ -1905,21 +1905,46 @@ function /s display_scanlines(nickname, alpha_lo, alpha_hi, m_theta, m_tilt, m_p
return graphname
end
static constant kProjScaleLinear = 2
/// @page PageProjections Projections
///
/// the functions of the anglescan package support the following map projections.
/// for a description of the different projections, see, for example,
/// https://en.wikipedia.org/wiki/Map_projection
///
/// | Selector | Projection | Function | Properties |
/// | :----: | :----: | :----: | :---- |
/// | kProjDist = 0 | azimuthal equidistant | r = c * theta | radius is proportional to polar angle. |
/// | kProjStereo = 1 | stereographic | r = c * tan theta/2 | circles on sphere map to circles. |
/// | kProjArea = 2 | azimuthal equal-area | r = c * sin theta/2 | preserves area measure. |
/// | kProjGnom = 3 | gnomonic | r = c * tan theta | great circles map to straight lines. |
/// | kProjOrtho = 4 | orthographic | r = c * sin theta | k-space mapping in ARPES and LEED. |
///
/// the projections in this package are defined for 0 <= theta < 90.
///
constant kProjDist = 0
constant kProjStereo = 1
constant kProjArea = 2
constant kProjGnom = 3
constant kProjOrtho = 4
static constant kProjScaleDist = 2
static constant kProjScaleStereo = 2
static constant kProjScaleAzim = 2
static constant kProjScaleArea = 2
// scaled so that radius(gnom) = radius(stereo) for polar = 88
static constant kProjScaleGnomonic = 0.06744519021
static constant kProjScaleGnom = 0.06744519021
static constant kProjScaleOrtho = 2
/// calculate the projected polar angle
///
/// @param polar polar angle in degrees
///
/// @param projection mapping function from polar to cartesian coordinates
/// @arg 0 linear
/// @arg 1 stereographic (default)
/// @arg 2 azimuthal
/// @arg 3 gnomonic (0 <= polar < 90)
/// @param projection mapping function from polar to cartesian coordinates.
/// see @ref PageProjections for details.
/// @arg kProjDist = 0 azimuthal equidistant
/// @arg kProjStereo = 1 stereographic (default)
/// @arg kProjArea = 2 azimuthal equal-area
/// @arg kProjGnom = 3 gnomonic (0 <= polar < 90)
/// @arg kProjOrtho = 4 orthographic
///
/// @return projected radius.
/// the radius is scaled such that grazing emission maps to 2.
@ -1933,17 +1958,20 @@ threadsafe function calc_graph_radius(polar, [projection])
variable radius
switch(projection)
case 1: // stereographic
case kProjStereo: // stereographic
radius = kProjScaleStereo * tan(polar / 2 * pi / 180)
break
case 2: // azimuthal
radius = kProjScaleAzim * cos((180 - polar) / 2 * pi / 180)
case kProjArea: // equal area
radius = kProjScaleArea * sin(polar / 2 * pi / 180)
break
case 3: // gnomonic
radius = polar < 90 ? kProjScaleGnomonic * tan(polar * pi / 180) : inf
case kProjGnom: // gnomonic
radius = polar < 90 ? kProjScaleGnom * tan(polar * pi / 180) : inf
break
default: // linear
radius = kProjScaleLinear * polar / 90
case kProjOrtho: // orthographic
radius = kProjScaleOrtho * sin(polar * pi / 180)
break
default: // equidistant
radius = kProjScaleDist * polar / 90
endswitch
return radius
@ -1955,10 +1983,13 @@ end
///
/// @param x, y projected Cartesian coordinate
///
/// @param projection mapping function from polar to cartesian coordinates
/// @arg 0 linear
/// @arg 1 stereographic (default)
/// @arg 2 azimuthal
/// @param projection mapping function from polar to cartesian coordinates.
/// see @ref PageProjections for details.
/// @arg kProjDist = 0 azimuthal equidistant
/// @arg kProjStereo = 1 stereographic (default)
/// @arg kProjArea = 2 azimuthal equal-area
/// @arg kProjGnom = 3 gnomonic (0 <= polar < 90)
/// @arg kProjOrtho = 4 orthographic
///
/// @returns polar angle in degrees
///
@ -1976,17 +2007,20 @@ threadsafe function calc_graph_polar(x, y, [projection])
radius = sqrt(x^2 + y^2)
switch(projection)
case 1: // stereographic
case kProjStereo: // stereographic
polar = 2 * atan(radius / kProjScaleStereo) * 180 / pi
break
case 2: // azimuthal
polar = 180 - 2 * acos(radius / kProjScaleAzim) * 180 / pi
case kProjArea: // equal area
polar = 2 * asin(radius / kProjScaleArea) * 180 / pi
break
case 3: // gnomonic
polar = atan(radius / kProjScaleGnomonic) * 180 / pi
case kProjGnom: // gnomonic
polar = atan(radius / kProjScaleGnom) * 180 / pi
break
default: // linear
polar = 90 * radius / kProjScaleLinear
case kProjOrtho: // orthographic
polar = asin(radius / kProjScaleOrtho) * 180 / pi
break
default: // equidistant
polar = 90 * radius / kProjScaleDist
endswitch
return polar
@ -1997,10 +2031,13 @@ end
/// @param x, y projected Cartesian coordinate
///
/// @param projection mapping function from polar to cartesian coordinates.
/// projections 0-2 have no effect on the azimuthal coordinate.
/// @arg 0 linear
/// @arg 1 stereographic (default)
/// @arg 2 azimuthal
/// all supported projections are azimuthal, they have no effect on the azimuthal coordinate.
/// see @ref PageProjections for details.
/// @arg kProjDist = 0 azimuthal equidistant
/// @arg kProjStereo = 1 stereographic (default)
/// @arg kProjArea = 2 azimuthal equal-area
/// @arg kProjGnom = 3 gnomonic (0 <= polar < 90)
/// @arg kProjOrtho = 4 orthographic
///
/// @param zeroAngle zeroAngleWhere parameter of polar graphs
/// @arg 0 (default) zero is at the 3 o'clock position
@ -2586,4 +2623,158 @@ function import_tpi_scan(nickname, theta, phi, intensity, [folding, npolar, nogr
endfor
display_hemi_scan(nickname)
end
end
/// extract a polar cut from a hemispherical scan.
///
/// for each polar angle, the function first extracts all azimuthal angles.
/// the intensity is then interpolated between the nearest neighbours of the given azimuth.
///
/// the hemi grid must have been created in the current data folder by the make_hemi_grid function.
/// correct ordering is required.
///
/// @param nickname name of the scan dataset.
/// can be empty if no prefix is used.
/// the dataset must be in the current datafolder.
///
/// @param azim azimuthal angle in degrees
///
/// @return reference of the created wave.
/// the wave has the same name as the intensity wave of the dataset
/// with the suffix "_azi" and the azimuthal angle rounded to integer.
/// it is created in the same datafolder as the original data.
///
function /wave hemi_polar_cut(nickname, azim)
string nickname
variable azim
if (strlen(nickname))
string s_prefix = nickname + "_"
string s_int = s_prefix + "i"
else
s_prefix = ""
s_int = "values"
endif
string s_totals = s_prefix + "tot"
string s_weights = s_prefix + "wt"
string s_polar = s_prefix + "pol"
string s_azim = s_prefix + "az"
string s_index = s_prefix + "index"
string s_theta = s_prefix + "th"
string s_dphi = s_prefix + "dphi"
string s_nphis = s_prefix + "nphis"
string s_cut
sprintf s_cut, "%s_azi%03u", s_int, round(azim)
wave w_polar = $s_polar
wave w_azim = $s_azim
wave w_values = $s_int
wave w_totals = $s_totals
wave w_weights = $s_weights
wave w_index = $s_index
wave w_theta = $s_theta
wave w_dphi = $s_dphi
wave w_nphis = $s_nphis
variable npol = numpnts(w_theta)
variable ipol
variable pol_st = abs(w_theta[1] - w_theta[0])
variable pol
variable pol1, pol2
variable nsel
make /n=(npol) /o $s_cut
wave w_cut = $s_cut
setscale /i x w_theta[0], w_theta[numpnts(w_theta)-1], "deg", w_cut
make /n=1 /free azi_slice
make /n=1 /free values_slice
for (ipol = 0; ipol < npol; ipol += 1)
pol = w_theta[ipol]
pol1 = pol - pol_st / 2
pol2 = pol + pol_st / 2
extract /free /indx w_polar, sel, (pol1 < w_polar) && (w_polar <= pol2)
nsel = numpnts(sel)
if (nsel > 0)
redimension /n=(nsel+2) azi_slice, values_slice
azi_slice[1, nsel] = w_azim[sel[p-1]]
azi_slice[0] = azi_slice[nsel] - 360
azi_slice[nsel+1] = azi_slice[1] + 360
values_slice[1, nsel] = w_values[sel[p-1]]
values_slice[0] = values_slice[nsel]
values_slice[nsel+1] = values_slice[1]
w_cut[ipol] = interp(azim, azi_slice, values_slice)
else
w_cut[ipol] = nan
endif
endfor
return w_cut
end
/// extract an azimuthal cut from a hemispherical scan
///
/// the function extracts all azimuthal angles that are present for the given polar angle.
///
/// the hemi grid must have been created in the current data folder by the make_hemi_grid function.
/// correct ordering is required.
///
/// @param nickname name of the scan dataset.
/// can be empty if no prefix is used.
/// the dataset must be in the current datafolder.
///
/// @param pol polar angle in degrees
///
/// @return reference of the created wave.
/// the wave has the same name as the intensity wave of the dataset
/// with the suffix "_azi" and the azimuthal angle rounded to integer.
/// it is created in the same datafolder as the original data.
///
function /wave hemi_azi_cut(nickname, pol)
string nickname
variable pol
if (strlen(nickname))
string s_prefix = nickname + "_"
string s_int = s_prefix + "i"
else
s_prefix = ""
s_int = "values"
endif
string s_totals = s_prefix + "tot"
string s_weights = s_prefix + "wt"
string s_polar = s_prefix + "pol"
string s_azim = s_prefix + "az"
string s_index = s_prefix + "index"
string s_theta = s_prefix + "th"
string s_dphi = s_prefix + "dphi"
string s_nphis = s_prefix + "nphis"
string s_cut
sprintf s_cut, "%s_pol%03u", s_int, round(pol)
wave w_polar = $s_polar
wave w_azim = $s_azim
wave w_values = $s_int
wave w_totals = $s_totals
wave w_weights = $s_weights
wave w_index = $s_index
wave w_theta = $s_theta
wave w_dphi = $s_dphi
wave w_nphis = $s_nphis
variable pol_st = abs(w_theta[1] - w_theta[0])
variable pol1, pol2
variable nsel
pol1 = pol - pol_st / 2
pol2 = pol + pol_st / 2
extract /free /indx w_polar, sel, (pol1 < w_polar) && (w_polar <= pol2)
nsel = numpnts(sel)
if (nsel > 0)
make /n=(nsel) /o $s_cut
wave w_cut = $s_cut
w_cut = w_values[sel]
setscale /i x w_azim[sel[0]], w_azim[sel[nsel-1]], "deg", w_cut
return w_cut
else
return $""
endif
end

View File

@ -454,7 +454,7 @@ function Au4f(w, x): fitfunc
vc1 = w[ip] / sqrt(pi) * vc2
vc3 = w[ip+1]
vc4 = vc2 * w[ip+2] / 2
v += vc1 * Voigt(vc2 * (x - vc3), vc4)
v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
endfor
return v
@ -491,25 +491,25 @@ function Au4f_2p2(w, x): fitfunc
vc1 = w[3] / sqrt(pi) * vc2
vc3 = w[4]
vc4 = vc2 * w[5] / 2
v += vc1 * Voigt(vc2 * (x - vc3), vc4)
v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
// 5/2 surface
vc1 = w[3] / sqrt(pi) * vc2 * w[9]
vc3 = w[4] + w[10]
vc4 = vc2 * w[5] / 2
v += vc1 * Voigt(vc2 * (x - vc3), vc4)
v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
// 7/2 bulk
vc1 = w[6] / sqrt(pi) * vc2
vc3 = w[7]
vc4 = vc2 * w[8] / 2
v += vc1 * Voigt(vc2 * (x - vc3), vc4)
v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
// 7/2 surface
vc1 = w[6] / sqrt(pi) * vc2 * w[9]
vc3 = w[7] + w[10]
vc4 = vc2 * w[8] / 2
v += vc1 * Voigt(vc2 * (x - vc3), vc4)
v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
return v
@ -578,37 +578,37 @@ function Au4f_2p3(w, x): fitfunc
vc1 = w[3] / sqrt(pi) * vc2
vc3 = w[4]
vc4 = vc2 * w[5] / 2
v += vc1 * Voigt(vc2 * (x - vc3), vc4)
v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
// 5/2 surface
vc1 = w[3] / sqrt(pi) * vc2 * w[9]
vc3 = w[4] + w[10]
vc4 = vc2 * w[5] / 2
v += vc1 * Voigt(vc2 * (x - vc3), vc4)
v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
// 5/2 2nd layer
vc1 = w[3] / sqrt(pi) * vc2 * w[11]
vc3 = w[4] + w[12]
vc4 = vc2 * w[5] / 2
v += vc1 * Voigt(vc2 * (x - vc3), vc4)
v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
// 7/2 bulk
vc1 = w[6] / sqrt(pi) * vc2
vc3 = w[7]
vc4 = vc2 * w[8] / 2
v += vc1 * Voigt(vc2 * (x - vc3), vc4)
v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
// 7/2 surface
vc1 = w[6] / sqrt(pi) * vc2 * w[9]
vc3 = w[7] + w[10]
vc4 = vc2 * w[8] / 2
v += vc1 * Voigt(vc2 * (x - vc3), vc4)
v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
// 7/2 2nd layer
vc1 = w[6] / sqrt(pi) * vc2 * w[11]
vc3 = w[7] + w[12]
vc4 = vc2 * w[8] / 2
v += vc1 * Voigt(vc2 * (x - vc3), vc4)
v += vc1 * VoigtFunc(vc2 * (x - vc3), vc4)
return v

View File

@ -676,6 +676,56 @@ function /s psh5_load_dataset(fileID, scanpath, datasetname, [set_scale])
return dataname
end
/// select the preferred dataset from a list of available datasets.
///
/// @param file_datasets semicolon-separated list of datasets that are available in the file.
///
/// @param pref_datasets semicolon-separated list of preferred datasets.
/// the items of the list are match strings for the Igor StringMatch function.
/// the first matching dataset is loaded from the file.
/// if no match is found, the first file dataset is selected.
///
/// @return name of selected dataset.
///
static function /s select_dataset(file_datasets, pref_datasets)
string file_datasets
string pref_datasets
variable index
variable nds = ItemsInList(file_datasets)
variable ids
string sds = ""
variable np = ItemsInList(pref_datasets)
variable ip
string sp
variable found = 0
if (nds > 0)
for (ip = 0; ip < np; ip += 1)
for (ids = 0; ids < nds; ids += 1)
sds = StringFromList(ids, file_datasets)
index = ItemsInList(sds, "/") - 1
sds = StringFromList(index, sds, "/")
sp = StringFromList(ip, pref_datasets)
if (StringMatch(sds, sp))
found = 1
break
endif
endfor
if (found)
break
endif
endfor
if (!found)
ids = 0
sds = StringFromList(ids, file_datasets)
index = ItemsInList(sds, "/") - 1
sds = StringFromList(index, sds, "/")
endif
endif
return sds
end
/// load a preview dataset from an open PShell HDF5 file.
///
/// if the dataset has a maximum of two dimensions, the function loads it at once.
@ -686,9 +736,6 @@ end
///
/// @param scanpath path to the scan group in the HDF5 file, e.g. "/scan 1".
///
/// @param dataset name of the dataset.
/// the name of the loaded wave is a cleaned up version of the dataset name.
///
/// @param set_scale by default, the function tries to set the wave scaling if the attributes have been loaded.
/// if multiple datasets are loaded from a file,
/// it is more efficient to set the scaling of all loaded datasets at the end by calling ps_scale_datasets().
@ -720,45 +767,11 @@ function /s psh5_load_scan_preview(fileID, scanpath, [set_scale, pref_datasets])
dfref dataDF = saveDF
string datasets = psh5_list_scan_datasets(fileID, scanpath)
variable index
variable nds = ItemsInList(datasets)
variable ids
string sds
variable np = ItemsInList(pref_datasets)
variable ip
string sp
variable found = 0
if (nds > 0)
for (ip = 0; ip < np; ip += 1)
for (ids = 0; ids < nds; ids += 1)
sds = StringFromList(ids, datasets)
index = ItemsInList(sds, "/") - 1
sds = StringFromList(index, sds, "/")
sp = StringFromList(ip, pref_datasets)
if (StringMatch(sds, sp))
found = 1
break
endif
endfor
if (found)
break
endif
endfor
if (!found)
ids = 0
sds = StringFromList(ids, datasets)
index = ItemsInList(sds, "/") - 1
sds = StringFromList(index, sds, "/")
endif
else
return ""
endif
string datasetname = select_dataset(datasets, pref_datasets)
string datasetpath
string datasetname = sds
datasetpath = scanpath + "/" + datasetname
datasetpath = ReplaceString("//", datasetpath, "/")
STRUCT HDF5DataInfo di // Defined in HDF5 Browser.ipf.
InitHDF5DataInfo(di)
variable err = HDF5DatasetInfo(fileID, datasetpath, 0, di)
@ -824,6 +837,153 @@ function /s psh5_load_scan_preview(fileID, scanpath, [set_scale, pref_datasets])
return dataname
end
/// load a longitudinal section of a scan from an open PShell HDF5 file.
///
/// the dataset must have three dimensions.
///
///
/// @param fileID ID of open HDF5 file from psh5_open_file().
///
/// @param scanpath path to the scan group in the HDF5 file, e.g. "/scan 1".
///
/// @param dim reserved, must be 0.
///
/// @param set_scale by default, the function tries to set the wave scaling if the attributes have been loaded.
/// if multiple datasets are loaded from a file,
/// it is more efficient to set the scaling of all loaded datasets at the end by calling ps_scale_datasets().
/// @arg 1 (default) set the wave scaling.
/// @arg 0 do not set the wave scaling.
///
/// @param pref_datasets semicolon-separated list of preferred datasets.
/// the items of the list are match strings for the Igor StringMatch function.
/// the first matching dataset is loaded from the file.
/// if no match is found, the first dataset listed in the file is loaded.
/// if empty, a hard-coded default preference list is used.
///
/// @return name of loaded wave if successful. empty string otherwise.
///
/// @warning EXPERIMENTAL: this function is under development.
///
function /s psh5_load_scan_section(fileID, scanpath, dim, [set_scale, pref_datasets])
variable fileID
string scanpath
variable dim
variable set_scale
string pref_datasets
// select first dimension (future argument)
// 0 = first dimension is x axis (energy of scienta image)
dim = 0
if (ParamIsDefault(set_scale))
set_scale = 1
endif
if (ParamIsDefault(pref_datasets) || (strlen(pref_datasets) == 0))
pref_datasets = kPreviewDatasets
endif
dfref saveDF = GetDataFolderDFR()
dfref dataDF = saveDF
string datasets = psh5_list_scan_datasets(fileID, scanpath)
string datasetname = select_dataset(datasets, pref_datasets)
string datasetpath
datasetpath = scanpath + "/" + datasetname
datasetpath = ReplaceString("//", datasetpath, "/")
string dataname = StringFromList(ItemsInList(datasetpath, "/") - 1, datasetpath, "/")
string destname = dataname[0,29] + num2str(dim)
STRUCT HDF5DataInfo di // Defined in HDF5 Browser.ipf.
InitHDF5DataInfo(di)
variable err = HDF5DatasetInfo(fileID, datasetpath, 0, di)
if (err != 0)
print "error accessing detector/data"
return ""
elseif (di.ndims != 3)
print "error: rank of dataset != 3"
return ""
endif
variable idx, idy, idz, idt
variable transpose = WhichListItem(dataname, kTransposedDatasets) >= 0
if (transpose)
idx = 1
idy = 0
else
idx = 0
idy = 1
endif
idz = 2
idt = 3
variable nx, ny, nz
nx = di.dims[idx]
ny = di.dims[idy]
nz = di.dims[idz]
HDF5MakeHyperslabWave(GetDataFolder(1) + "slab", max(di.ndims, 4))
wave slab
slab[][%Start] = 0
slab[][%Stride] = 1
slab[][%Count] = 1
slab[][%Block] = 1
if (dim == 0)
slab[idy][%Start] = floor(ny / 2)
slab[idx][%Block] = nx
make /n=(nx,nz) /o $destname
else
slab[idx][%Start] = floor(nx / 2)
slab[idy][%Block] = ny
make /n=(ny,nz) /o $destname
endif
slab[idz][%Block] = nz
wave data = $destname
data = 0
HDF5LoadData /O /Q /Z /SLAB=slab /N=slabdata fileID, datasetpath
if (!v_flag)
wave slabdata
if (transpose)
data += slabdata[0][p][q][0]
else
data += slabdata[p][0][q][0]
endif
endif
killwaves /z slab, slabdata
if (set_scale)
make /n=(1,1,1) /free dummy
ps_set_dimlabels2(dummy, dataname)
setdimlabel 0, -1, $GetDimLabel(dummy, dim, -1), data
setdimlabel 1, -1, $kScanDimLabel, data
setdatafolder dataDF
string positioners
string positioner
string positionerpath
positioners = psh5_load_scan_meta(fileID, scanpath)
wave /t /z ScanWritables
if (waveexists(ScanWritables) && (numpnts(ScanWritables) >= 1))
positioner = ScanWritables[0]
if (strlen(positioner) > 0)
positionerpath = scanpath + "/" + positioner
positionerpath = ReplaceString("//", positionerpath, "/")
HDF5LoadData /O /Q /Z fileID, positionerpath
endif
endif
setdatafolder dataDF
newdatafolder /o/s attr
killwaves /a/z
psh5_load_scan_attrs(fileID, scanpath, attr_sets=2)
setdatafolder dataDF
ps_scale_dataset(data)
endif
return destname
end
/// load metadata of a PShell dataset.
///
/// "metadata" are the HDF5 attributes attached to the scan dataset.
@ -1136,9 +1296,27 @@ end
function ps_set_dimlabels(data)
wave data
string name = NameOfWave(data)
variable dummy
ps_set_dimlabels2(data, NameOfWave(data))
end
/// set dimension labels according to the axis type
///
/// same as ps_set_dimlabels() except that the dimension labels are set
/// according to a separate name argument instead of the wave name.
///
/// @param data data wave as loaded from PShell file.
///
/// @param name original name of the dataset in the PShell file.
///
/// @return @arg 0 all labels set successfully.
/// @arg 1 unidentified data source.
/// @arg 2 wave does not contain data.
///
function ps_set_dimlabels2(data, name)
wave data
string name
variable dummy
try
// intrinsic dimensions
strswitch(name)