#%TITLE% BLCAL.MAC
#%NAME% %B%blcal.mac%B% - do a calibration of the beamline flux.%BR%
#%CATEGORY% MX
#%END%

global FLUX_CALCSERVER
global FLUX_EGYMOT
global FLUX_SPECEH
global FLUX_TRANSMISSION
global FLUX_DIODETYPE
global FLUX_CALCULATION_FLAG
global FLUX_T[]
global FLUX_D[]
global FLUX_ELEMENT
global CALIBRATED_DIODE[]
global FLUX_CT_I0
global FLUX_CT_I0FLUX
global FLUX_CT_I1
global FLUX_CT_I1FLUX


#%UU% 
#%MDESC%
def flux_setup '{ 
  global FLUX_CALCSERVER
  global FLUX_EGYMOT
  global FLUX_SPECEH
  global FLUX_TRANSMISSION
  global FLUX_ELEMENT
  global FLUX_DIODETYPE

  FLUX_CALCSERVER="$1"
  #"id23/calc/1"
  FLUX_EGYMOT="$2"
  #"energy"
  FLUX_SPECEH="$3"
  #"lid232:eh1"
  FLUX_TRANSMISSION="$4"
  #"100 80 60 40 20"
  FLUX_ELEMENT="$5"
  #"Mn Fe Co Ni Cu Zn As Se Br Kr" 
  FLUX_DIODETYPE="$6"
  if (FLUX_DIODETYPE == 0)
    FLUX_DIODETYPE=3
}'

#%UU%
#%MDESC%
def flux_counters_setup '{
global FLUX_EGYMOT
global FLUX_CT_I0
global FLUX_CT_I0FLUX
global FLUX_CT_I1
global FLUX_CT_I1FLUX
global FLUX_COUNTERS

  FLUX_EGYMOT="$1"
  FLUX_CT_I0="$2"
  FLUX_CT_I0FLUX="$3"
  FLUX_CT_I1="$4"
  FLUX_CT_I1FLUX="$5"

  if (whatis("local_get_calibration") == 0)
    eval("def local_get_calibration() \'{ }\'")

  cdef("user_getcounts", "flux_counters_calculation()\n", "__flux__") 
}'

#%IU% (diode_type)
#%MDESC% Get the flux from calculation of a calibrated diode. The
#%B%diode_type%B% is 3 (AL diode) or 4 (Mylar diode).
def flux_from_diode(diode_type) '{
  local cmd
  
  if (diode_type==0) {
    diode_type=FLUX_DIODETYPE
  }

  getangles
  wm energy

  ct

  cmd = sprintf("flux_from_diode.start(%d, %f, 3, %g)", diode_type, A[motor_num(FLUX_EGYMOT)], S[diode])
  printf("\nCalculate diode flux using parameters: \n\tdiode type %d\n\tenergy %6.2f\n\tdiode thickness 3 (given as a constant)\n\tCALIBRATED DIODE reads %6.6g \n", diode_type, A[motor_num(FLUX_EGYMOT)], S[diode]) 
  return esrf_io(FLUX_CALCSERVER, "DevRemoteDouble", cmd)
}'

#%IU%
#%MDESC% Open the safety and fast shutters.
def flux_openshutters '{
  shopen
  remote_eval(FLUX_SPECEH, "msopen")
  remote_async(FLUX_SPECEH, "backout")
  wait(0x8)
}'

#%IU%
#%MDESC% Close the safety and fast shutters
def flux_closeshutters '{
  remote_eval(FLUX_SPECEH, "msclose")
  shclose
}'

#%IU% diode_type
#%MDESC% Get the flux procedure. The diode type is 3 (AL) or 4 (Mylar).
def get_flux '{
  # input parameter : diode type
  flux_openshutters

  print flux_from_diode($1)

  flux_closeshutters
}'

#%IU% (t)
#%MDESC% Set the transmission %B%t%B%.
def flux_set_transmission(t) '{
  print "setting transmission to " t
  remote_eval(FLUX_SPECEH, sprintf("transmission %d", t))
}'

#%IU% diode_type
#%MDESC% Get the flux for different transmissions (as defined in flux_setup).
#The %B%diode_type%B% is 3 (AL diode) or 4 (Mylar diode).
def flux_scan_transmission '{
local transmissions[]
local f
local t
local result
local r
local i0_value
local i1_value
local data_str

unglobal FLUX_T
unglobal FLUX_D
global FLUX_T[]
global FLUX_D[]
unglobal FLUX
global FLUX[]

  if  (($1 != 3) && ($1 != 4)) {
     eprintf ("diode type not specified, exiting...\n")
     exit
  }

  cdef("cleanup_once", "flux_closeshutters;", "__flux__")

  flux_openshutters

  split(FLUX_TRANSMISSION, transmissions)

  for (i in transmissions) {
    t = transmissions[i] 
    flux_set_transmission(t)

    f=flux_from_diode($1)
    i0_value = S[i0]	  
    i1_value = S[i1]

    FLUX_T["i0"][t]=f
    FLUX_T["i1"][t]=f
    FLUX_D["i0"][t]=i0_value
    FLUX_D["i1"][t]=i1_value
  }

  flux_closeshutters

  data_str="["
  ptctr=0
  for (t in FLUX_T["i0"]) {
    data_str=sprintf("%s(%f,%f),", data_str, FLUX_D["i0"][t], FLUX_T["i0"][t])
    printf ("Data point %d: i0 = %6.1f, diode = %6.6g\n",ptctr,FLUX_D["i0"][t],FLUX_T["i0"][t])
    ptctr++
  }
  data_str=sprintf("%s]", data_str)
  result = esrf_io(FLUX_CALCSERVER, "DevRemoteString", sprintf("flux_from_diode.leastsq(\"%s\")", data_str))
  split(result, r)

  FLUX["i0"]["value"]=r[0]
  FLUX["i0"]["err"]=r[1]
  printf("Calibration factor for i0 is %6.6g with error value %f.\n\n",FLUX["i0"]["value"]=r[0],  FLUX["i0"]["err"]=r[1])
  
  data_str="["
  ptctr=0
  for (t in FLUX_T["i1"]) {
    data_str=sprintf("%s(%f,%f),", data_str, FLUX_D["i1"][t], FLUX_T["i1"][t])
    printf ("Data point %d: i1 = %6.1f, diode = %6.6g\n",ptctr,FLUX_D["i1"][t],FLUX_T["i1"][t])
    ptctr++
  }
  data_str=sprintf("%s]", data_str)
  print data_str
  result = esrf_io(FLUX_CALCSERVER, "DevRemoteString", sprintf("flux_from_diode.leastsq(\"%s\")", data_str))
  split(result, r)

  FLUX["i1"]["value"]=r[0]
  FLUX["i1"]["err"]=r[1] 

  printf("Calibration factor for i1 is %6.6g with error value %f.\n",FLUX["i1"]["value"]=r[0],  FLUX["i1"]["err"]=r[1])

}'

#%IU%
#%MDESC% do the calibration for different energies (as defined by flus_setup).
def flux_scan_energy '{
  local elements[]
  local el
  local edge 
  local id
  unglobal CALIBRATED_DIODE
  global CALIBRATED_DIODE[]
  
  split(FLUX_ELEMENT, elements)

  open("/users/blissadm/local/spec/userconf/calibrated_diodes.dat")

  for (i in elements) {
    el = elements[i]
    if (isnum(el) == 0) {
      # get edge energy
      printf ("Element: %s\n", el)
      if (remote_eval(FLUX_SPECEH, sprintf("_ae_getedge(\"%s\", \"K\")", el))==-1) {
        if (remote_eval(FLUX_SPECEH, sprintf("_ae_getedge(\"%s\", \"L\")", el))==-1) {
          eprintf ("Invalid element, cannot find edge %s, please remove from setup\n", el)
          exit  
        }
      }
      edge = prop_get(FLUX_SPECEH, "var/AE_E")
      printf("Edge: %g eV\n", edge)
    } else if (isnum(el) > 0) {
      if (el>500) {
        edge = el
      } else {
        edge = el * 1000
      }
      printf("Energy: %g eV\n", edge)
    } else {
      eprintf ("Invalid energy %g, please remove from setup\n", el)
      exit  
    }

    # have to do mono optimisation for energy
    id=remote_async(FLUX_SPECEH, sprintf("prodc_set_mono %g", edge/1000))
    wait(0x8)

    flux_scan_transmission $1

    CALIBRATED_DIODE[i]["energy"]=edge/1000
    CALIBRATED_DIODE[i]["i0"]=FLUX["i0"]["value"]
    CALIBRATED_DIODE[i]["i1"]=FLUX["i1"]["value"]
    on("/users/blissadm/local/spec/userconf/calibrated_diodes.dat")
    printf("%g %g %g\n", edge, FLUX["i0"]["value"], FLUX["i1"]["value"])
    off("/users/blissadm/local/spec/userconf/calibrated_diodes.dat")
  }

  print CALIBRATED_DIODE
  print "****************************/n/n"
}'

#%UU%
#%MDESC% The calibration procedure.
def flux_do_calibration '{
  close("/users/blissadm/local/spec/userconf/calibrated_diodes.dat")
  unix("rm -f /users/blissadm/local/spec/userconf/calibrated_diodes.dat")

  if (FLUX_ELEMENT=="") {
    # single wavelength beamline
    flux_scan_transmission $1
    getangles
    CALIBRATED_DIODE[i]["energy"]=A[motor_num(FLUX_EGYMOT)]*1000
    CALIBRATED_DIODE[i]["i0"]=FLUX["i0"]["value"]
    CALIBRATED_DIODE[i]["i1"]=FLUX["i1"]["value"]
    on("/users/blissadm/local/spec/userconf/calibrated_diodes.dat")
    printf("%g %g %g\n", CALIBRATED_DIODE[i]["energy"], FLUX["i0"]["value"], FLUX["i1"]["value"])
    off("/users/blissadm/local/spec/userconf/calibrated_diodes.dat")
  } else {
    flux_scan_energy $1
  }

  close("/users/blissadm/local/spec/userconf/calibrated_diodes.dat")
}'

#%IU%
#%MDESC% Load the calibration values from a file.
def flux_load_calibration '{
  unglobal CALIBRATED_DIODE
  global CALIBRATED_DIODE[]
  global CALIBRATED_TABLE
  local line
  local values[]
  local n i

  i = 0
  CALIBRATED_TABLE["ok"] = 0
  if (getline("/users/blissadm/local/spec/userconf/calibrated_diodes.dat", "open")!=0)
  {
    print "cannot open calibrated diodes file"
    exit
  }

  line = getline("/users/blissadm/local/spec/userconf/calibrated_diodes.dat") 

  while (line != -1) {
    n = split(line, values)
    CALIBRATED_DIODE[i]["energy"]=values[0]
    CALIBRATED_DIODE[i]["i0"]=values[1]
    CALIBRATED_DIODE[i]["i1"]=values[2]+0.0        
    i = i + 1
    line = getline("/users/blissadm/local/spec/userconf/calibrated_diodes.dat")
  }

  getline("/users/blissadm/local/spec/userconf/calibrated_diodes.dat", "close")
  CALIBRATED_TABLE["ok"] = 1
}'

#%IU% ()
#%MDESC% return the calibration values for the i0 and i1 counters at the
#current energy.
def flux_get_calibration() '{
global CALIBRATED_TABLE APERTURE_COEF
local i size egy cal

  APERTURE_COEF = 1

  if (CALIBRATED_TABLE["ok"] == 0) {
    printf ("loading the calibration table\n")
    flux_load_calibration
  }
  getangles
  egy = A[motor_num(FLUX_EGYMOT)]
  #the energy should be calculated in eV
  egy *= 1000

  size = asso_len(CALIBRATED_DIODE)
  local float array tmparr[size]

  if (size > 1) {
    for (i=0; i<size; i++)
      tmparr[i] = CALIBRATED_DIODE[i]["energy"]
    idx = _findidx(egy,tmparr,size-1)
    if (idx == -1) {
      printf ("could not find an index for energy %g: %g\n", egy, tmparr[i-1])
      return[0:-1, 1:-1]
    }
    cal_i0 = _interpol(CALIBRATED_DIODE, egy, idx, "i0")
    cal_i1 = _interpol(CALIBRATED_DIODE, egy, idx, "i1")
  } else {
    cal_i0 = CALIBRATED_DIODE["0"]["i0"]
    cal_i1 = CALIBRATED_DIODE["0"]["i1"]
  }
  cal = local_get_calibration()

  if (cal != 0) {
    APERTURE_COEF = cal
    cal_i0 *= cal
    cal_i1 *= cal
  }
  return [0:cal_i0, 1:cal_i1]
}'

#%IU% ()
#%MDESC% Calculate the flux and update the flux pseudo counters.
def flux_counters_calculation() '{
local cal[]

  if (FLUX_COUNTERS == 0)
    return(0)
  FLUX_CALCULATION_FLAG = FLUX_CALCULATION_FLAG+1
  cal = flux_get_calibration()
  S[cnt_num(FLUX_CT_I0FLUX)]=cal[0]*S[cnt_num(FLUX_CT_I0)]
  S[cnt_num(FLUX_CT_I1FLUX)]=cal[1]*S[cnt_num(FLUX_CT_I1)]
  if (FLUX_DEBUG) {
    printf("-- flux on I0 = %g (I0 value) x %g (calibration) = %g\n", S[cnt_num(FLUX_CT_I0)], cal[0], S[cnt_num(FLUX_CT_I0FLUX)])
    printf("-- flux on I1 = %g (I1 value) x %g (calibration) = %g\n", S[cnt_num(FLUX_CT_I1)], cal[1], S[cnt_num(FLUX_CT_I1FLUX)])
  }
}'

#%UU%
#%MDESC% Deactivate the flux pseudo counters.
def flux_counters_on '{
global FLUX_COUNTERS
  FLUX_COUNTERS = 1
  cdef("user_getcounts", "flux_counters_calculation()\n", "__flux__")
}'

#%UU%
#%MDESC% Activate the flux pseudo counters.
def flux_counters_off '{
  FLUX_COUNTERS = 0
  cdef("user_getcounts", "", "__flux__", "delete") 
}'

#%IU% (arr, val)
#%MDESC% Do linear interpolation between two points of the array %B%arr%B%
#for the x value %B%val%B%. Return the y value, or -1 if value out of bounds.
def _interpol(arr, val, idx, name) '{
local a b x1 x2 y1 y2 size dif

  size = asso_len(arr)

  x1 = arr[idx]["energy"]
  y1 = arr[idx][name]

  if (fabs(x1 - val) < 0.001)
    return(y1)

  if (x1 < val) {
    if (idx == 0)
      return(-1)
    x2 = arr[idx-1]["energy"]
    y2 = arr[idx-1][name]
  } else {
    if (idx == size)
      return(-1)
    x2 = arr[idx+1]["energy"]
    y2 = arr[idx+1][name]
  }
  dif = (x2-x1)
  if (dif == 0)
    dif = 1
  b = (y2-y1)/dif
  a = y1 - b*x1

  return(a+ b*val)
}'

#%IU% (val, arr, size)
#%MDESC% return the index of the closest to a value %B%val%B% element in the
#array %B%arr%B% of size %B%size%B%.
def _findidx(val,arr,size) '{
local i pos
local float array tmparr[size+1]

  pos = -1

  if (arr[0] > arr[size]) {
    if ((val > arr[0]) || (val < arr[size]))
      return(pos)
  } else {
    print "val " val " arr[0] " arr[0] " arr[size] " arr[size]
    if ((val < arr[0]) || (val > arr[size]))
      return(pos)
  }

  for (i=0; i<size+1; i++) {
    tmparr[i] = fabs(arr[i] - val)
    if (tmparr[i] < 0.0001)
      pos = i
  }
    
  if (pos == -1)
    pos = array_op("i_at_min", tmparr)
	
  return(pos)
}'

#%UU% (variable)
#%MDESC% Return 1 if %B%variable%B% is a number, 0 otherwise.
def isnum(variab) '{
  if ((variab+0) == (variab)) {
    return(1)
  } else {
    return(0)
  }
}'

#%MACROS%
#%IMACROS%
#%TOC%
#%AUTHOR% BLISS%BR%
#$Revision: 1.14 $$Date: 2011/04/26 14:45:20 $
#%END%
#%LOG%
#$Log: blcal.mac,v $
#Revision 1.14  2011/04/26 14:45:20  beteva
#flux_scan_energy accepts energies in eV or KeV as well.
#
#Revision 1.13  2011/01/26 16:42:39  beteva
#corrected a typo; added passing diode type when flux_do_calibration executed
#
#Revision 1.12  2010/12/14 14:32:10  guijarro
#added FLUX_CALCULATION_FLAG global
#
#Revision 1.11  2010/11/02 15:49:42  beteva
#changed flux_counters_setup due to changes in SPEC
#
#Revision 1.10  2010/10/29 11:41:17  beteva
#added local_get_calibration, FLUX_COUNTERS to check if counters on/off
#defined CALIBRATED_TABLE["ok"] to prevent loading the table at every count
#
#Revision 1.9  2009/09/08 16:01:22  spruce
#if there is only one energy, no need to interpolate.
#This caused problems when the energy value has a decimal place resulting in a -1 result
#
#Revision 1.8  2009/09/07 09:39:17  spruce
#Improvements to the messages printed during diode calibration to help the
#less softwarey people have an idea whats going on.
#
#Revision 1.7  2009/04/23 15:30:06  beteva
#changed interpolation routine
#