PBSwissMX/python/helicalscan.py

#!/usr/bin/env python
# *-----------------------------------------------------------------------*
# |                                                                       |
# |  Copyright (c) 2017 by Paul Scherrer Institute (http://www.psi.ch)    |
# |                                                                       |
# |              Author Thierry Zamofing (thierry.zamofing@psi.ch)        |
# *-----------------------------------------------------------------------*
'''
tools to setup and execute a helical scan of a cristal
#THIS IS JUST TESTING CODE TO SOLVE FINDING THE ROTATION CENTER
'''

import os, sys, json
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from utilities import *

class HelicalScan:
    def __init__(self,args):
      if args.cfg:
        fh=open(args.cfg,'r')
        s=fh.read()
        cfg=json.loads(s, object_hook=ConvUtf8)
        s=json.dumps(cfg, indent=2, separators=(',', ': '));print(s)
      else:
        fn='/tmp/shapepath4'
        #fn='/home/zamofing_t/Documents/prj/SwissFEL/epics_ioc_modules/ESB_MX/python/data/'+time.strftime('%y-%m-%d-%H_%M_%S')
        #cfg = {"sequencer": ['gen_grid_points(w=5,h=5,pitch=100,rnd=0.4)', 'sort_points()','gen_prog(file="'+fn+'.prg",host="SAR-CPPM-EXPMX1",mode=1,pt2pt_time=10,cnt=1)', 'plot_gather("'+fn+'.npz")']}
        #cfg = {"sequencer": ['test_find_rot_ctr()']}
        cfg = {"sequencer": ['test_find_rot_ctr(n=5. ,per=1. ,phi=24.6  ,bias=2.31,ampl=4.12)']}

      self.cfg=dotdict(cfg)
      self.args=args

    def run(self):
      print('args='+str(self.args))
      print('cfg='+str(self.cfg))
      #try:
      #  self.points=np.array(self.cfg.points)
      #except AttributeError:
      #  pass
      try:
        sequencer= self.cfg.pop('sequencer')
      except KeyError:
        print('no command sequence to execute')
      else:
        dryrun=self.args.dryrun
        for cmd in sequencer:
          print('>'*5+' '+cmd+' '+'<'*5)
          if not dryrun:
            eval('self.' + cmd)

    def meas_rot_ctr(self,y,per=1):
      # find the amplitude bias and phase of an equidistant sampled sinus
      # it needs at least 3 measurements e.g. at 0,120 240 deg or 0 90 180 270 deg
      # per is the number of persiods, default is 1 period =360 deg
      n=len(y)
      f = np.fft.fft(y)
      idx=int(per)
      bias=np.absolute(f[0]/n)
      phase=np.angle(f[idx])
      ampl=np.absolute(f[idx])*2/n
      return (bias,phase,ampl)

    def test_find_rot_ctr(self,n=3. ,per=1. ,phi=37  ,bias=4.1,ampl=2.4):
      # find the rotation center, amplitude out of n (niminum 3) measurements
      # n    number of equidistant measurements
      # per  number of periods (full rotation of all measurements nut be a interger value for precise measurements)
      # phi  phase
      # bias bias value
      # ampl amplitude

      t = np.arange(n)
      y=ampl*np.cos(2*np.pi*(per/n*t+phi/360.))+bias
      sp = np.fft.fft(y)
      freq = np.fft.fftfreq(t.shape[-1])
      plt.figure(1)
      plt.subplot(311)
      plt.plot(t,y,'b.-')
      plt.subplot(312)
      #plt.plot(t, sp.real,'b.-', t, sp.imag,'r.-')
      plt.step(t, sp.real,'b.-', t, sp.imag,'r.-', where='mid')
      #plt.stem(t, sp.real,'b-')
      #plt.plot(freq, sp.real,'b.-', freq, sp.imag,'r.-')

      idx=int(per)
      bias=np.absolute(sp[0]/n)
      phase=np.angle(sp[idx])
      ampl=np.absolute(sp[idx]) * 2 / n
      print('bias: '+str(bias))
      print('phase: '+str(phase*360./2/np.pi))
      print('amplitude: '+str(ampl))

      plt.subplot(313)
      t2 = np.linspace(0,2*np.pi,64)
      y2=ampl*np.cos(t2+phase)+bias
      plt.plot(t2,y2,'g-')
      plt.stem(t/n*2*np.pi*per,y,'b-')


      plt.show()
      pass


if __name__=='__main__':
  from optparse import OptionParser, IndentedHelpFormatter
  class MyFormatter(IndentedHelpFormatter):
    'helper class for formating the OptionParser'

    def __init__(self):
      IndentedHelpFormatter.__init__(self)

    def format_epilog(self, epilog):
      if epilog:
        return epilog
      else:
        return ""

  def parse_args():
    'main command line interpreter function'
    #usage: gpasciiCommunicator.py --host=PPMACZT84 myPowerBRICK.cfg
    (h, t)=os.path.split(sys.argv[0]);cmd='\n  '+(t if len(h)>3 else sys.argv[0])+' '
    exampleCmd=('-n',
                '-v15'
               )
    epilog=__doc__+'''
Examples:'''+''.join(map(lambda s:cmd+s, exampleCmd))+'\n '

    fmt=MyFormatter()
    parser=OptionParser(epilog=epilog, formatter=fmt)

    parser.add_option('-v', '--verbose', type="int", dest='verbose', help='verbosity bits (see below)', default=0)
    parser.add_option('-n', '--dryrun', action='store_true', help='dryrun to stdout')
    parser.add_option('--xy', action='store_true', help='sort x,y instead y,x')
    parser.add_option('--cfg', help='config file containing json configuration structure')

    (args, other)=parser.parse_args()
    args.other=other

    sp=HelicalScan(args)
    sp.run()
#------------------ Main Code ----------------------------------
  #ssh_test()
  ret=parse_args()
  exit(ret)