pearl-anglescan-process.ipf

Go to the documentation of this file.

#pragma rtGlobals=3// Use modern global access method and strict wave access.
#pragma version = 1.8
#pragma IgorVersion = 6.2
#pragma ModuleName = PearlAnglescanProcess
#include "pearl-vector-operations"
#include "pearl-polar-coordinates"
#include <New Polar Graphs>

// copyright (c) 2013-17 Paul Scherrer Institut
//
// Licensed under the Apache License, Version 2.0 (the "License");
//   you may not use this file except in compliance with the License.
//   You may obtain a copy of the License at
//   http:///www.apache.org/licenses/LICENSE-2.0
//
// Please acknowledge the use of this code.




variable strip_delete_frames(wave strip, variable qlo, variable qhi, wave theta, wave tilt, wave phi){
    wave strip// 2D data, X-axis = analyser angle, Y-axis = arbitrary manipulator scan
    variable qlo
    variable qhi
    wave theta
    wave tilt
    wave phi
    
    if (qlo > qhi)
        return -1
    endif
    
    // source indices
    variable snx = dimsize(strip, 0)
    variable sny = dimsize(strip, 1)
    variable sq1lo = 0
    variable sq1hi = max(qlo-1, 0)
    variable sq2lo = min(qhi+1, sny - 1)
    variable sq2hi = dimsize(strip, 1) - 1
    
    // dest indices
    variable dnx = snx
    variable dny = sny - (sq2lo - sq1hi + 1)
    variable dq1lo = 0
    variable dq1hi = sq1hi
    variable dq2lo = dq1hi + 1
    variable dq2hi = dny - 1
    variable q1ofs = sq1lo - dq1lo
    variable q2ofs = sq2lo - dq2lo
    
    duplicate /free strip, strip_copy
    redimension /n=(dnx,dny) strip
    strip[][dq1lo,dq1hi] = strip_copy[p][q + q1ofs]
    strip[][dq2lo,dq2hi] = strip_copy[p][q + q2ofs]
    
    duplicate /free theta, theta_copy
    redimension /n=(dny) theta
    theta[dq1lo,dq1hi] = theta_copy[p + q1ofs]
    theta[dq2lo,dq2hi] = theta_copy[p + q2ofs]
    
    duplicate /free tilt, tilt_copy
    redimension /n=(dny) tilt
    tilt[dq1lo,dq1hi] = tilt_copy[p + q1ofs]
    tilt[dq2lo,dq2hi] = tilt_copy[p + q2ofs]
    
    duplicate /free phi, phi_copy
    redimension /n=(dny) phi
    phi[dq1lo,dq1hi] = phi_copy[p + q1ofs]
    phi[dq2lo,dq2hi] = phi_copy[p + q2ofs]
    
    return 0
};

variable normalize_strip_x(wave strip, variable smooth_method = defaultValue, variable smooth_factor = defaultValue, variable check = defaultValue){
    wave strip
    variable smooth_method
    variable smooth_factor
    variable check
    
    if (ParamIsDefault(smooth_method))
        smooth_method = 1
    endif
    if (ParamIsDefault(smooth_factor))
        switch(smooth_method)
            case 4:
                smooth_factor = 0.5
                break
            default:
                smooth_factor = 2
        endswitch
    endif
    if (ParamIsDefault(check))
        check = 0
    endif
    
    // average over all scan positions
    wave dist = ad_profile_x(strip, -inf, inf, "")
    variable div = mean(dist)
    dist /= div

    if (check)
        duplicate /o dist, check_dist
    endif
    
    // smooth distribution function
    switch(smooth_method)
        case 1:
            Smooth /B /E=3 smooth_factor, dist
            break
        case 2:
            Smooth /E=3 smooth_factor, dist
            break
        case 3:
            make /n=1 /d /free fit_params
            fit_scienta_ang_transm(dist, fit_params)
            dist = scienta_ang_transm(fit_params, x)
            break
        case 4:
            loess /smth=(smooth_factor) srcWave=dist
            break
    endswitch
    
    if (check)
        duplicate /o dist, check_smoo
    endif

    // divide
    if (check != 2)
        strip /= dist[p]
    endif
};

variable normalize_strip_phi(wave strip, wave theta, wave phi, variable theta_offset = defaultValue, variable theta_range = defaultValue, variable check = defaultValue){
    wave strip
    wave theta
    wave phi
    variable theta_offset
    variable theta_range
    variable check
    
    if (ParamIsDefault(check))
        check = 0
    endif
    if (ParamIsDefault(theta_offset))
        theta_offset = 0
    endif
    if (ParamIsDefault(theta_range))
        theta_offset = 10
    endif

    // average over analyser angles
    wave dist = ad_profile_y(strip, -inf, inf, "")
    
    // smooth distribution function
    duplicate /free dist, dist_smoo
    duplicate /free theta, theta_int
    theta_int = theta - theta_offset
    duplicate /free phi, phi_int
    setscale /p x phi_int[0], phi_int[1] - phi_int[0], waveunits(phi, -1), dist, dist_smoo

    extract /free /indx dist, red_idx, theta_int < theta_range
    duplicate /free red_idx, red_dist, red_phi
    red_dist = dist[red_idx]
    red_phi = phi_int[red_idx]
    
    variable wavg = mean(red_dist)
    make /n=4 /d /free coef
    coef[0] = {wavg, wavg/100, pi/180, 0}
    CurveFit /q /h="0010" /g /w=2 sin, kwcWave=coef, red_dist /x=red_phi
    dist_smoo = coef[0] + coef[1] * sin(coef[2] * phi_int[p] + coef[3])
    
    // divide
    if (check != 2)
        strip = strip / dist_smoo[q] * coef[0]
    endif
    
    // check
    if (check)
        duplicate /o dist, check_dist
        duplicate /o dist_smoo, check_smoo
        setscale /p x dimoffset(dist,0), dimdelta(dist,0), waveunits(dist,0), check_dist, check_smoo
    endif
};

variable normalize_strip_theta(wave strip, wave theta, variable theta_offset = defaultValue, variable smooth_method = defaultValue, variable smooth_factor = defaultValue, variable check = defaultValue){
    wave strip
    wave theta
    variable theta_offset
    variable smooth_method
    variable smooth_factor
    variable check
    
    if (ParamIsDefault(check))
        check = 0
    endif
    if (ParamIsDefault(theta_offset))
        theta_offset = 0
    endif
    if (ParamIsDefault(smooth_method))
        smooth_method = 4
    endif
    if (ParamIsDefault(smooth_factor))
        smooth_factor = 0.5
    endif

    // average over analyser angles
    wave dist = ad_profile_y(strip, -inf, inf, "")
    
    // smooth distribution function
    duplicate /free dist, dist_smoo
    duplicate /free theta, theta_int
    theta_int = theta - theta_offset
    setscale /p x theta_int[0], theta_int[1] - theta_int[0], waveunits(theta,-1), dist, dist_smoo
    variable nx = dimsize(strip, 0)
    variable ix
    
    switch(smooth_method)
        case 1:
            Smooth /B /E=3 smooth_factor, dist_smoo
            break
        case 2:
            Smooth /E=3 smooth_factor, dist_smoo
            break
        case 4:
            loess /dest=dist_smoo /smth=(smooth_factor) srcWave=dist, factors={theta_int}
            break
        case 3:
            for (ix = 0; ix < nx; ix += 1)
                dist = strip[ix][p]
                if (smooth_factor > 1)
                    CurveFit /nthr=0 /q /w=2 poly smooth_factor+1, dist /x=theta_int /d=dist_smoo
                else
                    CurveFit /nthr=0 /q /w=2 line, dist /x=theta_int /d=dist_smoo
                endif
                strip[ix,ix][] /= dist_smoo[q]
            endfor
            dist_smoo = 1
            break
    endswitch
    
    // divide
    if (check != 2)
        strip /= dist_smoo[q]
    endif
    
    // check
    if (check)
        duplicate /o dist, check_dist
        duplicate /o dist_smoo, check_smoo
        setscale /p x dimoffset(dist,0), dimdelta(dist,0), waveunits(dist,0), check_dist, check_smoo
    endif
};

variable normalize_strip_2d(wave strip, wave theta, variable theta_offset = defaultValue, variable smooth_method = defaultValue, variable smooth_factor = defaultValue, variable check = defaultValue){
    wave strip
    wave theta
    variable theta_offset
    variable smooth_method
    variable smooth_factor
    variable check
    
    if (ParamIsDefault(check))
        check = 0
    endif
    if (ParamIsDefault(theta_offset))
        theta_offset = 0
    endif
    if (ParamIsDefault(smooth_method))
        smooth_method = 4
    endif
    if (ParamIsDefault(smooth_factor))
        smooth_factor = 0.5
    endif

    variable nx = dimsize(strip, 0)
    variable ny = dimsize(strip, 1)
    
    duplicate /free strip, dist, alpha_int, theta_int
    theta_int = theta[q] - theta_offset
    alpha_int = dimoffset(strip, 0) + p * dimdelta(strip, 0)
    redimension /n=(nx * ny) dist, alpha_int, theta_int
    
    switch(smooth_method)
        case 4:
            loess /dest=dist_smoo /smth=(smooth_factor) srcWave=dist, factors={alpha_int, theta_int}
            redimension /n=(nx, ny) dist_smoo
            break
        default:
            Abort "undefined smooth method"
            break
    endswitch
    
    // divide
    if (check != 2)
        strip /= dist_smoo
    endif
    
    // check
    if (check)
        //duplicate /o dist, check_dist
        duplicate /o dist_smoo, check_smoo
    endif
};

variable crop_strip(wave strip, variable xlo, variable xhi){
    wave strip
    variable xlo
    variable xhi
    
    variable plo = round((xlo - dimoffset(strip, 0)) / dimdelta(strip, 0))
    variable phi = round((xhi - dimoffset(strip, 0)) / dimdelta(strip, 0))
    xlo = plo * dimdelta(strip, 0) + dimoffset(strip, 0)
    xhi = phi * dimdelta(strip, 0) + dimoffset(strip, 0)
    variable nx = phi - plo + 1
    variable ny = dimsize(strip, 1)
    
    duplicate /free strip, strip_copy
    redimension /n=(nx,ny) strip
    strip = strip_copy[p + plo][q]
    setscale /i x xlo, xhi, waveunits(strip, 0), strip
};

variable pizza_service(wave data, string nickname, variable theta_offset, variable tilt_offset, variable phi_offset, variable npolar = defaultValue, variable nograph = defaultValue, variable folding = defaultValue, variable xpdplot = defaultValue){
    wave data
    string nickname
    variable theta_offset
    variable tilt_offset
    variable phi_offset
    variable npolar
    variable nograph
    variable folding
    variable xpdplot
        
    if (ParamIsDefault(npolar))
        npolar = 91
    endif
    if (ParamIsDefault(nograph))
        nograph = 0
    endif
    if (ParamIsDefault(folding))
        folding = 1
    endif
    if (ParamIsDefault(xpdplot))
        xpdplot = 0
    endif
    
    // sort out data folder structure
    dfref saveDF = GetDataFolderDFR()
    dfref dataDF = GetWavesDataFolderDFR(data)
    setdatafolder dataDF
    if (DataFolderExists(":attr"))
        setdatafolder :attr
    endif
    dfref attrDF = GetDataFolderDFR()
    
    wave /sdfr=attrDF ManipulatorTheta
    wave /sdfr=attrDF ManipulatorTilt
    wave /sdfr=attrDF ManipulatorPhi

    if ((dimsize(ManipulatorTheta, 0) != dimsize(data, 1)) || (dimsize(ManipulatorTilt, 0) != dimsize(data, 1)) || (dimsize(ManipulatorPhi, 0) != dimsize(data, 1)))
        Abort "Warning: The dimension size of the manipulator waves does not match the Y dimension of the data wave!\rIf you restructured the data wave, please use pizza_service_2 with properly scaled manipulator waves."
    endif

    duplicate /free ManipulatorTheta, m_theta
    duplicate /free ManipulatorTilt, m_tilt
    duplicate /free ManipulatorPhi, m_phi

    m_theta -= theta_offset
    m_tilt -= tilt_offset
    m_phi -= phi_offset
    
    pizza_service_2(data, nickname, m_theta, m_tilt, m_phi, npolar=npolar, nograph=nograph, folding=folding, xpdplot=xpdplot)
    
    setdatafolder saveDF
};

variable pizza_service_2(wave data, string nickname, wave m_theta, wave m_tilt, wave m_phi, variable npolar = defaultValue, variable nograph = defaultValue, variable folding = defaultValue, variable xpdplot = defaultValue){
    wave data
    string nickname
    wave m_theta
    wave m_tilt
    wave m_phi
    variable npolar
    variable nograph
    variable folding
    variable xpdplot
        
    if (ParamIsDefault(npolar))
        npolar = 91
    endif
    if (ParamIsDefault(nograph))
        nograph = 0
    endif
    if (ParamIsDefault(folding))
        folding = 1
    endif
    if (ParamIsDefault(xpdplot))
        xpdplot = 0
    endif
    
    if ((dimsize(m_theta, 0) != dimsize(data, 1)) || (dimsize(m_tilt, 0) != dimsize(data, 1)) || (dimsize(m_phi, 0) != dimsize(data, 1)))
        Abort "Warning: The dimension size of the manipulator waves does not match the Y dimension of the data wave!"
    endif
    
    string graphname = "graph_" + nickname
    string outprefix = nickname
    
    // sort out data folder structure
    dfref saveDF = GetDataFolderDFR()
    dfref dataDF = GetWavesDataFolderDFR(data)
    setdatafolder dataDF

    if (xpdplot)
        setdatafolder root:
        outprefix = nickname
    else
        setdatafolder dataDF
        newdatafolder /s/o $nickname
        outprefix = ""
    endif
    dfref destDF = GetDataFolderDFR()
    
    // performance monitoring
    variable timerRefNum
    variable /g pol_perf_secs
    timerRefNum = startMSTimer
    
    duplicate /free m_tilt, corr_tilt
    duplicate /free m_phi, corr_phi
    corr_tilt = -m_tilt// checked 140702
    corr_phi = m_phi// checked 140702
    
    make /n=1/d/free d_polar, d_azi
    
    convert_angles_ttpd2polar(m_theta, corr_tilt, corr_phi, data, d_polar, d_azi)
    d_azi += 180// changed 151030 (v1.6)
    make_hemi_grid(npolar, outprefix, xpdplot=xpdplot)
    variable ifold
    for (ifold = 0; ifold < folding; ifold += 1)
        d_azi = d_azi >= 360 ? d_azi - 360 : d_azi
        hemi_add_anglescan(outprefix, data, d_polar, d_azi)
        d_azi += 360 / folding
    endfor
    
    // normalize folding
    if (strlen(outprefix))
        string s_prefix = outprefix + "_"
        string s_int = s_prefix + "i"
    else
        s_prefix = ""
        s_int = "values"
    endif
    if (folding > 1)
        wave values = $s_int
        values /= folding
    endif
    
    if (!nograph)
        display_hemi_scan(outprefix, graphname = graphname)
    endif
    
    if (timerRefNum >= 0)
        pol_perf_secs = stopMSTimer(timerRefNum) / 1e6
    endif

    setdatafolder saveDF
};

variable show_analyser_line(variable theta, variable tilt, variable phi, variable theta_offset, variable tilt_offset, variable phi_offset, variable npolar = defaultValue, variable nograph = defaultValue, variable xpdplot = defaultValue){
    variable theta
    variable tilt
    variable phi
    variable theta_offset
    variable tilt_offset
    variable phi_offset
    variable npolar
    variable nograph
    variable xpdplot
        
    string nickname = "analyser"

    if (ParamIsDefault(npolar))
        npolar = 91
    endif
    if (ParamIsDefault(nograph))
        nograph = 0
    endif
    if (ParamIsDefault(xpdplot))
        xpdplot = 0
    endif
    string graphname = "graph_" + nickname
    string outprefix = nickname
    
    // sort out data folder structure
    dfref saveDF = GetDataFolderDFR()
    dfref dataDF = saveDF
    if (xpdplot)
        setdatafolder root:
        outprefix = nickname
    else
        setdatafolder dataDF
        newdatafolder /s/o $nickname
        outprefix = ""
    endif
    dfref destDF = GetDataFolderDFR()
    
    make /n=1 /free m_theta
    make /n=1 /free m_tilt
    make /n=1 /free m_phi
    m_theta = theta - theta_offset
    m_tilt = tilt - tilt_offset
    m_tilt *= -1// checked 140702
    m_phi = phi - phi_offset
    //m_phi *= -1 // checked 140702
    
    make /n=60 /free data
    setscale /i x -30, 30, data
    data = x
    make /n=1/d/free d_polar, d_azi
    
    convert_angles_ttpa2polar(m_theta, m_tilt, m_phi, data, d_polar, d_azi)
    d_azi += 180// changed 151030 (v1.6)
    d_azi = d_azi >= 360 ? d_azi - 360 : d_azi
    make_hemi_grid(npolar, outprefix, xpdplot=xpdplot)
    hemi_add_anglescan(outprefix, data, d_polar, d_azi)
    
    if (!nograph)
        display_hemi_scan(outprefix, graphname = graphname)
    endif
    
    setdatafolder saveDF
};


variable convert_angles_ttpd2polar(wave theta, wave tilt, wave phi, wave data, wave polar, wave azi){
    wave theta, tilt, phi// see convert_angles_ttpa2polar
    wave data// in, 1D or 2D
        // X-scale must be set to analyser angle scale
    wave polar, azi// see convert_angles_ttpa2polar
    
    make /n=(dimsize(data, 0)) /d /free ana
    setscale /p x dimoffset(data, 0), dimdelta(data, 0), waveunits(data, 0), ana
    ana = x
    convert_angles_ttpa2polar(theta, tilt, phi, ana, polar, azi)
};

variable convert_angles_ttpa2polar(wave theta, wave tilt, wave phi, wave analyser, wave polar, wave azi){
    wave theta
    wave tilt
    wave phi
    wave analyser
    wave polar, azi
        
    variable nn = numpnts(theta)
    variable na = numpnts(analyser)
    redimension /n=(na, nn) polar, azi

    variable radius = 1// don't need to specify - everything is scalable
    
    // step 1: calculate cartesian detection vectors at normal emission
    // 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
    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)

    variable ii
    for (ii = 0; ii < nn; ii += 1)
        // step 2: rotate the detection vectors according to the manipulator angles
        // the order of rotations is important because we rotate about fixed axes
        // y-axis = tilt rotation axis
        // x-axis = polar rotation axis
        // z-axis = normal emission = azimuthal rotation axis
        w_rot_cart = w_orig_cart
        rotate_y_wave(w_rot_cart, -tilt[ii])
        rotate_x_wave(w_rot_cart, -theta[ii])
        rotate_z_wave(w_rot_cart, -phi[ii] - 90)
        // map the vectors back to the sample coordinate system
        cart2polar_wave(w_rot_cart, w_rot_polar)
        // copy to output
        polar[][ii] = w_rot_polar[1][p]
        azi[][ii] = w_rot_polar[2][p]
    endfor
};

static variable line_average(wave source, wave dest){
    // is this function used?
    wave source
    wave dest
    
    variable ii
    variable nn = dimsize(source, 1)
    make /n=(dimsize(source, 0))/d/free line
    for (ii = 0; ii < nn; ii += 1)
        line = source[p][ii]
        wavestats /q line
        dest[][ii] = line[p] / v_max
    endfor
};

static variable calc_nth(variable Theta_st, variable Theta_in, variable th, variable Phi_ran, variable Phi_ref, string Holomode){
    Variable Theta_st, Theta_in, th, Phi_ran, Phi_ref
       String Holomode
    Variable The_step
    Variable deg2rad=0.01745329
        
    if ( cmpstr(Holomode, "Stereographic") == 0)
       The_step =trunc( Phi_ran*sin(th*deg2rad)*Phi_ref/Theta_st )
       if(th==90)
        The_step =trunc( Phi_ran*sin(th*pi/180)*Phi_ref/Theta_st )
       endif
    else
        if (cmpstr(Holomode, "Parallel") == 0)
            The_step=trunc( Phi_ran*tan(th*deg2rad)*Phi_ref/Theta_st )
        else
            if ( cmpstr(Holomode, "h") == 0)
                    The_step=trunc( th/Theta_in*Phi_ran/Theta_st )
            else
                //altro
            endif
        endif
    endif
    
    return(The_step)
};

static variable calc_phi_step(variable Theta_in, variable th, variable Theta_st, variable Phi_ran, variable Phi_ref, string Holomode){
    Variable Theta_in, th, Theta_st, Phi_ran, Phi_ref
    String Holomode
    
    Variable Phi_st
    Variable deg2rad=0.01745329

    if ( cmpstr(Holomode, "Stereographic") == 0 )
        if ((th < 0.5) || (trunc(Phi_ran*sin(th*deg2rad)*Phi_ref/Theta_st) == 0))
            Phi_st=0.0
        else
            Phi_st=Phi_ran/(trunc(Phi_ran*sin(th*deg2rad)*Phi_ref/Theta_st))
        endif
        if(th==90)
            Phi_st=2.0
        endif
    endif

    if ( cmpstr(Holomode, "Parallel") == 0 )
        if((th < 0.5) || (trunc(Phi_ran*tan(th*deg2rad)*Phi_ref/Theta_st) == 0))
            Phi_st=0.0
        else
            Phi_st=Phi_ran/(trunc(Phi_ran*tan(th*deg2rad)*Phi_ref/Theta_st))
        endif
    endif

    if ( cmpstr(Holomode, "h") == 0 )
        if((th < 0.5) || (trunc(Phi_ran*sin(th*deg2rad)*Phi_ref/Theta_st) == 0))
            Phi_st=0.0
        else
            Phi_st=Phi_ran/trunc(th/Theta_in*Phi_ran/Theta_st)
        endif
    endif
    
    if (Phi_st==0)
        Phi_st=360
    endif
    
    return(Phi_st)
};

static variable Calc_The_step(variable th, variable Theta_st, string Holomode){
    String Holomode
    Variable th, Theta_st
    
    Variable deg2rad=0.01745329, dt_loc,The_step

    if ( (cmpstr(Holomode, "Stereographic")) ==0 )
        The_step=Theta_st
    endif
    
    if ( (cmpstr(Holomode, "h")) ==0 )
        The_step=Theta_st
    endif
 
    if ( cmpstr(Holomode, "Parallel") == 0 )
        if(th < 89.5)
            dt_loc = Theta_st/cos(th*deg2rad)
            if(dt_loc > 10)
                dt_loc=10
            endif
            The_step=dt_loc
        else
            The_step=10
        endif
    endif
    return(The_step)
};

static variable CalcN_Theta(string HoloMode, variable Theta_in, variable Theta_ran, variable Theta_st){
    String HoloMode
    Variable Theta_in,Theta_ran,Theta_st
    Variable n_theta, aux, aux1,ii

    aux = Theta_in
    aux1= Theta_in - Theta_ran
    ii = 0
    do
        aux = aux - Calc_The_step(aux, Theta_st, HoloMode)
        if(aux<=Theta_in-Theta_ran)
            aux=Theta_in-Theta_ran
        endif
        ii = ii+1
    while((aux>aux1)%&(Theta_in-aux<=Theta_ran))//
    n_theta=ii+1
    Return(n_theta)
};

variable make_hemi_grid(variable npol, string nickname, variable xpdplot = defaultValue){
    variable npol
    string nickname
    variable xpdplot
    
    if (ParamIsDefault(xpdplot))
        xpdplot = 0
    endif
    
    string HoloMode = "h"
    variable Theta_in = 90
    variable Theta_ran = 90
    variable Theta_st = 90 / (npol - 1)
    variable Phi_ran = 360
    variable Phi_ref = 1
    variable Phi_in = 0
    
    variable n_theta = CalcN_Theta(HoloMode, Theta_in, Theta_ran, Theta_st)

    // wave names
    if (strlen(nickname))
        string s_prefix = nickname + "_"
        string s_int = s_prefix + "i"// Intensity wave (counts/sec)
    else
        s_prefix = ""
        s_int = "values"// "i" is not a valid wave name
    endif
    string s_polar = s_prefix + "pol"// thetas for each int-point of the holo
    string s_azim = s_prefix + "az"// phis for each int-point of the holo

    string s_index = s_prefix + "index"// starting index for each theta
    string s_theta = s_prefix + "th"// theta values
    string s_dphi = s_prefix + "dphi"// delta phis at each theta
    string s_nphis = s_prefix + "nphis"// number of phis at each theta

    string s_HoloData = s_prefix + "data"// All holo exp.- parameter information
    string s_HoloInfo = s_prefix + "info"
    
    // the following two waves are used by the pearl-anglescan procedures but not by XPDplot
    string s_tot = s_prefix + "tot"// accumulated counts at each point
    string s_weight = s_prefix + "wt"// total accumulation time at each point (arb. units)

    make /O/D/n=(n_theta) $s_index /wave=index
    make /O/D/n=(n_theta) $s_theta /wave=theta
    make /O/D/n=(n_theta) $s_dphi /wave=dphi
    make /O/D/n=(n_theta) $s_nphis /wave=nphis
    
    //---------- calculate phi-step-size for this theta:
    variable aux = Calc_The_step(Theta_in, Theta_st, HoloMode)
    dphi[0] =  calc_phi_step(Theta_in, Theta_in, aux, Phi_ran, Phi_ref, HoloMode)
    Theta[0] = Theta_in
    nphis[0] = calc_nth(aux, Theta_in, Theta_in, Phi_ran, Phi_ref, HoloMode)
    Index[0] = nphis[0]
    
    //---------- calculate number of phis, phi-step, and starting-index for each theta:
    variable ii = 1
    do
        Theta[ii] = Theta[ii-1] - aux
        if(Theta[ii] <= Theta_in-Theta_ran)
            Theta[ii] = Theta_in-Theta_ran
        endif
        aux =  Calc_The_step(Theta[ii], Theta_st, HoloMode)
        dphi[ii] = calc_phi_step(Theta_in, Theta[ii],  aux, Phi_ran, Phi_ref, HoloMode)
        nphis[ii] = calc_nth(aux, Theta_in, Theta[ii], Phi_ran, Phi_ref, HoloMode)
        Index[ii] = Index[ii-1] + nphis[ii]
        ii=ii+1
    while(ii < n_theta)

    if (Index[n_theta-1]==Index[n_theta-2])
        Index[n_theta-1]=Index[n_theta-2]+1
        nphis[n_theta-1]=1
    endif

    variable NumPoints = sum(nphis, 0, numpnts(nphis))
    
    //---------- calculate theta and phi for each data point:
    make /O/D/N=(NumPoints) $s_polar /wave=polar, $s_azim /wave=azim
    note azim, "version=1.6"
    
    ii = 0
    variable StartIndex = 0
    variable EndIndex
    do
        EndIndex=Index[ii]
        Polar[StartIndex, EndIndex-1]=Theta[ii]
        Azim[StartIndex, EndIndex-1]= mod(Phi_ran+(x-StartIndex)*dphi[ii]+Phi_in,Phi_ran)
        ii = ii + 1
        StartIndex = EndIndex
    while(ii < n_theta)

    duplicate /o azim, $s_int /wave=values
    duplicate /o azim, $s_tot /wave=totals
    duplicate /o azim, $s_weight /wave=weights
    values = nan
    totals = 0
    weights = 0

    // XPDplot metadata
    if (xpdplot)
        string s_FileName = ""
        string s_Comment = "created by pearl-anglescan-process.ipf"
        string s_HoloMode = "Stereographic"
        variable /g gb_SpectraFile = 0
    
        Make/O/D/n=22 $s_HoloData /wave=HoloData
        HoloData[0] = NaN// v_StartKE
        HoloData[1] = NaN// v_StoppKE
        HoloData[6] = NumPoints
        HoloData[7] = Theta_in
        HoloData[8] = Theta_ran
        HoloData[9] = Theta_st
        HoloData[11] = Phi_in
        HoloData[12] = Phi_ran
        HoloData[13] = Theta_st
        HoloData[15] = Phi_ref
        HoloData[16] = Phi_ran
        HoloData[17] = 0// v_HoloBit (stereographic)

        Make/O/T/n=22 $s_HoloInfo /wave=HoloInfo
        HoloInfo[0] = s_FileName
        HoloInfo[1] = s_Comment
        HoloInfo[10] = s_HoloMode
        HoloInfo[11] = ""// s_MeasuringMode

        // notebook for XPDplot
        if (WinType(NickName) == 5)
            Notebook $NickName selection={startOfFile, endOfFile}
            Notebook $NickName text=""
        else
            NewNotebook /F=0 /K=1 /N=$NickName /W=(5,40,341,260)
            Notebook $NickName defaultTab=140
            Notebook $NickName statusWidth=300
            Notebook $NickName backRGB=(56797,56797,56797)
            Notebook $NickName pageMargins={80,80,80,80}
            Notebook $NickName fSize=10
            Notebook $NickName fStyle=0,textRGB=(65535,0,26214)
            Notebook $NickName textRGB=(65535,0,26214)
        endif
        Notebook $NickName text = "File:\t" + s_FileName + "\r"
        Notebook $NickName text = "*** " + s_Comment + " ***\r\r"
        Notebook $NickName text = "Angle-Mode:\t" + s_HoloMode + "\r"
        Notebook $NickName text = "XPDplot Nickname:\t" + NickName + "\r"
    endif
};

string get_hemi_nickname(wave w){
    wave w
    
    string prefix = get_hemi_prefix(w)
    string wname = nameofwave(w)
    string nickname
    
    if (strlen(prefix))
        nickname = prefix
    else
        string s_wave_df = GetWavesDataFolder(w, 1)
        dfref   parent_df = $(s_wave_df + "::")
        nickname = GetDataFolder(0, parent_df)
    endif
    
    return nickname
};

string get_hemi_prefix(wave w){
    wave w

    string wname = nameofwave(w)
    string prefix
    if (ItemsInList(wname, "_") >= 2)
        prefix = StringFromList(0, wname, "_")
    else
        prefix = ""
    endif
    
    return prefix
};

dfr find_hemi_data(string nickname, string* prefix, string* intwave){
    string nickname
    string &prefix
    string &intwave

    dfref datadf
    prefix = ""
    intwave = "values"
    if (strlen(nickname))
        if (DataFolderExists(nickname))
            datadf = $nickname
        else
            datadf = getdatafolderdfr()
            prefix = nickname + "_"
            intwave = prefix + "i"
            if (exists(intwave) != 1)
                datadf = root:
            endif
        endif
    else
        datadf = getdatafolderdfr()
        prefix = ""
        intwave = "values"
    endif
    return datadf
};

variable clear_hemi_grid(string nickname){
    string nickname

    dfref datadf
    string s_prefix
    string s_int
    datadf = find_hemi_data(nickname, s_prefix, s_int)

    string s_totals = s_prefix + "tot"
    string s_weights = s_prefix + "wt"
    
    wave /sdfr=datadf /z w_values = $s_int
    wave /sdfr=datadf /z w_totals = $s_totals
    wave /sdfr=datadf /z w_weights = $s_weights

    if (waveexists(w_totals))
        w_totals = 0
    endif
    if (waveexists(w_weights))
        w_weights = 0
    endif
    if (waveexists(w_values))
        w_values = nan
    endif
};

variable duplicate_hemi_scan(string source_nickname, dfref dest_folder, string dest_nickname, variable xpdplot = defaultValue){
    string source_nickname
    dfref dest_folder
    string dest_nickname
    variable xpdplot
    
    if (ParamIsDefault(xpdplot))
        xpdplot = 0
    endif
    
    dfref savedf = getdatafolderdfr()

    // source data
    string s_prefix = ""
    string s_int = "values"
    dfref source_df = find_hemi_data(source_nickname, s_prefix, s_int)
    string s_polar = s_prefix + "pol"
    string s_azim = s_prefix + "az"
    string s_theta = s_prefix + "th"
    string s_tot = s_prefix + "tot"
    string s_weight = s_prefix + "wt"
    string s_matrix = s_prefix + "matrix"

    wave /sdfr=source_df theta1 = $s_theta
    wave /sdfr=source_df polar1 = $s_polar
    wave /sdfr=source_df azim1 = $s_azim
    wave /sdfr=source_df tot1 = $s_tot
    wave /sdfr=source_df weight1 = $s_weight
    wave /sdfr=source_df values1 = $s_int
    wave /sdfr=source_df /z matrix1 = $s_matrix
    
    variable npol = numpnts(theta1)
    
    setdatafolder dest_folder
    make_hemi_grid(npol, dest_nickname, xpdplot=xpdplot)
    
    // dest data
    dfref dest_df = find_hemi_data(dest_nickname, s_prefix, s_int)
    s_polar = s_prefix + "pol"
    s_azim = s_prefix + "az"
    s_theta = s_prefix + "th"
    s_tot = s_prefix + "tot"
    s_weight = s_prefix + "wt"
    s_matrix = s_prefix + "matrix"

    wave /sdfr=dest_df theta2 = $s_theta
    wave /sdfr=dest_df polar2 = $s_polar
    wave /sdfr=dest_df azim2 = $s_azim
    wave /sdfr=dest_df tot2 = $s_tot
    wave /sdfr=dest_df weight2 = $s_weight
    wave /sdfr=dest_df values2 = $s_int

    tot2 = tot1
    weight2 = weight1
    values2 = values1
    if (waveexists(matrix1))
        setdatafolder dest_df
        duplicate /o matrix1, $s_matrix
    endif
    
    if (!(NumberByKey("version", note(azim1), "=", "\r") >= 1.6))
        azim2 += 180// changed 151030 (v1.6)
        azim2 = azim2 >= 360 ? azim2 - 360 : azim2
    endif
    
    setdatafolder saveDF
};

variable rotate_hemi_scan(string nickname, variable angle){
    string nickname
    variable angle
    
    dfref savedf = getdatafolderdfr()

    string s_prefix = ""
    string s_int = "values"
    dfref df = find_hemi_data(nickname, s_prefix, s_int)

    string s_polar = s_prefix + "pol"
    string s_azim = s_prefix + "az"
    string s_tot = s_prefix + "tot"
    string s_weight = s_prefix + "wt"
    
    wave /sdfr=df polar = $s_polar
    wave /sdfr=df azim = $s_azim
    wave /sdfr=df tot = $s_tot
    wave /sdfr=df weight = $s_weight
    wave /sdfr=df values = $s_int

    azim += angle
    azim = azim < 0 ? azim + 360 : azim
    azim = azim >= 360 ? azim - 360 : azim

    duplicate /free polar, neg_polar
    neg_polar = -polar
    sort {neg_polar, azim}, polar, azim, tot, weight, values
    
    setdatafolder saveDF
};

string display_hemi_scan(string nickname, variable projection = defaultValue, variable graphtype = defaultValue, variable do_ticks = defaultValue, variable do_grids = defaultValue, string graphname = defaultValue){
    string nickname
    variable projection
    variable graphtype
    variable do_ticks
    variable do_grids
    string graphname
    
    dfref savedf = getdatafolderdfr()

    if (ParamIsDefault(projection))
        projection = 1
    endif
    if (ParamIsDefault(graphtype))
        graphtype = 1
    endif
    if (ParamIsDefault(do_ticks))
        do_ticks = 3
    endif
    if (ParamIsDefault(do_grids))
        do_grids = 3
    endif
    if (ParamIsDefault(graphname))
        if (strlen(nickname) > 0)
            graphname = nickname
        else
            graphname = GetDataFolder(0)
        endif
    endif
    
    // hemi grid waves
    string s_prefix = ""
    string s_int = "values"
    dfref df = find_hemi_data(nickname, s_prefix, s_int)

    string s_polar = s_prefix + "pol"
    string s_azim = s_prefix + "az"
    string s_matrix = s_prefix + "matrix"
    
    wave /sdfr=df /z values = $s_int
    wave /sdfr=df /z azim = $s_azim
    wave /sdfr=df /z polar = $s_polar
    wave /sdfr=df /z matrix = $s_matrix

    setdatafolder df
    string s_ster_rad = s_prefix + "ster_rad"
    duplicate /o polar, $s_ster_rad /wave=ster_rad
    ster_rad = calc_graph_radius(polar, projection=projection)
    
    string s_ster_x = s_prefix + "ster_x"
    string s_ster_y = s_prefix + "ster_y"
    duplicate /o azim, $s_ster_x /wave=ster_x, $s_ster_y /wave=ster_y
    ster_x = ster_rad * cos(azim * pi / 180)
    ster_y = ster_rad * sin(azim * pi / 180)

    variable azim_offset = 0
    if (!(NumberByKey("version", note(azim), "=", "\r") >= 1.6))
        DoAlert /T="display hemi scan" 0, "your dataset doesn't include the version 1.6 flag. if it was created with an earlier version that might be okay. please check that the orientation is correct!"
        azim_offset = 180// changed 151030 (v1.6)
    endif

    string s_trace
    switch(graphtype)
        case 1:
            graphname = display_polar_graph(graphname, angle_offset=azim_offset, do_ticks=do_ticks)

            s_trace = WMPolarAppendTrace(graphname, ster_rad, azim, 360)
            ModifyGraph /W=$graphname mode($s_trace)=2, lsize($s_trace)=2
            ModifyGraph /W=$graphname zColor($s_trace)={values,*,*,BlueGreenOrange,0}
            
            ColorScale /W=$graphname /C /N=text0 /E=2 /F=0 /B=1 /A=RB /X=0.00 /Y=0.00 trace=polarY0
            ColorScale /W=$graphname /C /N=text0 side=2, width=5, heightPct=40, frame=0.50, lblMargin=0
            ColorScale /W=$graphname /C /N=text0 nticks=2, minor=1, tickLen=4.00, tickThick=0.50

            SetWindow $graphname, userdata(projection)=num2str(projection)
            draw_hemi_axes(graphname, do_grids=do_grids)
            break
        case 3:
            graphname = display_polar_graph(graphname, angle_offset=azim_offset, do_ticks=do_ticks)

            s_trace = WMPolarAppendTrace(graphname, ster_rad, azim, 360)
            ModifyGraph /W=$graphname mode($s_trace)=0, lsize($s_trace)=0
            AppendImage /L=VertCrossing /B=HorizCrossing matrix

            ColorScale /W=$graphname /C /N=text0 /E=2 /F=0 /B=1 /A=RB /X=0.00 /Y=0.00 image=$s_matrix
            ColorScale /W=$graphname /C /N=text0 side=2, width=5, heightPct=40, frame=0.50, lblMargin=0
            ColorScale /W=$graphname /C /N=text0 nticks=2, minor=1, tickLen=4.00, tickThick=0.50
        
            SetWindow $graphname, userdata(projection)=num2str(projection)
            draw_hemi_axes(graphname, do_grids=do_grids)
            break
    endswitch
    
    setdatafolder savedf
    return graphname
};

static string display_polar_graph(string graphname, variable angle_offset = defaultValue, variable do_ticks = defaultValue){
    
    string graphname
    variable angle_offset
    variable do_ticks
    
    dfref savedf = GetDataFolderDFR()
    
    if (ParamIsDefault(angle_offset))
        angle_offset = 0
    endif
    if (ParamIsDefault(do_ticks))
        do_ticks = 3
    endif
    
    if ((strlen(graphname) == 0) || (wintype(graphname) == 0))
        Display /k=1 /W=(10,45,360,345)
        DoWindow /C $graphname
        graphname = WMNewPolarGraph("", graphname)
        WMPolarGraphSetVar(graphname, "zeroAngleWhere", angle_offset)
        
        WMPolarGraphSetVar(graphname, "angleAxisThick", 0.5)
        WMPolarGraphSetStr(graphname, "doMajorAngleTicks", "manual")
        WMPolarGraphSetVar(graphname, "majorAngleInc", 30)// major ticks in 30 deg steps
        WMPolarGraphSetVar(graphname, "minorAngleTicks", 2)// minor ticks in 10 deg steps
        WMPolarGraphSetStr(graphname, "angleTicksLocation", "Outside")
        WMPolarGraphSetVar(graphname, "doAngleTickLabelSubRange", 1)
        WMPolarGraphSetVar(graphname, "angleTickLabelRangeStart", 0)
        WMPolarGraphSetVar(graphname, "angleTickLabelRangeExtent", 90)
        WMPolarGraphSetStr(graphname, "angleTickLabelNotation", "%g")
  
        WMPolarGraphSetVar(graphname, "doPolarGrids", 0)
        WMPolarGraphSetVar(graphname, "doRadiusTickLabels", 0)
        WMPolarGraphSetStr(graphname, "radiusAxesWhere", "  Off")// note the leading spaces, cf. WMPolarAnglesForRadiusAxes
        WMPolarGraphSetStr(graphname, "radiusTicksLocation", "Off")

        WMPolarGraphSetVar(graphname, "majorTickLength", 2)
        WMPolarGraphSetVar(graphname, "majorTickThick", 0.5)
        WMPolarGraphSetVar(graphname, "minorTickLength", 1)
        WMPolarGraphSetVar(graphname, "minorTickThick", 0.5)
        WMPolarGraphSetVar(graphname, "tickLabelOpaque", 0)
        WMPolarGraphSetVar(graphname, "tickLabelFontSize", 7)

        // changes
        if (do_ticks & 1)
            WMPolarGraphSetStr(graphname, "angleTicksLocation", "Outside")
        else
            WMPolarGraphSetStr(graphname, "angleTicksLocation", "Off")
        endif
        if (do_ticks & 2)
            WMPolarGraphSetVar(graphname, "doMinorAngleTicks", 1)
        else
            WMPolarGraphSetVar(graphname, "doMinorAngleTicks", 0)
        endif
        
        DoWindow /T $graphname, graphname

        // cursor info in angles
        string graphdf = "root:packages:WMPolarGraphs:" + graphname
        setdatafolder graphdf
        // current theta, phi coordinates are stored in global variables in the package folder of the graph
        variable /g csrA_theta
        variable /g csrA_phi
        variable /g csrB_theta
        variable /g csrB_phi
        // the text box is hidden initially. it shows up and hides with the cursor info box.
        string tb
        tb = "\\{"
        tb = tb + "\"A = (%.1f, %.1f)\","
        tb = tb + graphdf + ":csrA_theta,"
        tb = tb + graphdf + ":csrA_phi"
        tb = tb + "}"
        TextBox /W=$graphname /A=LT /B=1 /E=2 /F=0 /N=tb_angles /X=1 /Y=1 /V=0 tb
        tb = "\\{"
        tb = tb + "\"B = (%.1f, %.1f)\","
        tb = tb + graphdf + ":csrB_theta,"
        tb = tb + graphdf + ":csrB_phi"
        tb = tb + "}"
        AppendText /W=$graphname /N=tb_angles tb
        // updates are triggered by a window hook
        SetWindow $graphname, hook(polar_graph_hook)=PearlAnglescanProcess#polar_graph_hook
    else
        // graph window exists
        DoWindow /F $graphname
    endif

    setdatafolder savedf
    return graphname
};

static string draw_hemi_axes(string graphname, variable do_grids = defaultValue){
    string graphname
    variable do_grids
    
    if (ParamIsDefault(do_grids))
        do_grids = 3
    endif

    dfref savedf = GetDataFolderDFR()

    string sproj = GetUserData(graphname, "", "projection")
    variable projection = str2num("0" + sproj)

    SetDrawLayer /W=$graphname ProgFront

    // polar axis
    SetDrawEnv /W=$graphname xcoord=HorizCrossing, ycoord=VertCrossing
    SetDrawEnv /W=$graphname linethick= 0.5
    SetDrawEnv /W=$graphname dash=2
    SetDrawEnv /W=$graphname fillpat=0
    SetDrawEnv /W=$graphname fname="default", fsize=7
    SetDrawEnv /W=$graphname textxjust=1, textyjust=1
    //SetDrawEnv /W=$graphname linefgc=(65535,65535,65535)
    SetDrawEnv /W=$graphname save

    if (do_grids & 1)
        DrawLine /W=$graphname 0, -2, 0, 2
        DrawLine /W=$graphname -2, 0, 2, 0
    endif

    variable radi
    if (do_grids & 2)
        radi = calc_graph_radius(0.5, projection=projection)
        DrawOval /W=$graphname -radi, radi, radi, -radi
        radi = calc_graph_radius(30, projection=projection)
        DrawOval /W=$graphname -radi, radi, radi, -radi
        radi = calc_graph_radius(60, projection=projection)
        DrawOval /W=$graphname -radi, radi, radi, -radi

        SetDrawEnv /W=$graphname textxjust= 1,textyjust= 2
        SetDrawEnv /W=$graphname save
        radi = calc_graph_radius(30, projection=projection)
        DrawText /W=$graphname radi, -0.1, "30"
        radi = calc_graph_radius(60, projection=projection)
        DrawText /W=$graphname radi, -0.1, "60"
    endif
    
    setdatafolder savedf
};

variable draw_diffraction_cone(string graphname, string groupname, variable theta_axis, variable theta_inner, variable phi){
    string graphname
    string groupname
    
    variable theta_axis
    variable theta_inner
    variable phi
    
    variable r_axis = calc_graph_radius(theta_axis)
    variable r_inner = calc_graph_radius(theta_inner)
    variable r_outer = calc_graph_radius(2 * theta_axis - theta_inner)
    
    SetDrawEnv push
    SetDrawLayer UserFront
    DrawAction getgroup=$groupname, delete
    SetDrawEnv gstart, gname=$groupname
    variable xc, yc, xr, yr

    // cone periphery
    variable r_center = (r_outer + r_inner) / 2
    variable r_radius = (r_outer - r_inner) / 2
    xc = r_center * cos(phi * pi / 180)
    yc = r_center * sin(phi * pi / 180)
    xr = r_radius
    yr = r_radius
    SetDrawEnv xcoord=HorizCrossing, ycoord=VertCrossing
    SetDrawEnv dash=11, fillpat=0
    DrawOval xc - xr, yc - yr, xc + xr, yc + yr

    // cone axis    
    xc = r_axis * cos(phi * pi / 180)
    yc = r_axis * sin(phi * pi / 180)
    r_radius = calc_graph_radius(2)
    xr = r_radius
    yr = r_radius
    SetDrawEnv xcoord=HorizCrossing, ycoord=VertCrossing
    SetDrawEnv fillfgc=(0,0,0)
    DrawOval xc - xr, yc - yr, xc + xr, yc + yr
    
    SetDrawEnv gstop
    SetDrawEnv pop
};

string display_scanlines(string nickname, variable alpha_lo, variable alpha_hi, wave m_theta, wave m_tilt, wave m_phi, variable folding = defaultValue, variable projection = defaultValue){
    string nickname
    variable alpha_lo
    variable alpha_hi
    wave m_theta
    wave m_tilt
    wave m_phi
    variable folding
    variable projection
    
    if (ParamIsDefault(folding))
        folding = 1
    endif
    if (ParamIsDefault(projection))
        projection = 1
    endif
    
    // sort out data folder structure
    dfref saveDF = GetDataFolderDFR()
    newdatafolder /s/o $nickname
    string graphname = "graph_" + nickname
    
    duplicate /free m_tilt, loc_m_tilt
    loc_m_tilt = -m_tilt
    
    make /n=1 /d /free d_polar, d_azi
    variable n_alpha = round(alpha_hi - alpha_lo) + 1
    make /n=(n_alpha) /d /free analyser
    setscale /i x alpha_lo, alpha_hi, "", analyser
    analyser = x
    
    convert_angles_ttpa2polar(m_theta, loc_m_tilt, m_phi, analyser, d_polar, d_azi)
    duplicate /free d_polar, d_radius
    d_radius = calc_graph_radius(d_polar, projection=projection)
    d_azi += 180// changed 151030 (v1.6)
    
    graphname = display_polar_graph(graphname)
    SetWindow $graphname, userdata(projection)=num2str(projection)
    
    variable ifold
    variable iang
    variable nang = numpnts(m_theta)
    string s_rad
    string s_azi
    string s_trace
    for (ifold = 0; ifold < folding; ifold += 1)
        d_azi = d_azi >= 360 ? d_azi - 360 : d_azi
        for (iang = 0; iang < nang; iang += 1)
            sprintf s_rad, "rad_%d_%d", ifold, iang
            duplicate /o analyser, $s_rad
            wave w_rad = $s_rad
            w_rad = d_radius[p][iang]
            
            sprintf s_azi, "azi_%d_%d", ifold, iang
            duplicate /o analyser, $s_azi
            wave w_azi = $s_azi
            w_azi = d_azi[p][iang]
            
            if (numtype(sum(w_rad)) == 0)
                s_trace = WMPolarAppendTrace(graphname, w_rad, w_azi, 360)
                ModifyGraph /w=$graphname mode($s_trace)=0, lsize($s_trace)=0.5
            endif
        endfor
        d_azi += 360 / folding
    endfor

    draw_hemi_axes(graphname)

    setdatafolder saveDF
    return graphname
};

const variable kProjDist = 0;
const variable kProjStereo = 1;
const variable kProjArea = 2;
const variable kProjGnom = 3;
const variable kProjOrtho = 4;

static const variable kProjScaleDist = 2;
static const variable kProjScaleStereo = 2;
static const variable kProjScaleArea = 2;
// scaled so that radius(gnom) = radius(stereo) for polar = 88
static const variable kProjScaleGnom = 0.06744519021;
static const variable kProjScaleOrtho = 2;

threadsafe variable calc_graph_radius(variable polar, variable projection = defaultValue){
    variable polar
    variable projection
    
    if (ParamIsDefault(projection))
        projection = 1
    endif

    variable radius
    switch(projection)
        case kProjStereo:// stereographic
            radius = kProjScaleStereo * tan(polar / 2 * pi / 180)
            break
        case kProjArea:// equal area
            radius = kProjScaleArea * sin(polar / 2 * pi / 180)
            break
        case kProjGnom:// gnomonic
            radius = polar < 90 ? kProjScaleGnom * tan(polar * pi / 180) : inf
            break
        case kProjOrtho:// orthographic
            radius = kProjScaleOrtho * sin(polar * pi / 180)
            break
        default:// equidistant
            radius = kProjScaleDist * polar / 90
    endswitch
    
    return radius
};

threadsafe variable calc_graph_polar(variable x, variable y, variable projection = defaultValue){
    variable x
    variable y
    variable projection
    
    if (ParamIsDefault(projection))
        projection = 1
    endif

    variable radius
    variable polar
    
    radius = sqrt(x^2 + y^2)
    switch(projection)
        case kProjStereo:// stereographic
            polar = 2 * atan(radius / kProjScaleStereo) * 180 / pi
            break
        case kProjArea:// equal area
            polar = 2 * asin(radius / kProjScaleArea) * 180 / pi
            break
        case kProjGnom:// gnomonic
            polar = atan(radius / kProjScaleGnom) * 180 / pi
            break
        case kProjOrtho:// orthographic
            polar = asin(radius / kProjScaleOrtho) * 180 / pi
            break
        default:// equidistant
            polar = 90 * radius / kProjScaleDist
    endswitch

    return polar
};

threadsafe variable calc_graph_azi(variable x, variable y, variable projection = defaultValue, variable zeroAngle = defaultValue){
    variable x
    variable y
    variable projection
    variable zeroAngle
    
    if (ParamIsDefault(projection))
        projection = 1
    endif
    if (ParamIsDefault(zeroAngle))
        zeroAngle = 0
    endif

    variable azi
    if (x > 0)
        azi = atan(y / x) * 180 / pi
    else
        azi = atan(y / x) * 180 / pi + 180
    endif
    
    azi += zeroAngle
    if (azi < 0)
        azi += 360
    endif
    if (azi >= 360)
        azi -= 360
    endif
    if (numtype(azi) != 0)
        azi = 0
    endif
    
    return azi
};

static variable update_polar_info(string graphname){
    string graphname
    
    dfref savedf = GetDataFolderDFR()

    string graphdf = "root:packages:WMPolarGraphs:" + graphname
    setdatafolder graphdf
    
    nvar csrA_theta
    nvar csrA_phi
    nvar csrB_theta
    nvar csrB_phi

    string sproj = GetUserData(graphname, "", "projection")
    variable projection = str2num("0" + sproj)
    nvar zeroAngleWhere
    
    variable x = hcsr(A, graphname)
    variable y = vcsr(A, graphname)
    csrA_theta = calc_graph_polar(x, y, projection=projection)
    csrA_phi = calc_graph_azi(x, y, projection=projection, zeroAngle=zeroAngleWhere)

    x = hcsr(B, graphname)
    y = vcsr(B, graphname)
    csrB_theta = calc_graph_polar(x, y, projection=projection)
    csrB_phi = calc_graph_azi(x, y, projection=projection, zeroAngle=zeroAngleWhere)
    
    setdatafolder savedf
};

static variable polar_graph_hook(WMWinHookStruct* s){
    STRUCT WMWinHookStruct &s

    Variable hookResult = 0

    switch(s.eventCode)
        case 7:// cursor moved
            update_polar_info(s.winname)
            break
        case 20:// show info
            TextBox /W=$s.winname /N=tb_angles /C /V=1
            break
        case 21:// hide info
            TextBox /W=$s.winname /N=tb_angles /C /V=0
            break
    endswitch

    return hookResult// 0 if nothing done, else 1
};

variable set_polar_graph_cursor(string nickname, string cursorname, variable polar_angle, variable azim_angle, string graphname = defaultValue){
    string nickname
    string cursorname
    variable polar_angle
    variable azim_angle
    string graphname
    
    if (ParamIsDefault(graphname))
        if (strlen(nickname) > 0)
            graphname = nickname
        else
            graphname = GetDataFolder(0)
        endif
    endif
    
    string s_prefix = ""
    string s_int = "values"
    dfref df = find_hemi_data(nickname, s_prefix, s_int)

    string s_polar = s_prefix + "pol"
    string s_azim = s_prefix + "az"
    wave /sdfr=df /z azim = $s_azim
    wave /sdfr=df /z polar = $s_polar

    FindLevel /P /Q polar, polar_angle
    if (v_flag == 0)
        variable polar_level = floor(v_levelx)
        FindLevel /P /Q /R=[polar_level] azim, azim_angle
        if (v_flag == 0)
            variable azim_level = round(v_levelx)
            string tracename = "polarY0"
            Cursor /W=$graphname /P $cursorname $traceName azim_level
        endif
    endif
};

variable hemi_add_anglescan(string nickname, wave values, wave polar, wave azi, wave weights = defaultValue){
    string nickname// name prefix of holo waves.
        // may be empty.
    wave values// intensity values
        // the wave can be one- or two-dimensional.
        // no specific order required, the function sorts the arrays internally
    wave polar// polar coordinates. allowed range 0 <= theta <= 90
        // dimensions corresponding to value.
    wave azi// azimuthal coordinates. allowed range -360 <= phi < +360
        // dimensions corresponding to value.
    wave weights// total accumulation time of each point of values. default = 1

    if (ParamIsDefault(weights))
        duplicate /free values, weights
        weights = 1
    endif
    
    // quick check whether hemi grid is existing
    string s_prefix = ""
    string s_int = "values"
    dfref df = find_hemi_data(nickname, s_prefix, s_int)

    string s_polar = s_prefix + "pol"
    string s_azim = s_prefix + "az"
    string s_theta = s_prefix + "th"

    wave /sdfr=df /z w_values = $s_int
    wave /sdfr=df /z w_azim = $s_azim
    wave /sdfr=df /z w_polar = $s_polar
    wave /sdfr=df /z w_theta = $s_theta
    if (!waveexists(w_values) || !waveexists(w_azim) || !waveexists(w_polar))
        abort "Missing hemispherical scan grid. Please call make_hemi_grid() first."
    endif
    
    // make internal copies, one-dimensional, ordered in theta
    duplicate /free values, values_copy
    duplicate /free polar, polar_copy
    duplicate /free azi, azi_copy
    duplicate /free weights, weights_copy
    variable nn = dimsize(values, 0) * max(dimsize(values, 1), 1)
    redimension /n=(nn) values_copy, polar_copy, azi_copy, weights_copy
    sort /r polar_copy, polar_copy, azi_copy, values_copy, weights_copy
    
    variable pol
    variable pol_st = abs(w_theta[1] - w_theta[0])
    variable pol1, pol2
    
    duplicate /free azi_copy, azi_slice
    duplicate /free values_copy, values_slice
    duplicate /free weights_copy, weights_slice
    for (pol = 90; pol >= 0; pol -= pol_st)
        pol1 = pol - pol_st / 2
        pol2 = pol + pol_st / 2
        extract /free /indx polar_copy, sel, (pol1 < polar_copy) && (polar_copy <= pol2)
        if (numpnts(sel) > 0)
            redimension /n=(numpnts(sel)) azi_slice, values_slice, weights_slice
            azi_slice = azi_copy[sel]
            values_slice = values_copy[sel]
            weights_slice = weights_copy[sel]
            hemi_add_aziscan(nickname, values_slice, pol, azi_slice, weights=weights_slice)
        endif
    endfor
};

variable hemi_add_aziscan(string nickname, wave values, variable polar, wave azi, wave weights = defaultValue){
    string nickname// name prefix of holo waves.
        // may be empty.
    wave values// intensity values of the azimuthal scan at the positions given in the azi parameter
    variable polar// polar angle where to add the azi scan
    wave azi// angle positions of the azimuthal scan
        // acceptable range: >= -360 and < +360
        // no specific order required, the function sorts the array internally
    wave weights// total accumulation time of each point of values. default = 1

    if (ParamIsDefault(weights))
        duplicate /free values, weights
        weights = 1
    endif
        
    // hemi grid waves
    string s_prefix = ""
    string s_int = "values"
    dfref df = find_hemi_data(nickname, s_prefix, s_int)

    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"
    
    wave /sdfr=df w_polar = $s_polar
    wave /sdfr=df w_azim = $s_azim
    wave /sdfr=df w_values = $s_int
    wave /sdfr=df w_totals = $s_totals
    wave /sdfr=df w_weights = $s_weights
    wave /sdfr=df w_index = $s_index
    wave /sdfr=df w_theta = $s_theta
    wave /sdfr=df w_dphi = $s_dphi
    wave /sdfr=df w_nphis = $s_nphis
    
    // destination slice coordinates
    //polar = round(polar)
    //variable ipol = 90 - polar
    variable ipol = BinarySearch(w_theta, polar)
    if (ipol < 0)
        abort "assertion failed in hemi_add_aziscan(): polar angle not found in grid."
    endif
    
    variable d1, d2
    if (ipol >= 1)
        d1 = w_index[ipol - 1]
    else
        d1 = 0
    endif
    d2 = w_index[ipol] - 1
    variable nd = d2 - d1 + 1
    variable dphi = w_dphi[ipol]
    variable az1, az2
    
    // source slice coordinates 
    // order the slice from -dphi/2 to 360-dphi/2
    azi = azi < 0 ? azi + 360 : azi
    azi = azi >= 360 - dphi/2 ? azi - 360 : azi
    duplicate /free values, sel_values
    duplicate /free weights, sel_weights
    
    // loop over destination
    variable id
    variable v1, v2, w1, w2
    for (id = 0; id < nd; id += 1)
        az1 = (id - 0.5) * dphi
        az2 = (id + 0.5) * dphi
        extract /free /indx azi, sel, (az1 <= azi) && (azi < az2)
        if (numpnts(sel) > 0)
            redimension /n=(numpnts(sel)) sel_values, sel_weights
            sel_values = values[sel]
            sel_weights = weights[sel]
            v1 = w_totals[d1 + id]
            w1 = w_weights[d1 + id]
            if ((numtype(v1) == 2) || (w1 <= 0))
                v1 = 0
                w1 = 0
            endif
            v2 = sum(sel_values)
            w2 = sum(sel_weights)
            w_totals[d1 + id] = v1 + v2
            w_weights[d1 + id] = w1 + w2
        endif
    endfor
    w_values[d1, d1 + nd - 1] = w_totals[p] / w_weights[p]
};

variable interpolate_hemi_scan(string nickname){
    string nickname

    dfref savedf = GetDataFolderDFR()
    
    string s_prefix = ""
    string s_int = "values"
    dfref df = find_hemi_data(nickname, s_prefix, s_int)

    string s_polar = s_prefix + "pol"
    string s_azim = s_prefix + "az"
    string s_ster_x = s_prefix + "ster_x"
    string s_ster_y = s_prefix + "ster_y"

    wave /sdfr=df values = $s_int
    wave /sdfr=df azim = $s_azim
    wave /sdfr=df polar = $s_polar
    wave /sdfr=df ster_x = $s_ster_x
    wave /sdfr=df ster_y = $s_ster_y
    
    variable min_ster_x = wavemin(ster_x)
    variable max_ster_x = wavemax(ster_x)
    variable x0 = min_ster_x
    variable xn = 181
    variable dx = (max_ster_x - min_ster_x) / (xn - 1)
    make /n=(numpnts(ster_x), 3) /free triplet
    triplet[][0] = ster_x[p]
    triplet[][1] = ster_y[p]
    triplet[][2] = values[p]

    variable size = 181
    setdatafolder df
    make /n=(size, size) /d /o $(s_prefix + "matrix") /wave=matrix
    make /n=(size, size) /free mnorm
    ImageFromXYZ /as {ster_x, ster_y, values}, matrix, mnorm
    matrix /= mnorm
    matrixfilter NanZapMedian, matrix
    matrixfilter gauss, matrix
    
    duplicate /free values, ster_finite
    ster_finite = (numtype(values) == 0) * (ster_x^2 + ster_y^2)
    variable ster_max = wavemax(ster_finite)
    matrix = (x^2 + y^2) <= ster_max ? matrix : nan
    
    setdatafolder savedf
};

variable quick_pizza_image(wave data, string nickname, variable theta_offset, variable tilt_offset, variable phi_offset, variable npolar = defaultValue, variable nograph = defaultValue, variable folding = defaultValue){
    wave data// 2D intensity wave, see requirements above
    string nickname// nick name for output data
        // in default mode, this will be the name of a child folder containing the output
        // in XPDplot mode, this will be a prefix of the generated data in the root folder
    variable theta_offset// manipulator theta angle corresponding to normal emission
    variable tilt_offset// manipulator tilt angle corresponding to normal emission
    variable phi_offset// manipulator phi angle corresponding to phi_result = 0
    variable npolar// number of polar angles, determines polar and azimuthal step size
        // default = 91 (1 degree steps)
    variable nograph// 0 (default) = display a new polar graph
        // 1 = don't display a new graph (if a graph is existing from a previous call, it will update)
    variable folding// rotational averaging, default = 1
        
    if (ParamIsDefault(npolar))
        npolar = 91
    endif
    if (ParamIsDefault(nograph))
        nograph = 0
    endif
    if (ParamIsDefault(folding))
        folding = 1
    endif
    string graphname = "graph_" + nickname
    string s_prefix = ""
    
    // sort out data folder structure
    dfref saveDF = GetDataFolderDFR()
    dfref dataDF = GetWavesDataFolderDFR(data)
    setdatafolder dataDF
    if (DataFolderExists(":attr"))
        setdatafolder :attr
    endif
    dfref attrDF = GetDataFolderDFR()
    setdatafolder dataDF
    newdatafolder /s/o $nickname
    dfref destDF = GetDataFolderDFR()
    
    // performance monitoring
    variable timerRefNum
    variable /g xyz_perf_secs
    timerRefNum = startMSTimer

    wave /sdfr=attrDF ManipulatorTheta
    wave /sdfr=attrDF ManipulatorTilt
    wave /sdfr=attrDF ManipulatorPhi
    duplicate /free ManipulatorTheta, m_theta
    duplicate /free ManipulatorTilt, m_tilt
    duplicate /free ManipulatorPhi, m_phi
    m_theta -= theta_offset
    m_tilt -= tilt_offset
    m_tilt *= -1// checked 140702
    m_phi -= phi_offset
    //m_phi *= -1 // checked 140702
    
    make /n=1/d/free d_polar, d_azi
    convert_angles_ttpd2polar(m_theta, m_tilt, m_phi, data, d_polar, d_azi)
    d_azi += 180// changed 151030 (v1.6)
    d_azi = d_azi >= 360 ? d_azi - 360 : d_azi

    duplicate /free data, values
    variable nn = dimsize(values, 0) * max(dimsize(values, 1), 1)
    redimension /n=(nn) values, d_polar, d_azi
    duplicate /o d_polar, ster_rad, ster_x, ster_y

    variable projection = 1
    switch(projection)
        case 1:// stereographic
            ster_rad =  2 * tan(d_polar / 2 * pi / 180)
            break
        case 2:// azimuthal
            ster_rad = 2 * cos((180 - d_polar) / 2 * pi / 180)
            break
    endswitch
    string s_ster_x = s_prefix + "ster_x"
    string s_ster_y = s_prefix + "ster_y"
    
    nn = 401
    make /n=(nn, nn) /d /o matrix
    make /n=(nn, nn) /free mnorm
    setscale /i x -2, +2, matrix, mnorm
    setscale /i y -2, +2, matrix, mnorm
    matrix = 0
    mnorm = 0
    
    variable ifold
    for (ifold = 0; ifold < folding; ifold += 1)
        ster_x = ster_rad * cos(d_azi * pi / 180)
        ster_y = ster_rad * sin(d_azi * pi / 180)
        ImageFromXYZ {ster_x, ster_y, values}, matrix, mnorm
        d_azi = d_azi >= 180 ? d_azi + 360 / folding - 180 : d_azi + 360 / folding
    endfor
    
    matrix /= mnorm
    matrixfilter /n=5 NanZapMedian matrix
    matrixfilter /n=3 gauss matrix

    if (!nograph)
        display /k=1
        appendimage matrix
        modifygraph width={Plan,1,bottom,left}
    endif
    
    if (timerRefNum >= 0)
        xyz_perf_secs = stopMSTimer(timerRefNum) / 1e6
    endif

    setdatafolder saveDF
};

variable save_hemi_scan(string nickname, string pathname, string filename){
    string nickname
    string pathname
    string filename
    
    dfref savedf = getdatafolderdfr()

    // source data
    string s_prefix = ""
    string s_int = "values"
    dfref df = find_hemi_data(nickname, s_prefix, s_int)

    string s_polar = s_prefix + "pol"
    string s_azim = s_prefix + "az"
    string s_theta = s_prefix + "th"
    string s_tot = s_prefix + "tot"
    string s_weight = s_prefix + "wt"

    wave /sdfr=df theta1 = $s_theta
    wave /sdfr=df polar1 = $s_polar
    wave /sdfr=df azim1 = $s_azim
    wave /sdfr=df tot1 = $s_tot
    wave /sdfr=df weight1 = $s_weight
    wave /sdfr=df values1 = $s_int
    
    save /m="\r\n" /o /p=$pathname /t theta1, polar1, azim1, tot1, weight1, values1 as filename
    
    setdatafolder saveDF
};

variable load_hemi_scan(string nickname, string pathname, string filename){
    string nickname
    string pathname
    string filename
    
    dfref savedf = getdatafolderdfr()

    //loadwave /p=$pathname /t theta1, polar1, azim1, tot1, weight1, values1 as filename
    //LoadWave /t/p=pearl_explorer_filepath/q filename
    //svar waves = s_wavenames
    //if (v_flag > 0)
    //  string /g pearl_explorer_import = "load_itx_file"
    //endif
    
    setdatafolder saveDF
};

variable import_tpi_scan(string nickname, wave theta, wave phi, wave intensity, variable folding = defaultValue, variable npolar = defaultValue, variable nograph = defaultValue, variable xpdplot = defaultValue){
    string nickname
    wave theta
    wave phi
    wave intensity
    
    variable folding
    variable npolar
    variable nograph
    variable xpdplot
    
    if (ParamIsDefault(npolar))
        npolar = 91
    endif
    if (ParamIsDefault(nograph))
        nograph = 0
    endif
    if (ParamIsDefault(folding))
        folding = 1
    endif
    if (ParamIsDefault(xpdplot))
        xpdplot = 0
    endif
    
    make_hemi_grid(npolar, nickname, xpdplot=xpdplot)

    variable ifold
    duplicate /free phi, fold_phi
    for (ifold = 0; ifold < folding; ifold += 1)
        hemi_add_anglescan(nickname, intensity, theta, fold_phi)
        fold_phi = fold_phi >= 180 ? fold_phi + 360 / folding - fold_phi : fold_phi + 360 / folding
    endfor

    if (nograph==0)
        display_hemi_scan(nickname)
    endif
};

variable trim_hemi_scan(string nickname, variable theta_max){
    string nickname
    variable theta_max

    string s_prefix = ""
    string s_int = "values"
    dfref df = find_hemi_data(nickname, s_prefix, s_int)

    string s_totals = s_prefix + "tot"
    string s_weights = s_prefix + "wt"
    string s_polar = s_prefix + "pol"

    wave /sdfr=df w_polar = $s_polar
    wave /sdfr=df w_values = $s_int
    wave /sdfr=df w_totals = $s_totals
    wave /sdfr=df w_weights = $s_weights

    w_values = w_polar <= theta_max ? w_totals / w_weights : nan
};

wave hemi_polar_cut(string nickname, variable azim){
    string nickname
    variable azim

    dfref savedf = getdatafolderdfr()
    string s_prefix = ""
    string s_int = "values"
    dfref df = find_hemi_data(nickname, s_prefix, s_int)

    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 /sdfr=df w_polar = $s_polar
    wave /sdfr=df w_azim = $s_azim
    wave /sdfr=df w_values = $s_int
    wave /sdfr=df w_totals = $s_totals
    wave /sdfr=df w_weights = $s_weights
    wave /sdfr=df w_index = $s_index
    wave /sdfr=df w_theta = $s_theta
    wave /sdfr=df w_dphi = $s_dphi
    wave /sdfr=df 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
    
    setdatafolder df
    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
    
    setdatafolder savedf
    return w_cut
};

wave hemi_azi_cut(string nickname, variable pol){
    string nickname
    variable pol

    dfref savedf = getdatafolderdfr()
    string s_prefix = ""
    string s_int = "values"
    dfref df = find_hemi_data(nickname, s_prefix, s_int)

    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 /sdfr=df w_polar = $s_polar
    wave /sdfr=df w_azim = $s_azim
    wave /sdfr=df w_values = $s_int
    wave /sdfr=df w_totals = $s_totals
    wave /sdfr=df w_weights = $s_weights
    wave /sdfr=df w_index = $s_index
    wave /sdfr=df w_theta = $s_theta
    wave /sdfr=df w_dphi = $s_dphi
    wave /sdfr=df 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)
        setdatafolder df
        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]], "", w_cut
        setdatafolder savedf
        return w_cut
    else
        setdatafolder savedf
        return $""
    endif
    setdatafolder savedf
};

static variable check_contrast(wave values, variable pcmin, variable pcmax, variable* vmin, variable* vmax){
    wave values
    variable pcmin
    variable pcmax
    variable &vmin
    variable &vmax
    
    dfref save_df = GetDataFolderDFR()
    dfref dfr = NewFreeDataFolder()
    setdatafolder dfr
    StatsQuantiles /inan /iw /q /z values
    wave index = w_quantilesindex
    variable imin = round(numpnts(index) * pcmin / 100)
    variable imax = round(numpnts(index) * (100 - pcmax) / 100)
    vmin = values[index[imin]]
    vmax = values[index[imax]]
    KillDataFolder dfr
    setdatafolder save_df
};

variable set_contrast(variable pcmin, variable pcmax, string graphname = defaultValue, string colortable = defaultValue){
    variable pcmin
    variable pcmax
    string graphname
    string colortable
    
    if (ParamIsDefault(graphname))
        graphname = ""
    endif
    if (ParamIsDefault(colortable))
        colortable = ""
    endif
    
    dfref save_df = GetDataFolderDFR()
    
    string objname
    string info
    string wname
    string ctab
    variable rev
    variable n
    variable i
    variable vmin
    variable vmax

    string traces = TraceNameList(graphname, ";", 1+4)
    n = ItemsInList(traces, ";")
    for (i = 0; i < n; i += 1)
        objname = StringFromList(i, traces, ";")
        info = TraceInfo(graphname, objname, 0)
        if (strlen(info) > 0)
            info = StringByKey("RECREATION", info, ":", ";")
            info = StringByKey("zColor(x)", info, "=", ";")
            if (strlen(info) > 2)
                info = info[1,strlen(info)-2]
                wname = StringFromList(0, info, ",")
                wave w = $wname
                ctab = StringFromList(3, info, ",")
                rev = str2num("0" + StringFromList(4, info, ","))
                if (strlen(colortable) > 0)
                    ctab = colortable
                endif
                check_contrast(w, pcmin, pcmax, vmin, vmax)
                ModifyGraph /w=$graphname zColor($objname)={w, vmin, vmax, $ctab, rev}
            endif
        endif
    endfor
    
    string images = ImageNameList(graphname, ";")
    n = ItemsInList(images, ";")
    for (i = 0; i < n; i += 1)
        objname = StringFromList(i, images, ";")
        wave w = ImageNameToWaveRef(graphname, objname)
        info = ImageInfo(graphname, objname, 0)
        if (strlen(info) > 0)
            info = StringByKey("RECREATION", info, ":", ";")
            info = StringByKey("ctab", info, "=", ";")
            if (strlen(info) > 2)
                info = info[1,strlen(info)-2]
                ctab = StringFromList(2, info, ",")
                rev = str2num("0" + StringFromList(3, info, ","))
                if (strlen(colortable) > 0)
                    ctab = colortable
                endif
                check_contrast(w, pcmin, pcmax, vmin, vmax)
                ModifyImage /w=$graphname $objname ctab={vmin, vmax, $ctab, rev}
            endif
        endif
    endfor

    setdatafolder save_df
};