initial check-in

AU_Programs/2H/op_2h_exp.py

@@ -0,0 +1,635 @@
+# -*- coding: iso-8859-1 -*-
+
+TXEnableDelay = 2e-6
+TXEnableValue = 0b0001  # TTL line blanking RF amplifier (bit 0)
+TXPulseValue  = 0b0010  # TTL line triggering RF pulses  (bit 1)
+ADCSensitivity = 1	# voltage span for ADC
+
+def experiment(): # drives four experiments in a row: saturation-recovery, solid echo, spin-alignment, and Zeeman-order
+
+    # experiment toggles: 
+    satrec2_flag = True     # saturation-recovery on/off
+    solidecho_flag = True   # solid-echo on/off
+    spinal_flag = True      # spin-alignment on/off
+    zeeman_flag = False     # Zeeman-order on/off
+
+
+# ------------------ Saturation-recovery experiment settings ----------------------
+    
+    if satrec2_flag == True:  
+      # set up acquisition parameters:
+        pars = {}
+        pars['P90'] = 2.0e-6	# 90-degree pulse length (s)
+        pars['SF'] = 46.7e6	# spectrometer frequency (Hz)
+        pars['O1'] = -60e3	# offset from SF (Hz)
+        pars['SW'] = 10e6	# spectral window (Hz)
+        pars['SI'] = 1*512	# number of acquisition points
+        pars['NS'] = 8		# number of scans
+        pars['DS'] = 0		# number of dummy scans
+        pars['TAU'] = 1		# delay for recovery (s)
+        pars['D3'] = 20e-6	# echo delay (s)
+        pars['D4'] = 0		# echo pre-aquisition delay (s)
+        pars['PHA'] = 0		# receiver phase (degree)
+         # -*- these ain't variable: -*-
+        pars['NECH'] = 40       # number of saturation pulses
+        pars['D1'] = 100e-3	# starting interval in saturation sequence (s) 
+        pars['D2'] = 1e-4	# end interval in saturation sequence (s) 
+        pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+        pars['OUTFILE'] = None			# output file name
+
+      # specify a variable parameter (optional):
+        pars['VAR_PAR'] = 'TAU'	# variable parameter name (a string)   
+        start = 1e-3		# starting value
+        stop = 1		# end value
+        steps = 12		# number of values
+        log_scale = True	# log scale flag 
+        stag_range = False	# staggered range flag
+
+      # check parameters for safety: 
+        if pars['PHA'] < 0:
+            pars['PHA'] = 360 + pars['PHA']
+
+        if pars['P90'] > 20e-6:
+            raise Exception("pulse too long!!!")
+
+        # check whether a variable parameter is named: 
+        var_key = pars.get('VAR_PAR')
+        if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+            raise Exception("Pulse too long!!!")
+
+        if pars['NS']%8 != 0:
+            pars['NS'] = int(round(pars['NS'] / 8) + 1) * 8 
+            print 'Number of scans changed to ',pars['NS'],' due to phase cycling'
+
+        if pars['D1'] < pars['D2']:
+            raise Exception("D1 must be greater than D2!")
+
+        sat_length = sum(log_range(pars['D1'],pars['D2'],pars['NECH']))
+        if sat_length > 1.:
+            raise Exception("saturation sequence too long!!!")   
+
+        # start the experiment:
+        if var_key: 
+            # this is an arrayed experiment:
+            if log_scale:
+                array = log_range(start, stop, steps)
+            else:
+                array = lin_range(start, stop, steps)  
+
+            if stag_range:
+                array = staggered_range(array, size = 2)
+
+            # estimate the experiment time:
+            if var_key == 'TAU':
+                seconds = ((sat_length + pars['D3']*2) * steps + sum(array)) * (pars['NS'] + pars['DS'])
+            elif var_key == 'D3':
+                seconds = ((sat_length + pars['TAU']) * steps + sum(array)*2) * (pars['NS'] + pars['DS'])
+            else:
+                seconds = (sat_length + pars['TAU'] + pars['D3']*2) * steps * (pars['NS']+ pars['DS']) 
+            m, s = divmod(seconds, 60)
+            h, m = divmod(m, 60)
+            print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+            # loop for a variable parameter:
+            for index, pars[var_key] in enumerate(array):
+                print 'Arrayed saturation-recovery experiment for '+var_key+': run = '+str(index+1)+\
+                                 ' out of '+str(array.size)+': value = '+str(pars[var_key])
+                # loop for accumulation:
+                for run in xrange(pars['NS']+pars['DS']):
+                    yield satrec2_experiment(pars, run) 
+                synchronize()
+        else:
+           # estimate the experiment time:
+            seconds = (sat_length + pars['TAU'] + pars['D3']*2) * (pars['NS']+ pars['DS']) 
+            m, s = divmod(seconds, 60)
+            h, m = divmod(m, 60)
+            print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield satrec2_experiment(pars, run)
+
+# ------------------ Solid-echo experiment settings ----------------------
+    
+    if solidecho_flag == True: 
+      # set up acquisition parameters:
+        pars = {}
+        pars['P90'] = 2.0e-6	# 90-degree pulse length (s)
+        pars['SF'] = 46.7e6	# spectrometer frequency (Hz)
+        pars['O1'] = -60e3 	# offset from SF (Hz)
+        pars['SW'] = 500e3	# spectral window (Hz)
+        pars['SI'] = 1*1024	# number of acquisition points
+        pars['NS'] = 8		# number of scans
+        pars['DS'] = 0		# number of dummy scans
+        pars['RD'] = 1		# delay between scans (s)
+        pars['TAU'] = 20e-6	# echo delay (s)
+        pars['D4'] = 0e-6	# echo pre-acquisition delay (s)
+        pars['PHA'] = 0 	# receiver phase (degree)
+        pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+        pars['OUTFILE'] = None                  	# output file name
+
+      # specify a variable parameter (optional):
+        pars['VAR_PAR'] = None	# variable parameter name (a string)   
+        start = 10e-6		# starting value
+        stop = 30e-6		# end value
+        steps = 5		# number of values
+        log_scale = False	# log scale flag 
+        stag_range = False	# staggered range flag
+
+      # check parameters for safety: 
+        if pars['PHA'] < 0:
+            pars['PHA'] = 360 + pars['PHA']
+
+        if pars['P90'] > 20e-6:
+            raise Exception("Pulse too long!!!")
+
+        # check whether a variable parameter is named: 
+        var_key = pars.get('VAR_PAR')
+        if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+            raise Exception("Pulse too long!!!")
+
+        if pars['NS']%8 != 0:
+            pars['NS'] = int(round(pars['NS'] / 8) + 1) * 8 
+            print 'Number of scans changed to ',pars['NS'],' due to phase cycling' 
+
+      # start the experiment:
+        if var_key: 
+            # this is an arrayed experiment:
+            if log_scale:
+                array = log_range(start,stop,steps)
+            else:
+                array = lin_range(start,stop,steps)
+
+            if stag_range:
+                array = staggered_range(array, size = 2)
+
+            # estimate the experiment time:
+            if var_key == 'TAU':
+                seconds = (sum(array)*2 + pars['RD'] * steps) * (pars['NS'] + pars['DS'])
+            elif var_key == 'RD':
+                seconds = (sum(array) + pars['TAU']*2 * steps) * (pars['NS'] + pars['DS'])
+            else:
+                seconds = (pars['TAU']*2 + pars['RD']) * steps * (pars['NS']+ pars['DS'])
+            m, s = divmod(seconds, 60)
+            h, m = divmod(m, 60)
+            print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+            # loop for a variable parameter:
+            for index, pars[var_key] in enumerate(array):
+                print 'Arrayed solid-echo experiment for '+var_key+': run = '+str(index+1)+\
+                                 ' out of '+str(array.size)+': value = '+str(pars[var_key])
+                # loop for accumulation:
+                for run in xrange(pars['NS']+pars['DS']):
+                    yield solidecho_experiment(pars, run)
+                synchronize()
+        else:
+           # estimate the experiment time:
+            seconds = (pars['TAU']*2 + pars['RD']) * (pars['NS']+ pars['DS'])
+            m, s = divmod(seconds, 60)
+            h, m = divmod(m, 60)
+            print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield solidecho_experiment(pars, run)       
+
+
+# ---------------- Spin-alignment experiment settings ------------------
+
+    if spinal_flag == True:   
+      # set up acquisition parameters:         
+        pars = {}
+        pars['P90'] = 2.0e-6	# 90-degree pulse length (s)
+        pars['SF'] = 46.7e6	# spectrometer frequency (Hz)
+        pars['O1'] = -60e3 	# offset from SF (Hz)
+        pars['SW'] = 10e6    	# spectral window (Hz)
+        pars['SI'] = 1*512  	# number of acquisition points
+        pars['NS'] = 8		# number of scans
+        pars['DS'] = 0		# number of dummy scans
+        pars['RD'] = 1		# delay between scans (s)
+        pars['D1'] = 30e-6	# delay after first pulse, or tp (s)
+        pars['D2'] = 100e-6	# delay after second pulse, or tm (s)
+        pars['PHA'] = -36	# receiver phase (degree)
+        pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+        pars['OUTFILE'] = None			# output file name
+    
+      # specify a variable parameter (optional):
+        pars['VAR_PAR'] = 'D2'	# variable parameter name (a string)   
+        start = 30e-6		# starting value
+        stop = 1e-3		# end value
+        steps = 12		# number of values
+        log_scale = True	# log scale flag 
+        stag_range = False	# staggered range flag
+
+      # check parameters for safety:
+        if pars['PHA'] < 0:
+            pars['PHA'] = 360 + pars['PHA']
+        
+        if pars['P90'] > 20e-6:
+            raise Exception("Pulse too long!!!")
+                    
+        # check whether a variable parameter is named: 
+        var_key = pars.get('VAR_PAR')
+        if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+            raise Exception("Pulse too long!!!")
+
+        if pars['NS']%8 != 0:
+            pars['NS'] = int(round(pars['NS'] / 8) + 1) * 8 
+            print 'Number of scans changed to ', pars['NS'], ' due to phase cycling'
+
+        # start the experiment:
+        if var_key: 
+            # this is an arrayed experiment:
+            if log_scale:
+                array = log_range(start,stop,steps)
+            else:
+                array = lin_range(start,stop,steps) 
+        
+            if stag_range:
+                array = staggered_range(array, size = 2)
+
+           # estimate the experiment time:
+            if var_key == 'D1':
+                seconds = (sum(array)*2 + (pars['D2'] + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+            elif var_key == 'D2':
+                seconds = (sum(array) + (pars['D1']*2 + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+            elif var_key == 'RD':
+                seconds = (sum(array) + (pars['D1']*2 + pars['D2']) * steps) * (pars['NS'] + pars['DS'])
+            else:
+                seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * steps * (pars['NS']+ pars['DS'])
+            m, s = divmod(seconds, 60)
+            h, m = divmod(m, 60)
+            print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+            # loop for a variable parameter:
+            for index, pars[var_key] in enumerate(array):
+                print 'Arrayed spin-alignment experiment for '+var_key+': run = '+str(index+1)+\
+                                 ' out of '+str(array.size)+': value = '+str(pars[var_key])
+                # loop for accumulation:
+                for run in xrange(pars['NS']+pars['DS']):
+                    yield spinal_experiment(pars, run)
+                synchronize()
+
+        else:
+           # estimate the experiment time:
+            seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * (pars['NS']+ pars['DS'])
+            m, s = divmod(seconds, 60)
+            h, m = divmod(m, 60)
+            print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+           # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield spinal_experiment(pars, run)
+
+# ------------------ Zeeman-order experiment settings ----------------------
+
+    if zeeman_flag == True:
+        # set up acquisition parameters:         
+        pars = {}
+        pars['P90'] = 2.0e-6	# 90-degree pulse length (s)
+        pars['SF'] = 46.7e6 	# spectrometer frequency (Hz)
+        pars['O1'] = -60e3	# offset from SF (Hz)
+        pars['SW'] = 10e6	# spectral window (Hz)
+        pars['SI'] = 1*512	# number of acquisition points
+        pars['NS'] = 8		# number of scans
+        pars['DS'] = 0		# number of dummy scans
+        pars['RD'] = 1		# delay between scans (s)
+        pars['D1'] = 30e-6	# delay after first pulse, or 'short tau' (s)
+        pars['D2'] = 100e-6	# delay after second pulse, or 'long tau' (s)
+        pars['PHA'] = 0 	# receiver phase (degree)
+        pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+        pars['OUTFILE'] = None                  	# output file name
+
+        # specify a variable parameter (optional):
+        pars['VAR_PAR'] = 'D2'	# variable parameter name (a string)
+        start = 30e-6		# starting value
+        stop = 1e-3		# end value
+        steps = 12		# number of values
+        log_scale = True	# log scale flag 
+        stag_range = False	# staggered range flag
+
+        # check parameters for safety: 
+        if pars['PHA'] < 0:
+            pars['PHA'] = 360 + pars['PHA']
+
+        if pars['P90'] > 20e-6:
+            raise Exception("Pulse too long!!!")
+
+        # check whether a variable parameter is named: 
+        var_key = pars.get('VAR_PAR')
+        if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+            raise Exception("Pulse too long!!!")
+
+        if pars['NS']%8 != 0:
+            pars['NS'] = int(round(pars['NS'] / 8) + 1) * 8 
+            print 'Number of scans changed to ',pars['NS'],' due to phase cycling'
+
+        # start the experiment:
+        if var_key:
+            # this is an arrayed experiment:
+            if log_scale:
+                array = log_range(start,stop,steps)
+            else:
+                array = lin_range(start,stop,steps)
+
+            if stag_range:
+                array = staggered_range(array, size = 2)
+
+            # estimate the experiment time:
+            if var_key == 'D1':
+                seconds = (sum(array)*2 + (pars['D2'] + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+            elif var_key == 'D2':
+                seconds = (sum(array) + (pars['D1']*2 + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+            elif var_key == 'RD':
+                seconds = (sum(array) + (pars['D1']*2 + pars['D2']) * steps) * (pars['NS'] + pars['DS'])
+            else:
+                seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * steps * (pars['NS']+ pars['DS'])
+            m, s = divmod(seconds, 60)
+            h, m = divmod(m, 60)
+            print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+            # loop for a variable parameter:
+            for index, pars[var_key] in enumerate(array):
+                print 'Arrayed Zeeman-order experiment for '+var_key+': run = '+str(index+1)+\
+                                 ' out of '+str(array.size)+': value = '+str(pars[var_key])
+                # loop for accumulation:
+                for run in xrange(pars['NS']+pars['DS']):
+                    yield zeeman_experiment(pars, run) 
+                synchronize()
+        else:
+            # estimate the experiment time:
+            seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * (pars['NS']+ pars['DS'])
+            m, s = divmod(seconds, 60)
+            h, m = divmod(m, 60)
+            print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield zeeman_experiment(pars, run)
+			
+# the pulse programs:
+
+def satrec2_experiment(pars, run):
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+    
+    pars['PROG'] = 'satrec2_experiment'
+    
+    # phase lists:
+    pars['PH1'] = [  0] 					# saturation pulses
+    pars['PH3'] = [ 0, 180,   0, 180, 90, 270,  90, 270]	# 1st 90-degree pulse
+    pars['PH4'] = [90,  90, 270, 270,  0,   0, 180, 180]	# 2nd 90-degree pulse
+    pars['PH2'] = [ 0, 180,   0, 180, 90, 270,  90, 270]	# receiver
+
+    # read in variables:
+    P90 = pars['P90']
+    SF =  pars['SF']
+    O1 =  pars['O1']
+    NECH = pars['NECH']
+    D1 =  pars['D1']
+    D2 =  pars['D2']
+    D3 =  pars['D3']
+    D4 =  pars['D4']
+    TAU = pars['TAU']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH4 = pars['PH4'][run%len(pars['PH4'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']
+
+    # set variable delay list for saturation pulses: 
+    vdlist = log_range(D2, D1, NECH-1)
+
+    # set sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+
+    # the pulse sequence: 
+
+    # saturation:
+    e.set_frequency(SF+O1, phase=PH1)			# set frequency and phase for saturation pulses
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# apply 90-degree pulse
+    for delay in vdlist[::-1]:
+        e.wait(delay-P90-TXEnableDelay)			# wait for next saturation pulse
+        e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+        e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue) # apply 90-degree pulse
+
+    # recovery:
+    e.wait(TAU)						# wait for tau
+    e.set_phase(PH3)					# set phase for next pulse
+     
+    # echo detection:
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# apply 90-degree pulse
+    e.wait(D3-P90/2-TXEnableDelay)			# echo delay
+    e.set_phase(PH4)					# set phase for next pulse
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# apply 90-degree pulse 
+    e.set_phase(PHA)					# set phase for receiver
+    e.wait(D3-P90/2+D4)					# echo delay    
+    e.record(SI, SW, sensitivity=ADCSensitivity)	# acquisition
+
+    # write experiment attributes:  
+    for key in pars.keys():
+        e.set_description(key, pars[key])	# pulse sequence parameters
+    e.set_description('run', run)		# current scan
+    e.set_description('rec_phase', -PH2)	# current receiver phase
+
+    return e    
+        
+
+def solidecho_experiment(pars, run): 
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+
+    pars['PROG'] = 'solidecho_experiment'
+
+    # phase lists [from Tecmag's pulse sequence]:
+    pars['PH1'] = [ 0, 180,   0, 180, 90, 270,  90, 270]	# 1st 90-degree pulse
+    pars['PH3'] = [90,  90, 270, 270,  0,   0, 180, 180]	# 2nd 90-degree pulse
+    pars['PH2'] = [ 0, 180,   0, 180, 90, 270,  90, 270]	# receiver
+
+
+    # read in variables:
+    P90 = pars['P90']
+    SF = pars['SF']
+    O1 = pars['O1'] 
+    RD =  pars['RD']
+    TAU = pars['TAU']
+    D4 = pars['D4']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']
+
+    # set sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+
+    # run the pulse sequence: 
+    e.wait(RD)						# delay between scans
+    e.set_frequency(SF+O1, phase=PH1)			# set frequency and phase for 1st RF pulse
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# apply 1st 90-degree pulse
+    e.wait(TAU-P90/2-TXEnableDelay)			# wait for TAU
+    e.set_phase(PH3)					# set phase for 2nd 90-degree pulse
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enalble RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# apply 2nd 90-degree pulse
+    e.set_phase(PHA)					# set phase for receiver
+    e.wait(TAU-P90/2+D4)				# wait for TAU       
+    e.record(SI, SW, sensitivity=ADCSensitivity)	# acquire echo points
+
+    # write the experiment parameters:  
+    for key in pars.keys():
+        e.set_description(key, pars[key])	# pulse sequence parameters
+    e.set_description('run', run)		# current scan
+    e.set_description('rec_phase', -PH2)	# current receiver phase
+
+    return e
+    
+def spinal_experiment(pars, run):
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+                       
+    pars['PROG'] = 'spinal_experiment'                       
+                       
+   # 8-step phase cycle (1-14 modifided to deal with T1-recovery and extended for Re/Im imbalance)
+    pars['PH1'] = [0, 270,   0, 270,     90,  90, 180, 180 ] # 1st (90-degree) pulse
+    pars['PH3'] = [90,180,  90, 180,    180, 180,  90,  90 ] # 2nd (45-degree) pulse
+    pars['PH4'] = [90, 90, 270, 270,    180,   0,   0, 180 ] # 3rd (45-degree) pulse
+    pars['PH2'] = [0, 180, 180,   0,     90, 270,  90, 270 ] # receiver
+        
+   # read in variables:
+    P90 = pars['P90']
+    P45 = pars['P90']*0.5
+    P1  = pars['P90']*0.5
+    SF =  pars['SF']
+    O1 =  pars['O1']    
+    RD =  pars['RD']
+    D1 =  pars['D1']
+    D2 =  pars['D2']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH4 = pars['PH4'][run%len(pars['PH4'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']
+
+   # set sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+
+   # run the pulse sequence:
+    e.wait(RD)						# relaxation delay between scans
+    e.set_frequency(SF+O1, phase=PH1)
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse 
+
+    e.wait(D1-P90/2-TXEnableDelay)			# 'short tau'
+    e.set_phase(PH3)
+
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)
+    e.ttl_pulse(P45, value=TXEnableValue|TXPulseValue)	# 45-degree pulse 
+
+    e.wait(D2-P45/2-TXEnableDelay)			# 'long tau'
+    e.set_phase(PH4)
+
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)
+    e.ttl_pulse(P1, value=TXEnableValue|TXPulseValue)	# 45-degree pulse
+
+    e.wait(TXEnableDelay)
+    e.set_phase(PHA)
+    e.wait(5e-6)#D1-P45/2-TXEnableDelay)    	                # 'short tau'
+    e.record(SI, SW, sensitivity=ADCSensitivity)        # acquisition                                     
+
+   # write experiment parameters:
+    for key in pars.keys():
+        e.set_description(key, pars[key])	# acquisition parameters
+    e.set_description('run', run)		# current scan
+    e.set_description('rec_phase', -PH2)	# current receiver phase
+
+    return e
+
+def zeeman_experiment(pars, run):
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+        
+    pars['PROG'] = 'zeeman_experiment'
+
+   # 8-step phase cycle (1-21 modifided to deal with T1-recovery and extended for Re/Im imbalance)
+    pars['PH1'] = [0, 270,   0, 270,   180,  90, 180, 90] # 1st (90-degree) pulse
+    pars['PH3'] = [0,  90,   0,  90,     0,  90,   0, 90] # 2nd (90-degree) pulse
+    pars['PH4'] = [0,   0, 180, 180,   270, 270,  90, 90]  # 3rd (90-degree) pulse
+    pars['PH2'] = [0, 180, 180,   0,    90, 270, 270, 90]  # receiver
+
+   # read in variables:
+    P90 = pars['P90']
+    SF =  pars['SF']
+    O1 =  pars['O1']    
+    RD =  pars['RD']
+    D1 =  pars['D1']
+    D2 =  pars['D2']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH4 = pars['PH4'][run%len(pars['PH4'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']
+
+   # set sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+
+   # run the pulse sequence:
+    e.wait(RD)                                  	# relaxation delay between scans
+    e.set_frequency(SF+O1, phase=PH1)
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse
+
+    e.wait(D1-P90/2-TXEnableDelay)			# 'short tau'
+
+    e.set_phase(PH3)
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse 
+
+    e.wait(D2-P90/2-TXEnableDelay)			# 'long tau'
+
+    e.set_phase(PH4)
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse
+
+    e.wait(TXEnableDelay)
+    e.set_phase(PHA)
+    e.wait(D1-P90/2-TXEnableDelay)			# 'short tau'
+    e.record(SI, SW, sensitivity=ADCSensitivity)	# acquisition
+
+   # write experiment parameters:
+    for key in pars.keys():
+        e.set_description(key, pars[key])	# acquisition parameters
+    e.set_description('run', run)		# current scan
+    e.set_description('rec_phase', -PH2)	# current receiver phase
+
+    return e

AU_Programs/2H/op_2h_res.py

@@ -0,0 +1,282 @@
+# -*- coding: iso-8859-1 -*-
+
+from numpy import *
+from scipy.signal import *
+from scipy.optimize import *
+from os import path, rename
+
+def result():
+
+    the_experiment = None	# current experiment's name
+    
+    measurements = {'satrec2_experiment': MeasurementResult('Saturation Recovery'),
+                  'solidecho_experiment': MeasurementResult('Solid Echo'),
+                     'spinal_experiment': MeasurementResult('Spin Alignment'),
+                     'zeeman_experiment': MeasurementResult('Zeeman Order')}
+                               
+    measurement_ranges = {'satrec2_experiment': [0.5e-6, 4.5e-6],
+                        'solidecho_experiment': [0.5e-6, 4.5e-6],
+                           'spinal_experiment': [0.5e-6, 4.5e-6],
+                           'zeeman_experiment': [0.5e-6, 4.5e-6]}
+    measurement_ranging = True
+                             
+    # loop over the incoming results: 
+    for timesignal in results:
+        if not isinstance(timesignal,ADC_Result): 
+            continue     
+
+        # read experiment parameters:
+        pars = timesignal.get_description_dictionary()
+                    
+        # keep track of the actual experiment's name:
+        if the_experiment != pars.get('PROG'):
+           the_experiment = pars.get('PROG')
+           suffix = ''	# output file name's suffix
+           counter = 1
+
+        # ---------------- digital filter ------------------ 
+             
+        # get actual sampling rate of timesignal:
+        sampling_rate = timesignal.get_sampling_rate()
+
+        # get user-defined spectrum width:
+        spec_width = pars['SW']
+        
+        # specify cutoff frequency, in relative units:
+        cutoff = spec_width / sampling_rate
+                    
+        # number of filter's coefficients:
+        numtaps = 29
+                
+        if cutoff < 1: # otherwise no filter applied
+            
+            # use firwin to create a lowpass FIR filter:
+            fir_coeff = firwin(numtaps, cutoff)
+            
+            # downsize x according to user-defined spectral window:
+            skip = int(sampling_rate / spec_width)
+            timesignal.x = timesignal.x[::skip]
+        
+            for i in range(2):
+                # apply the filter to ith channel:
+                timesignal.y[i] = lfilter(fir_coeff, 1.0, timesignal.y[i])
+
+                # zeroize first N-1 "corrupted" samples:
+                timesignal.y[i][:numtaps-1] = 0.0
+
+                # circular left shift of y:
+                timesignal.y[i] = roll(timesignal.y[i], -(numtaps-1))    
+
+                # downsize y to user-defined number of samples (SI):
+                timesignal.y[i] = timesignal.y[i][::skip]
+
+            # update the sampling_rate attribute of the signal's:
+            timesignal.set_sampling_rate(spec_width)
+
+        # ----------------------------------------------------
+
+        # rotate timesignal according to current receiver's phase:
+        timesignal.phase(pars['rec_phase'])
+
+        # provide timesignal to the display tab:
+        data['Current scan'] = timesignal
+
+        # accumulate...
+        if not locals().get('accu'):
+            accu = Accumulation()
+
+        # skip dummy scans, if any:
+        if pars['run'] < 0: continue
+
+        # add up:
+        accu += timesignal
+
+        # provide accumulation to the display tab:
+        data['Accumulation'] = accu
+
+        # check how many scans are done:
+        if accu.n == pars['NS']: # accumulation is complete
+           
+	    # make a copy:
+            echo = accu + 0
+
+            # compute the signal's phase:
+            #phi0 = arctan2(accu.y[1][0], accu.y[0][0]) * 180 / pi
+            #if not locals().get('ref'): ref = phi0   
+            #print 'phi0 = ', phi0
+
+            # rotate the signal to maximize Re (optional):
+            #echo.phase(-phi0)
+
+            # check whether it is an arrayed experiment: 
+            var_key = pars.get('VAR_PAR')
+            if var_key:
+                # get variable parameter's value:
+                var_value = pars.get(var_key)	
+
+                # provide signal recorded with this var_value to the display tab:
+                data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu
+
+		# estimate noise level:
+                if not locals().get('noise'):
+                    n = int(0.1*accu.x.size)
+                    noise = 3*std(accu.y[0][-n-29:-30])  				
+					
+                # measure signal intensity vs. var_value:	
+                if the_experiment in measurements.keys():
+				
+                    # option a: the sum of samples within the given range:
+                    if measurement_ranging == True:
+                        [start, stop] = echo.get_sampling_rate() * array(measurement_ranges[the_experiment])
+                        measurements[the_experiment][var_value] = sum(echo.y[0][int(start):int(stop)]) 
+						
+                    # option b: the sum of all samples above noise: 
+                    else:
+                        measurements[the_experiment][var_value] = sum(echo.y[0][echo.y[0]>noise])
+                    
+                    # store a measurement:
+                    data[measurements[the_experiment].get_title()] = measurements[the_experiment]
+                    
+                    # update the file name suffix:
+                    suffix = '_' + str(counter)
+                    counter += 1
+					
+                else: 
+                    print "Cannot recognize experiment: continue without measuring"				
+
+            # save accu if required:   
+            outfile = pars.get('OUTFILE')
+            if outfile: 
+                datadir = pars.get('DATADIR')
+
+                # write raw data in Simpson format:
+                filename = datadir+outfile+suffix+'.dat'
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.dat')
+                accu.write_to_simpson(filename)
+
+               # write parameters in a text file:
+                filename = datadir+outfile+suffix+'.par' 
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.par')
+
+                fileobject = open(filename, 'w')
+                for key in sorted(pars.iterkeys()):
+                    if key=='run': continue
+                    if key=='rec_phase': continue
+                    fileobject.write('%s%s%s'%(key,'=', pars[key]))
+                    fileobject.write('\n')
+                fileobject.close()
+
+            # reset accumulation:
+            del accu
+
+    if var_key == 'TAU' or var_key == 'D2':
+       # mono-exponential saturation-recovery fit:
+        try:
+            xdata = measurements['satrec2_experiment'].get_xdata()
+            ydata = measurements['satrec2_experiment'].get_ydata()
+            [amplitude, rate, offset] = fitfunc_recovery(xdata, ydata)
+            print 'Mono-exponential fit to saturation-recovery data:'
+            print '%s%02g' % ('Amplitude = ', amplitude)
+            print '%s%02g' % ('T1 [s] = ', 1./rate)
+
+            # update display for fit:      
+            measurements['satrec2_experiment'].y = func_recovery([amplitude, rate, offset], xdata)
+            data[measurements['satrec2_experiment'].get_title()] = measurements['satrec2_experiment']
+        except:
+            pass
+		
+        # KWW fit to spin-alignment echoes:
+        try:
+            xdata = measurements['spinal_experiment'].get_xdata()
+            ydata = measurements['spinal_experiment'].get_ydata()
+            [amplitude, rate, beta] = fitfunc_kww(xdata, ydata)
+            print 'KWW fit to spin-alignment echoes:'
+            print '%s%02g' % ('Amplitude = ', amplitude)
+            print '%s%02g' % ('T2 [s] = ', 1./rate)
+            print '%s%01g' % ('Beta = ', beta)
+
+            # update display for the fit:
+            measurements['spinal_experiment'].y = func_kww([amplitude, rate, beta], xdata)
+            data[measurements['spinal_experiment'].get_title()] = measurements['spinal_experiment']
+        except:
+            pass
+		
+        # KWW fit to Zeeman-order echoes:
+        try:
+            xdata = measurements['zeeman_experiment'].get_xdata()
+            ydata = measurements['zeeman_experiment'].get_ydata()
+            [amplitude, rate, beta] = fitfunc_kww(xdata, ydata)
+            print 'KWW fit to Zeeman-order echoes:'
+            print '%s%02g' % ('Amplitude = ', amplitude)
+            print '%s%02g' % ('T2 [s] = ', 1./rate)
+            print '%s%01g' % ('Beta = ', beta)
+
+            # update display for the fit:
+            measurements['zeeman_experiment'].y = func_kww([amplitude, rate, beta], xdata)
+            data[measurements['zeeman_experiment'].get_title()] = measurements['zeeman_experiment']
+        except:
+            pass	
+
+# the fitting procedure for satrec_experiment:
+def fitfunc_recovery(xdata, ydata):
+
+    # initialize variable parameters:
+    try: 
+        # solve Az = b:
+        A = array((ones(xdata.size/2), xdata[0:xdata.size/2]))    
+        b = log(abs(mean(ydata[-2:]) - ydata[0:xdata.size/2]))     
+        z = linalg.lstsq(transpose(A), b)
+        amplitude = exp(z[0][0])
+        rate = -z[0][1]
+    except:
+        amplitude = abs(ydata[-1] - ydata[0])
+        rate = 1./(xdata[-1] - xdata[0])
+    offset = min(ydata)
+    p0 = [amplitude, rate, offset]
+
+    # run least-squares optimization:
+    plsq = leastsq(residuals_recovery, p0, args=(xdata, ydata))
+
+    return plsq[0] # best-fit parameters
+
+def residuals_recovery(p, xdata, ydata):
+    return ydata - func_recovery(p, xdata)
+
+# here is the function to fit
+def func_recovery(p, xdata):
+    return p[0]*(1-exp(-p[1]*xdata)) + p[2]
+
+			
+# the fitting procedure for spinal_experiment and zeeman_experiment:
+def fitfunc_kww(xdata, ydata): 
+
+    # initialize variable parameters:
+    try: 
+        # solve Az = b:
+        A = array((ones(xdata.size/2), xdata[0:xdata.size/2]))    
+        b = log(abs(ydata[0:xdata.size/2]))     
+        z = linalg.lstsq(transpose(A), b)
+        amplitude = exp(z[0][0])
+        rate = -z[0][1]
+    except:
+        amplitude = abs(ydata[0])
+        rate = 1./(xdata[-1] - xdata[0])
+    beta = 1.
+    p0 = [amplitude, rate, beta]
+
+    # run least-squares optimization:
+    plsq = leastsq(residuals_kww, p0, args=(xdata, ydata))
+
+    return plsq[0] # best-fit parameters
+
+def residuals_kww(p, xdata, ydata):
+    return ydata - func_kww(p, xdata)
+
+# here is the function to fit:
+def func_kww(p, xdata):
+    return p[0]*exp(-(p[1]*xdata)**p[2])
+
+
+pass

AU_Programs/Diffusiometry/op_diff_exp.py

@@ -0,0 +1,483 @@
+# -*- coding: iso-8859-1 -*-
+
+TXEnableDelay = 2e-6
+TXEnableValue = 0b0001  # TTL line blanking RF amplifier (bit 0)
+TXPulseValue  = 0b0010  # TTL line triggering RF pulses  (bit 1)
+ADCSensitivity = 1	# voltage span for ADC
+
+def experiment(): # drives three experiments in a row: saturation-recovery, stimulated echo, and Hanh echo
+
+    # experiments switches: 
+    satrec2_flag = True     # saturation-recovery switch
+    ste_flag = True         # stimulated-echo switch
+    hahn_flag = False       # Hahn-echo switch
+
+
+# ------------------ Saturation-recovery experiment settings ----------------------
+    
+    if satrec2_flag == True:  
+      # set up acquisition parameters:
+        pars = {}
+        pars['P90'] = 1.7e-6	# 90-degree pulse length (s)
+        pars['SF'] = 62.92e6	# spectrometer frequency (Hz)
+        pars['O1'] = 0e3	# offset from SF (Hz)
+        pars['SW'] = 20e6	# spectral window (Hz)
+        pars['SI'] = 1*512	# number of acquisition points
+        pars['NS'] = 512	# number of scans
+        pars['DS'] = 0		# number of dummy scans
+        pars['TAU'] = 1		# delay for recovery (s)
+        pars['D3'] = 20e-6	# echo delay (s)
+        pars['D4'] = 0		# echo pre-aquisition delay (s)
+        pars['PHA'] = 50	# receiver phase (degree)
+         # -*- these ain't variable: -*-
+        pars['NECH'] = 20       # number of saturation pulses
+        pars['D1'] = 100e-3	# starting interval in saturation sequence (s) 
+        pars['D2'] = 1e-4	# end interval in saturation sequence (s) 
+        pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+        pars['OUTFILE'] = '360K'			# output file name
+
+      # specify a variable parameter (optional):
+        pars['VAR_PAR'] = 'TAU'	# variable parameter name (a string)   
+        start = 1e-3		# starting value
+        stop = 1e-0		# end value
+        steps = 12		# number of values
+        log_scale = True	# log scale flag 
+        stag_range = False	# staggered range flag
+
+      # check parameters for safety: 
+        if pars['PHA'] < 0:
+            pars['PHA'] = 360 + pars['PHA']
+
+        if pars['P90'] > 20e-6:
+            raise Exception("pulse too long!!!")
+
+        # check whether a variable parameter is named: 
+        var_key = pars.get('VAR_PAR')
+        if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+            raise Exception("Pulse too long!!!")
+
+        if pars['NS']%8 != 0:
+            pars['NS'] = int(round(pars['NS'] / 8) + 1) * 8 
+            print 'Number of scans changed to ',pars['NS'],' due to phase cycling'
+
+        if pars['D1'] < pars['D2']:
+            raise Exception("D1 must be greater than D2!")
+
+        sat_length = sum(log_range(pars['D1'],pars['D2'],pars['NECH']))
+        if sat_length > 1.:
+            raise Exception("saturation sequence too long!!!")   
+
+        # start the experiment:
+        if var_key: 
+            # this is an arrayed experiment:
+            if log_scale:
+                array = log_range(start, stop, steps)
+            else:
+                array = lin_range(start, stop, steps)  
+
+            if stag_range:
+                array = staggered_range(array, size = 2)
+
+            # estimate the experiment time:
+            if var_key == 'TAU':
+                seconds = ((sat_length + pars['D3']*2) * steps + sum(array)) * (pars['NS'] + pars['DS'])
+            elif var_key == 'D3':
+                seconds = ((sat_length + pars['TAU']) * steps + sum(array)*2) * (pars['NS'] + pars['DS'])
+            else:
+                seconds = (sat_length + pars['TAU'] + pars['D3']*2) * steps * (pars['NS']+ pars['DS']) 
+            m, s = divmod(seconds, 60)
+            h, m = divmod(m, 60)
+            print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+            # loop for a variable parameter:
+            for index, pars[var_key] in enumerate(array):
+                print 'Arrayed saturation-recovery experiment for '+var_key+': run = '+str(index+1)+\
+                                 ' out of '+str(array.size)+': value = '+str(pars[var_key])
+                # loop for accumulation:
+                for run in xrange(pars['NS']+pars['DS']):
+                    yield satrec2_experiment(pars, run) 
+                synchronize()
+        else:
+           # estimate the experiment time:
+            seconds = (sat_length + pars['TAU'] + pars['D3']*2) * (pars['NS']+ pars['DS']) 
+            m, s = divmod(seconds, 60)
+            h, m = divmod(m, 60)
+            print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield satrec2_experiment(pars, run)
+
+
+# ---------------- Stimulated-echo experiment settings ------------------
+
+    if ste_flag == True:   
+      # set up acquisition parameters:
+        pars = {}
+        pars['P90'] = 1.6e-6	# 90-degree pulse length (s)
+        pars['SF'] = 62.92e6	# spectrometer frequency (Hz)
+        pars['O1'] = 0e3	# offset from SF (Hz)
+        pars['SW'] = 20e6	# spectral window (Hz)
+        pars['SI'] = 1*512	# number of acquisition points
+        pars['NS'] = 512	# number of scans
+        pars['DS'] = 0		# number of dummy scans
+        pars['RD'] = 1		# delay between scans (s)
+        pars['D1'] = 400e-6	# delay after first pulse (short tau) (s)
+        pars['D2'] = 20e-6	# delay after second pulse (long tau) (s) 
+        pars['PHA'] = 240	# receiver phase (degree)    
+        pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+        pars['OUTFILE'] = '210K'  			# output file name
+
+      # specify a variable parameter (optional):
+        pars['VAR_PAR'] = 'D2'	# variable parameter name (a string)   
+        start = 350e-6		# starting value
+        stop = 1e-0			# end value
+        steps = 16			# number of values
+        log_scale = True		# log scale flag 
+        stag_range = False		# staggered range flag
+
+      # check parameters for safety: 
+        if pars['PHA'] < 0:
+            pars['PHA'] = 360 + pars['PHA']
+
+        if pars['P90'] > 20e-6:
+            raise Exception("Pulse too long!!!")
+
+        # check whether a variable parameter is named: 
+        var_key = pars.get('VAR_PAR')
+        if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+            raise Exception("Pulse too long!!!")
+
+        if pars['NS']%16 != 0:
+            pars['NS'] = int(round(pars['NS'] / 16) + 1) * 16 
+            print 'Number of scans changed to ',pars['NS'],' due to phase cycling'
+
+        # start the experiment:
+        if var_key: 
+            # this is an arrayed experiment:
+            if log_scale:
+                array = log_range(start,stop,steps)
+            else:
+                array = lin_range(start,stop,steps) 
+
+            if stag_range:
+                array = staggered_range(array, size = 2)
+
+            # estimate the experiment time:
+            if var_key == 'D1':
+                seconds = (sum(array)*2 + (pars['D2'] + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+            elif var_key == 'D2':
+                seconds = (sum(array) + (pars['D1']*2 + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+            elif var_key == 'RD':
+                seconds = (sum(array) + (pars['D1']*2 + pars['D2']) * steps) * (pars['NS'] + pars['DS'])
+            else:
+                seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * steps * (pars['NS']+ pars['DS'])
+            m, s = divmod(seconds, 60)
+            h, m = divmod(m, 60)
+            print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+            # loop for a variable parameter:
+            for index, pars[var_key] in enumerate(array):
+                print 'Arrayed stimulated-echo experiment for '+var_key+': run = '+str(index+1)+\
+                                 ' out of '+str(array.size)+': value = '+str(pars[var_key])
+                # loop for accumulation:
+                for run in xrange(pars['NS']+pars['DS']):
+                    yield ste_experiment(pars, run) 
+                synchronize()
+        else:
+           # estimate the experiment time:
+            seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * (pars['NS']+ pars['DS'])
+            m, s = divmod(seconds, 60)
+            h, m = divmod(m, 60)
+            print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield ste_experiment(pars, run)       
+                
+                
+# ------------------ Hahn-echo experiment settings ----------------------
+    if hahn_flag == True:
+     # set up acquisition parameters:
+        pars = {}
+        pars['P90'] = 1.6e-6	# 90-degree pulse length (s)
+        pars['SF'] = 62.92e6	# spectrometer frequency (Hz)
+        pars['O1'] = 0e3	# offset from SF (Hz)
+        pars['SW'] = 20e6	# spectral window (Hz)
+        pars['SI'] = 1*512	# number of acquisition points
+        pars['NS'] = 512	# number of scans
+        pars['DS'] = 0		# number of dummy scans
+        pars['RD'] = 1		# delay between scans (s)
+        pars['TAU'] = 10e-6	# echo delay (s)
+        pars['D4'] = 0e-6	# echo pre-acquisition delay (s)
+        pars['PHA'] = 170        # receiver phase (degree)
+        pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+        pars['OUTFILE'] = 'test'			# output file name
+
+      # specify a variable parameter (optional):
+        pars['VAR_PAR'] = 'TAU'	# variable parameter name (a string)   
+        start = 20e-6		# starting value
+        stop = 1e-3			# end value
+        steps = 16			# number of values
+        log_scale = True		# log scale flag 
+        stag_range = False		# staggered range flag
+
+      # check parameters for safety: 
+        if pars['PHA'] < 0:
+            pars['PHA'] = 360 + pars['PHA']
+
+        if pars['P90'] > 20e-6:
+            raise Exception("Pulse too long!!!")
+
+        # check whether a variable parameter is named: 
+        var_key = pars.get('VAR_PAR')
+        if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+            raise Exception("Pulse too long!!!")
+
+        if pars['NS']%8 != 0:
+            pars['NS'] = int(round(pars['NS'] / 8) + 1) * 8 
+            print 'Number of scans changed to ',pars['NS'],' due to phase cycling' 
+
+      # start the experiment:
+        if var_key: 
+            # this is an arrayed experiment:
+            if log_scale:
+                array = log_range(start,stop,steps)
+            else:
+                array = lin_range(start,stop,steps)
+
+            if stag_range:
+                array = staggered_range(array, size = 2)
+
+            # estimate the experiment time:
+            if var_key == 'TAU':
+                seconds = (sum(array)*2 + pars['RD'] * steps) * (pars['NS'] + pars['DS'])
+            elif var_key == 'RD':
+                seconds = (sum(array) + pars['TAU']*2 * steps) * (pars['NS'] + pars['DS'])
+            else:
+                seconds = (pars['TAU']*2 + pars['RD']) * steps * (pars['NS']+ pars['DS'])
+            m, s = divmod(seconds, 60)
+            h, m = divmod(m, 60)
+            print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+            # loop for a variable parameter:
+            for index, pars[var_key] in enumerate(array):
+                print 'Arrayed Hahn-echo experiment for '+var_key+': run = '+str(index+1)+\
+                                 ' out of '+str(array.size)+': value = '+str(pars[var_key])
+                # loop for accumulation:
+                for run in xrange(pars['NS']+pars['DS']):
+                    yield hahn_experiment(pars, run)
+                synchronize()
+        else:
+           # estimate the experiment time:
+            seconds = (pars['TAU']*2 + pars['RD']) * (pars['NS']+ pars['DS'])
+            m, s = divmod(seconds, 60)
+            h, m = divmod(m, 60)
+            print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield hahn_experiment(pars, run)        
+
+                 
+			
+# the pulse programs:
+
+def satrec2_experiment(pars, run):
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+    
+    pars['PROG'] = 'satrec2_experiment'
+    
+    # phase lists:
+    pars['PH1'] = [  0] 					# saturation pulses
+    pars['PH3'] = [ 0, 180,   0, 180, 90, 270,  90, 270]	# 1st 90-degree pulse
+    pars['PH4'] = [90,  90, 270, 270,  0,   0, 180, 180]	# 2nd 90-degree pulse
+    pars['PH2'] = [ 0, 180,   0, 180, 90, 270,  90, 270]	# receiver
+
+    # read in variables:
+    P90 = pars['P90']
+    SF =  pars['SF']
+    O1 =  pars['O1']
+    NECH = pars['NECH']
+    D1 =  pars['D1']
+    D2 =  pars['D2']
+    D3 =  pars['D3']
+    D4 =  pars['D4']
+    TAU = pars['TAU']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH4 = pars['PH4'][run%len(pars['PH4'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']
+
+    # set variable delay list for saturation pulses: 
+    vdlist = log_range(D2, D1, NECH-1)
+
+    # set sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+
+    # the pulse sequence: 
+
+    # saturation:
+    e.set_frequency(SF+O1, phase=PH1)			# set frequency and phase for saturation pulses
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# apply 90-degree pulse
+    for delay in vdlist[::-1]:
+        e.wait(delay-P90-TXEnableDelay)			# wait for next saturation pulse
+        e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+        e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue) # apply 90-degree pulse
+
+    # recovery:
+    e.wait(TAU)						# wait for tau
+    e.set_phase(PH3)					# set phase for next pulse
+     
+    # echo detection:
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# apply 90-degree pulse
+    e.wait(D3-P90/2-TXEnableDelay)			# echo delay
+    e.set_phase(PH4)					# set phase for next pulse
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# apply 90-degree pulse 
+    e.set_phase(PHA)					# set phase for receiver
+    e.wait(D3-P90/2+D4)					# echo delay    
+    e.record(SI, SW, sensitivity=ADCSensitivity)	# acquisition
+
+    # write experiment attributes:  
+    for key in pars.keys():
+        e.set_description(key, pars[key])	# pulse sequence parameters
+    e.set_description('run', run)		# current scan
+    e.set_description('rec_phase', -PH2)	# current receiver phase
+
+    return e    
+        
+
+def ste_experiment(pars, run):
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+		
+    pars['PROG'] = 'ste_experiment'
+    
+    # phase lists [16-phase cycle from JMR 157, 31 (2002)]:
+    pars['PH1'] = [0, 180,   0, 180,   0, 180,   0, 180,  90, 270,  90, 270,  90, 270,  90, 270] # 1st 90-degree pulse
+    pars['PH3'] = [0,   0, 180, 180,   0,   0, 180, 180,   0,   0, 180, 180,   0,   0, 180, 180] # 2nd 90-degree pulse
+    pars['PH4'] = [0,   0,   0,   0, 180, 180, 180, 180,   0,   0,   0,   0, 180, 180, 180, 180] # 3nd 90-degree pulse
+    pars['PH2'] = [0, 180, 180,   0, 180,   0,   0, 180, 270,  90,  90, 270,  90, 270, 270,  90] # receiver
+
+
+   # read in variables:
+    P90 = pars['P90']
+    SF =  pars['SF']
+    O1 =  pars['O1']    
+    RD =  pars['RD']
+    D1 =  pars['D1']
+    D2 =  pars['D2']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH4 = pars['PH4'][run%len(pars['PH4'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']
+
+   # set sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+
+   # run the pulse sequence:
+    e.wait(RD)                                  	# relaxation delay between scans
+    e.set_frequency(SF+O1, phase=PH1)
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse
+
+    e.wait(D1-P90/2-TXEnableDelay)			# 'short tau or tp'
+
+    e.set_phase(PH3)
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse 
+
+    e.wait(D2-P90/2-TXEnableDelay)			# 'long tau or tm'
+
+    e.set_phase(PH4)
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse
+
+    e.wait(TXEnableDelay)
+    e.set_phase(PHA)
+    e.wait(D1-P90/2-TXEnableDelay)                      # 'short tau or tp'
+    e.record(SI, SW, sensitivity=ADCSensitivity)	# acquisition
+
+   # write experiment parameters:
+    for key in pars.keys():
+        e.set_description(key, pars[key])	# acquisition parameters
+    e.set_description('run', run)		# current scan
+    e.set_description('rec_phase', -PH2)	# current receiver phase
+
+    return e
+    
+    
+def hahn_experiment(pars, run): 
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+		
+    pars['PROG'] = 'hahn_experiment'
+	
+    # phase lists [from Tecmag's pulse sequence]:
+    pars['PH1'] = [ 0, 180,   0, 180,  90, 270,  90, 270]	# 90-degree pulse
+    pars['PH3'] = [ 0,   0, 180, 180, 270, 270,  90,  90]	# 180-degree pulse
+    pars['PH2'] = [ 0, 180,   0, 180,  90, 270,  90, 270]	# receiver
+    
+    # read in variables:
+    P90 = pars['P90']
+    P180 = pars['P90']*2
+    SF = pars['SF']
+    O1 = pars['O1'] 
+    RD =  pars['RD']
+    TAU = pars['TAU']
+    D4 = pars['D4']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']
+
+    # set sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+
+    # run the pulse sequence: 
+    e.wait(RD)						# delay between scans
+    e.set_frequency(SF+O1, phase=PH1)			# set frequency and phase for 1st RF pulse
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# apply 1st 90-degree pulse
+    e.wait(TAU-P90/2-TXEnableDelay)			# wait for TAU
+    e.set_phase(PH3)					# set phase for 2nd 90-degree pulse
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enalble RF amplifier
+    e.ttl_pulse(P180, value=TXEnableValue|TXPulseValue)	# apply 2nd 90-degree pulse
+    e.set_phase(PHA)					# set phase for receiver
+    e.wait(TAU-P180/2+D4)				# wait for TAU       
+    e.record(SI, SW, sensitivity=ADCSensitivity)	# acquire echo points
+
+    # write the experiment parameters:  
+    for key in pars.keys():
+        e.set_description(key, pars[key])	# pulse sequence parameters
+    e.set_description('run', run)		# current scan
+    e.set_description('rec_phase', -PH2)	# current receiver phase
+
+    return e

AU_Programs/Diffusiometry/op_diff_res.py

@@ -0,0 +1,275 @@
+# -*- coding: iso-8859-1 -*-
+
+from numpy import *
+from scipy.signal import *
+from scipy.optimize import *
+from os import path, rename
+
+def result():
+
+    the_experiment = None	# current experiment's name
+    
+    measurements = {'satrec2_experiment': MeasurementResult('Saturation Recovery'),
+                        'ste_experiment': MeasurementResult('Stimulated Echo'),
+                       'hahn_experiment': MeasurementResult('Hahn Echo')}
+                               
+    measurement_ranges = {'satrec2_experiment': [0.5e-6, 4.5e-6],
+                              'ste_experiment': [1.5e-6, 4.5e-6],
+                             'hahn_experiment': [2.5e-6, 4.5e-6]}
+    measurement_ranging = True
+                             
+    # loop over the incoming results: 
+    for timesignal in results:
+        if not isinstance(timesignal,ADC_Result): 
+            continue     
+
+        # read experiment parameters:
+        pars = timesignal.get_description_dictionary()
+                    
+        # catch the actual experiment's name:
+        if the_experiment != pars.get('PROG'):
+           the_experiment = pars.get('PROG')
+           suffix = '' 	# output file name's suffix
+           counter = 1
+
+        # ---------------- digital filter ------------------ 
+             
+        # get actual sampling rate of timesignal:
+        sampling_rate = timesignal.get_sampling_rate()
+
+        # get user-defined spectrum width:
+        spec_width = pars['SW']
+        
+        # specify cutoff frequency, in relative units:
+        cutoff = spec_width / sampling_rate
+                    
+        # number of filter's coefficients:
+        numtaps = 29
+                
+        if cutoff < 1: # otherwise no filter applied
+            
+            # use firwin to create a lowpass FIR filter:
+            fir_coeff = firwin(numtaps, cutoff)
+            
+            # downsize x according to user-defined spectral window:
+            skip = int(sampling_rate / spec_width)
+            timesignal.x = timesignal.x[::skip]
+        
+            for i in range(2):
+                # apply the filter to ith channel:
+                timesignal.y[i] = lfilter(fir_coeff, 1.0, timesignal.y[i])
+
+                # zeroize first N-1 "corrupted" samples:
+                timesignal.y[i][:numtaps-1] = 0.0
+
+                # circular left shift of y:
+                timesignal.y[i] = roll(timesignal.y[i], -(numtaps-1))    
+
+                # downsize y to user-defined number of samples (SI):
+                timesignal.y[i] = timesignal.y[i][::skip]
+
+            # update the sampling_rate attribute of the signal's:
+            timesignal.set_sampling_rate(spec_width)
+
+        # ----------------------------------------------------
+
+        # rotate timesignal according to current receiver's phase:
+        timesignal.phase(pars['rec_phase'])
+
+        # provide timesignal to the display tab:
+        data['Current scan'] = timesignal
+
+        # accumulate...
+        if not locals().get('accu'):
+            accu = Accumulation()
+
+        # skip dummy scans, if any:
+        if pars['run'] < 0: continue
+
+        # add up:
+        accu += timesignal
+
+        # provide accumulation to the display tab:
+        data['Accumulation'] = accu
+
+        # check how many scans are done:
+        if accu.n == pars['NS']: # accumulation is complete
+           
+	    # make a copy:
+            echo = accu + 0
+
+            # compute the signal's phase:
+            #phi0 = arctan2(echo.y[1][0], echo.y[0][0]) * 180 / pi
+            #if not locals().get('ref'): ref = phi0   
+            #print 'phi0 = ', phi0
+
+            # rotate the signal to maximize Re (optional):
+            #echo.phase(-phi0)
+
+            # check whether it is an arrayed experiment: 
+            var_key = pars.get('VAR_PAR')
+            if var_key:
+                # get variable parameter's value:
+                var_value = pars.get(var_key)	
+
+                # store signals recorded for different var_values:
+                data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu
+
+		# estimate noise level:
+                if not locals().get('noise'):
+                    n = int(0.1*echo.x.size)
+                    noise = 3*std(echo.y[0][-n-numtaps:-1-numtaps])  				
+					
+                # measure echo intensity vs. var_value:	
+                if the_experiment in measurements.keys():
+				
+                    # option a: sum over the time interval specified:
+                    if measurement_ranging == True:
+                        [start, stop] = echo.get_sampling_rate() * array(measurement_ranges[the_experiment])
+                        measurements[the_experiment][var_value] = sum(echo.y[0][int(start):int(stop)]) 
+						
+                    # option b: sum of all samples above noise: 
+                    else:
+                        measurements[the_experiment][var_value] = sum(echo.y[0][echo.y[0]>noise])
+                    
+                    # store a measurement:
+                    data[measurements[the_experiment].get_title()] = measurements[the_experiment]
+                    
+                    # update the file name suffix:
+                    suffix = '_' + str(counter)
+                    counter += 1
+					
+                else: 
+                    print "Cannot recognize experiment: continue without measuring"				
+
+            # save accu if required:   
+            outfile = pars.get('OUTFILE')
+            if outfile: 
+                datadir = pars.get('DATADIR')
+
+                # write raw data in Simpson format:
+                filename = datadir+outfile+suffix+'.dat'
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.dat')
+                accu.write_to_simpson(filename)
+
+               # write parameters in a text file:
+                filename = datadir+outfile+suffix+'.par' 
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.par')
+
+                fileobject = open(filename, 'w')
+                for key in sorted(pars.iterkeys()):
+                    if key=='run': continue
+                    if key=='rec_phase': continue
+                    fileobject.write('%s%s%s'%(key,'=', pars[key]))
+                    fileobject.write('\n')
+                fileobject.close()
+
+            # reset accumulation:
+            del accu
+
+    if var_key == 'TAU' or var_key == 'D2':
+       # mono-exponential recovery fit:
+        try:
+            xdata = measurements['satrec2_experiment'].get_xdata()
+            ydata = measurements['satrec2_experiment'].get_ydata()
+            [amplitude, rate, offset] = fitfunc_recovery(xdata, ydata)
+            print '%s%02g' % ('Amplitude = ', amplitude)
+            print '%s%02g' % ('T1 [s] = ', 1./rate)
+
+            # update display for fit:      
+            measurements['satrec2_experiment'].y = func_recovery([amplitude, rate, offset], xdata)
+            data[measurements['satrec2_experiment'].get_title()] = measurements['satrec2_experiment']
+        except:
+            pass
+		
+        # mono-exponential decay fit to Hahn echoes:
+        try:
+            xdata = measurements['hahn_experiment'].get_xdata()
+            ydata = measurements['hahn_experiment'].get_ydata()
+            [amplitude, rate] = fitfunc_decay(xdata, ydata)
+            print 'Mono-exponential fit to Hahn echo decay:'
+            print '%s%02g' % ('Amplitude = ', amplitude)
+            print '%s%02g' % ('T2 [s] = ', 1./rate)
+
+            # update display for the fit:
+            measurements['hahn_experiment'].y = func_decay([amplitude, rate], xdata)
+            data[measurements['hahn_experiment'].get_title()] = measurements['hahn_experiment']
+        except:
+            pass
+		
+        try:		
+            # mono-exponential decay fit to stimulated echoes:
+            xdata = measurements['ste_experiment'].get_xdata()
+            ydata = measurements['ste_experiment'].get_ydata()
+            [amplitude, rate] = fitfunc_decay(xdata, ydata)
+            print 'Mono-exponential fit to stimulated echo decay:'
+            print '%s%02g' % ('Amplitude = ', amplitude)
+            print '%s%02g' % ('T2 [s] = ', 1./rate)
+
+            # update display for the fit:
+            measurements['ste_experiment'].y = func_decay([amplitude, rate], xdata)
+            data[measurements['ste_experiment'].get_title()] = measurements['ste_experiment']
+        except:
+            pass	
+
+# the fitting procedure for satrec:
+def fitfunc_recovery(xdata, ydata):
+
+    # initialize variable parameters:
+    try: 
+        # solve Az = b:
+        A = array((ones(xdata.size/2), xdata[0:xdata.size/2]))    
+        b = log(abs(mean(ydata[-2:]) - ydata[0:xdata.size/2]))     
+        z = linalg.lstsq(transpose(A), b)
+        amplitude = exp(z[0][0])
+        rate = -z[0][1]
+    except:
+        amplitude = abs(ydata[-1] - ydata[0])
+        rate = 1./(xdata[-1] - xdata[0])
+    offset = min(ydata)
+    p0 = [amplitude, rate, offset]
+
+    # run least-squares optimization:
+    plsq = leastsq(residuals_recovery, p0, args=(xdata, ydata))
+
+    return plsq[0] # best-fit parameters
+
+def residuals_recovery(p, xdata, ydata):
+    return ydata - func_recovery(p, xdata)
+
+# here is the function to fit
+def func_recovery(p, xdata):
+    return p[0]*(1-exp(-p[1]*xdata)) + p[2]
+
+			
+# the fitting procedure for hahn and ste:
+def fitfunc_decay(xdata, ydata): 
+
+    # initialize variable parameters:
+    try: 
+        # solve Az = b:
+        A = array((ones(xdata.size/2), xdata[0:xdata.size/2]))
+        b = log(abs(ydata[0:xdata.size/2]))
+        z = linalg.lstsq(transpose(A), b)
+        amplitude = exp(z[0][0])
+        rate = -z[0][1]
+    except:
+        amplitude = abs(ydata[0])
+        rate = 1./(xdata[-1] - xdata[0])
+    p0 = [amplitude, rate]
+
+    # run least-squares optimization:
+    plsq = leastsq(residuals_decay, p0, args=(xdata, ydata))
+
+    return plsq[0] # best-fit parameters
+
+def residuals_decay(p, xdata, ydata):
+    return ydata - func_decay(p, xdata)
+
+# here is the function to fit:
+def func_decay(p, xdata):
+    return p[0]*exp(-p[1]*xdata)
+
+pass

Scripts/CPMG/op_cpmg_exp.py

@@ -0,0 +1,156 @@
+# -*- coding: iso-8859-1 -*-
+
+TXEnableDelay = 2e-6
+TXEnableValue = 0b0001		# TTL line blanking RF amplifier (bit 0)
+TXPulseValue  = 0b0010		# TTL line triggering RF pulses  (bit 1)
+ADCSensitivity = 2		# voltage span for ADC
+
+def experiment(): # Carr-Purcell-Meiboom-Gill (CPMG) experiment
+    
+  # set up acquisition parameters:         
+    pars = {}
+    pars['P90'] = 1.7e-6	# 90-degree pulse length (s)
+    pars['SF'] = 338.7e6	# spectrometer frequency (Hz)
+    pars['O1'] = -60e3		# offset from SF (Hz)
+    pars['NS'] = 8		# number of scans
+    pars['DS'] = 0		# number of dummy scans
+    pars['RD'] = 3		# delay between scans (s)
+    pars['NECH'] = 16		# number of 180-degree pulses 
+    pars['TAU'] = 40e-6 	# half pulse period (s)
+    pars['PHA'] = -127		# receiver phase (degree)
+    pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+    pars['OUTFILE'] = None              	# output file name
+
+  # specify a variable parameter (optional):
+    pars['VAR_PAR'] = None	# variable parameter name (a string)   
+    start = 40e-6		# starting value
+    stop = 100e-6		# end value
+    steps = 10			# number of values
+    log_scale = False		# log-scale flag 
+    stag_range = False		# staggered range flag
+
+  # check parameters for safety: 
+    if pars['PHA'] < 0:
+        pars['PHA'] = 360 + pars['PHA']
+    
+    if pars['P90'] > 20e-6:
+        raise Exception("Pulse too long!!!")
+
+    # check whether a variable parameter is named: 
+    var_key = pars.get('VAR_PAR')
+    if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+        raise Exception("Pulse too long!!!")
+
+    if pars['NS']%4 != 0:
+        pars['NS'] = int(round(pars['NS'] / 4) + 1) * 4 
+        print 'Number of scans changed to ',pars['NS'],'due to phase cycling'
+
+  # start the experiment:
+    if var_key: 
+        # this is an arrayed experiment:
+        if log_scale:
+            array = log_range(start,stop,steps)
+        else:
+            array = lin_range(start,stop,steps) 
+			
+        if stag_range:
+            array = staggered_range(array, size = 2) 			
+
+        # estimate experiment time:
+        if var_key == 'TAU':
+            seconds = (sum(array)* 2* pars['NECH'] + pars['RD'] * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'NECH':
+            seconds = (pars['TAU']* 2* sum(array) + pars['RD'] * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'RD':
+            seconds = (pars['TAU']* 2* pars['NECH'] + sum(array)) * (pars['NS'] + pars['DS'])
+        else:
+            seconds = (pars['TAU']* 2* pars['NECH'] + pars['RD']) * steps * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for a variable parameter:
+        for index, pars[var_key] in enumerate(array):
+            print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\
+                             ' out of '+str(array.size)+': value = '+str(pars[var_key])							 
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield cpmg_experiment(pars, run)
+            synchronize()				
+				
+    else:
+        # estimate the experiment time:
+        seconds = (pars['TAU']* 2* pars['NECH'] + pars['RD']) * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for accumulation:
+        for run in xrange(pars['NS']+pars['DS']):
+            yield cpmg_experiment(pars, run)
+
+# the pulse program:
+
+def cpmg_experiment(pars, run):
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+
+    pars['PROG'] = 'cpmg_experiment'
+    
+    # phase lists:
+    pars['PH1'] = [0, 180, 90, 270]     # 90-degree pulse
+    pars['PH3'] = [90, 90, 180, 180]    # 180-degree pulse
+    pars['PH2'] = [0, 180, 90, 270]     # receiver
+
+    # read in variables:
+    P90 = pars['P90']
+    P180 = pars['P90']*2
+    SF =  pars['SF']
+    O1 =  pars['O1']
+    RD =  pars['RD']
+    NECH =  pars['NECH']
+    TAU = pars['TAU']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']
+
+    # set sampling parameters:
+    SI = 128    	# number of samples
+    SW = 20e6   	# sampling rate
+    AQ = (SI+6)/SW	# acquisition window
+
+    if TAU < (P90+P180)/2+TXEnableDelay or TAU < (P180+TXEnableDelay+AQ)/2:
+        raise Exception('pulse period is too short!')
+
+    if 2*TAU < P180+TXEnableDelay+SI/SW:
+        raise Exception('pulse period too short!')    
+
+    # run the pulse sequence: 
+    e.wait(RD)                                  	# delay between scans  
+    e.set_frequency(SF+O1, phase=PH1)           	# set frequency and phase for 90-degree pulse
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# apply 90-degree pulse
+    e.wait(TAU-P90/2-P180/2-TXEnableDelay)      	# wait for tau
+    e.set_phase(PH3)                            	# change phase for 180-degree pulse    
+
+    e.loop_start(NECH)                          	# ----- loop for echoes: -----
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P180, value=TXEnableValue|TXPulseValue)	# apply 180-degree pulse
+    e.set_phase(PHA)                            	# set phase for receiver     
+    e.wait(TAU-(P180+TXEnableDelay+AQ)/2)       	# pre-acquisition delay
+    e.record(SI, SW, timelength=AQ, sensitivity=ADCSensitivity)	# acquire echo samples
+    e.wait(TAU-(P180+TXEnableDelay+AQ)/2)       	# post-acquisition delay
+    e.set_phase(PH3)                            	# set phase for theta-degree pulse
+    e.loop_end()                                	# ----------------------------                             
+
+    # write experiment parameters:  
+    for key in pars.keys():
+        e.set_description(key, pars[key])	# acquisition parameters
+    e.set_description('run', run)       	# current scan
+    e.set_description('rec_phase', -PH2)	# current data route
+
+    return e

Scripts/CPMG/op_cpmg_res.py

@@ -0,0 +1,176 @@
+# -*- coding: iso-8859-1 -*-
+
+from numpy import *
+from scipy.signal import *
+from scipy.optimize import *
+from os import path, rename
+
+def result():
+  
+    measurement = MeasurementResult('Magnetization')
+  
+    suffix = '' 	# output file name's suffix and...
+    counter = 1		# counter for arrayed experiments
+        
+    # loop over the incoming results: 
+    for timesignal in results:
+        if not isinstance(timesignal,ADC_Result): 
+            continue   
+        
+        # read experiment parameters:
+        pars = timesignal.get_description_dictionary()
+        
+        # rotate timesignal by current receiver's phase:
+        timesignal.phase(pars['rec_phase'])
+
+        # provide timesignal to the display tab:
+        data['Current scan'] = timesignal
+
+        # accumulate...
+        if not locals().get('accu'):
+            accu = Accumulation()
+
+        # skip dummy scans, if any:
+        if pars['run'] < 0: continue
+
+        # add up:
+        accu += timesignal
+
+        # provide accumulation to the display tab:
+        data['Accumulation'] = accu
+
+        # check how many scans are done:
+        if accu.n == pars['NS']: # accumulation is complete
+
+            # get number of echoes:
+            num_echoes = pars['NECH'] 
+
+            # downsize accu to one point per echo:   
+            echodecay = accu + 0  		
+            echodecay.x = resize(echodecay.x, int(num_echoes))
+            echodecay.y[0] = resize(echodecay.y[0], int(num_echoes))
+            echodecay.y[1] = resize(echodecay.y[1], int(num_echoes))
+
+            # specify noise level: 
+            if not locals().get('noise'):
+                echo = accu.get_accu_by_index(0)
+                noise = 0.1*max(abs(echo.y[0]))
+                samples = abs(echo.y[0]) > noise
+
+            # set echo times and intensities:
+            for i in range(num_echoes):
+               # get ith echo:
+                echo = accu.get_accu_by_index(i)                          
+               # set echo timing:
+                echodecay.x[i] = i*2*pars['TAU']
+               # set echo value:
+                echodecay.y[0][i] = sum(echo.y[0][samples])	# the sum of echo points that exeed noise
+                echodecay.y[1][i] = sum(echo.y[1][samples]) 
+                #echodecay.y[0][i] = sum(echo.y[0])		# the sum of all echo points
+                #echodecay.y[1][i] = sum(echo.y[1])
+                #echodecay.y[0][i] = echo.y[0][echo.x.size/2]	# a middle echo point
+                #echodecay.y[1][i] = echo.y[1][echo.x.size/2]
+                
+            # compute a signal's phase:
+            phi0 = arctan2(echodecay.y[1][0], echodecay.y[0][0]) * 180 / pi
+            if not locals().get('ref'): ref = phi0   
+            print 'phi0 = ', phi0
+
+            # rotate signal to maximize Re (optional): 
+            #echodecay.phase(-phi0)
+
+            # provide echo decay to the display tab:
+            data['Echo Decay'] = echodecay
+
+            # fit a mono-exponential function to the echo decay:
+            [amplitude, rate] = fitfunc(echodecay.x, echodecay.y[0])
+            print '%s%02g' % ('Amplitude = ', amplitude)
+            print '%s%02g' % ('T2 [s] = ', 1./rate)
+
+            # provide the fit to the display tab:
+            fit = MeasurementResult('Mono-Exponential Fit')
+            for i, key in enumerate(echodecay.x):
+                fit[key] = echodecay.y[0][i]
+            fit.y = func([amplitude, rate], echodecay.x)
+            data[fit.get_title()] = fit
+                 
+            # check whether it is an arrayed experiment: 
+            var_key = pars.get('VAR_PAR')
+            if var_key:
+                # get variable parameter's value:
+                var_value = pars.get(var_key)
+                
+                # provide data recorded with different var_value's to the display tab:
+                data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu
+                data['Echo Decay'+"/"+var_key+"=%e"%(var_value)] = echodecay
+                data[fit.get_title()+"/"+var_key+"=%e"%(var_value)] = fit
+
+                # measure a signal parameter vs. var_value:
+                measurement[var_value] = amplitude
+                #measurement[var_value] = sum(echodecay.y[0][:])
+                #measurement[var_value] = 1./rate
+
+                # provide measurement to the display tab:
+                data[measurement.get_title()] = measurement
+
+            # save accu if required:   
+            outfile = pars.get('OUTFILE')
+            if outfile: 
+                datadir = pars.get('DATADIR')
+                                
+                # write data in Simpson format:
+                filename = datadir+outfile+suffix+'.dat'
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.dat')
+                accu.write_to_simpson(filename)
+
+                # write data in Tecmag format:
+                # filename = datadir+outfile+'.tnt'
+                # accu.write_to_tecmag(filename, nrecords=20)
+
+               # write parameters in a text file:
+                filename = datadir+outfile+suffix+'.par' 
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.par')
+
+                fileobject = open(filename, 'w')
+                for key in sorted(pars.iterkeys()):
+                    if key=='run': continue
+                    if key=='rec_phase': continue
+                    fileobject.write('%s%s%s'%(key,'=', pars[key]))
+                    fileobject.write('\n')
+                fileobject.close()
+
+            # reset accumulation:
+            del accu
+
+# the fitting procedure:
+def fitfunc(xdata, ydata): 
+
+    # initialize variable parameters:
+    try: 
+        # solve Az = b:
+        A = array((ones(xdata.size/2), xdata[0:xdata.size/2]))    
+        b = log(abs(ydata[0:xdata.size/2]))     
+        z = linalg.lstsq(transpose(A), b)
+        amplitude = exp(z[0][0])
+        rate = -z[0][1]
+    except:
+        amplitude = abs(ydata[0])
+        rate = 1./(xdata[-1] - xdata[0])
+    p0 = [amplitude, rate]
+
+    # run least-squares optimization:
+    plsq = leastsq(residuals, p0, args=(xdata, ydata))
+
+    return plsq[0] # best-fit parameters
+
+def residuals(p, xdata, ydata):
+    return ydata - func(p, xdata)
+
+# here is the function to fit:
+def func(p, xdata):
+    return p[0]*exp(-p[1]*xdata)
+
+
+pass

Scripts/FID/op_fid_exp.py

@@ -0,0 +1,137 @@
+# -*- coding: iso-8859-1 -*-
+
+TXEnableDelay = 2e-6
+TXEnableValue = 0b0001		# TTL line blanking RF amplifier (bit 0)
+TXPulseValue  = 0b0010		# TTL line triggering RF pulses  (bit 1)
+ADCSensitivity = 2		# voltage span for ADC
+
+def experiment(): # the basic pulse-acquire experiment
+
+  # set up acqusition parameters:
+    pars = {}
+    pars['P90'] = 1.7e-6	# 90-degree pulse length (s)
+    pars['SF'] = 338.7e6 	# spectrometer frequency (Hz)
+    pars['O1'] = -60e3 		# offset from SF (Hz)
+    pars['SW'] = 200e3		# spectral window (Hz)
+    pars['SI'] = 1*256		# number of acquisition points
+    pars['NS'] = 4		# number of scans
+    pars['DS'] = 0		# number of dummy scans
+    pars['RD'] = 3		# delay between scans (s)
+    pars['DEAD1'] = 5e-6	# receiver dead time (s)  
+    pars['PHA'] = 30		# receiver phase (degree)
+    pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+    pars['OUTFILE'] = None			# output file name
+
+   # specify a variable parameter (optional):
+    pars['VAR_PAR'] = None	# variable parameter name (a string)
+    start = 1e-6		# starting value
+    stop = 3e-6			# end value
+    steps = 4			# number of values
+    log_scale = False		# log scale flag 
+    stag_range = False		# staggered range flag
+
+  # check parameters for safety:
+    if pars['PHA'] < 0:
+        pars['PHA'] = 360 + pars['PHA']
+
+    if pars['P90'] > 20e-6:
+        raise Exception("Pulse too long!!!")
+
+    # check whether a variable parameter is named: 
+    var_key = pars.get('VAR_PAR')
+    if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+        raise Exception("Pulse too long!!!")
+
+    if pars['NS']%4 != 0:
+        pars['NS'] = int(round(pars['NS'] / 4) + 1) * 4 
+        print 'Number of scans changed to ', pars['NS'], 'due to phase cycling'
+
+    # start the experiment:
+    if var_key: 
+        # this is an arrayed experiment:
+   
+        if log_scale:
+            array = log_range(start,stop,steps)
+        else:
+            array = lin_range(start,stop,steps) 
+
+        if stag_range:
+            array = staggered_range(array, size = 2)
+
+        # estimate the experiment time:
+        if var_key == 'RD':
+            seconds = sum(array) * (pars['NS'] + pars['DS'])
+        else:
+            seconds = pars['RD'] * steps * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for a variable parameter:
+        for index, pars[var_key] in enumerate(array):
+            print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\
+                             ' out of '+str(array.size)+': value = '+str(pars[var_key])
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield fid_experiment(pars, run) 
+            synchronize()	
+    else:
+        # estimate the experiment time:
+        seconds = pars['RD'] * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for accumulation:
+        for run in xrange(pars['NS']+pars['DS']):
+            yield fid_experiment(pars, run)
+
+
+# the pulse program: 
+
+def fid_experiment(pars, run):
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+    
+    pars['PROG'] = 'fid_experiment'
+
+    # phase lists:
+    pars['PH1'] = [0, 180, 90, 270]	# 90-degree pulse
+    pars['PH2'] = [0, 180, 90, 270]	# receiver
+
+    # read in variables:
+    P90 = pars['P90'] 
+    SF = pars['SF']
+    O1 = pars['O1']   
+    RD = pars['RD']
+    DEAD1 = pars['DEAD1'] 
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']
+
+    # set sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+
+    # run the pulse sequence: 
+    e.wait(RD)                                          # delay between scans
+    e.set_frequency(SF+O1, phase=PH1)                   # set frequency and phase for RF pulse
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# apply RF pulse
+    e.set_phase(PHA)                                    # set phase for receiver 
+    e.wait(DEAD1)                                       # wait for coil ringdown
+    e.record(SI, SW, sensitivity=ADCSensitivity)	# acquire signal
+
+    # write experiment attributes:  
+    for key in pars.keys():
+        e.set_description(key, pars[key])		# acqusition parameters
+    e.set_description('run', run)			# current scan
+    e.set_description('rec_phase', -PH2)		# current receiver phase
+
+    return e

Scripts/FID/op_fid_res.py

@@ -0,0 +1,167 @@
+# -*- coding: iso-8859-1 -*-
+
+from numpy import *
+from scipy.signal import *
+from os import path, rename
+
+def result():
+
+    measurement = MeasurementResult('Magnetization')
+    
+    measurement_range = [0.0, 10e-6]
+    measurement_ranging = False
+
+    suffix = '' 	# output file name's suffix and...
+    counter = 1		# counter for arrayed experiments
+
+    # loop over the incoming results: 
+    for timesignal in results:
+        if not isinstance(timesignal,ADC_Result): 
+            continue     
+
+        # read experiment parameters:
+        pars = timesignal.get_description_dictionary()
+
+        # ---------------- digital filter ------------------ 
+             
+        # get actual sampling rate of timesignal:
+        sampling_rate = timesignal.get_sampling_rate()
+
+        # get user-defined spectrum width:
+        spec_width = pars['SW']
+        
+        # specify cutoff frequency, in relative units:
+        cutoff = spec_width / sampling_rate
+                
+        if cutoff < 1: # otherwise no filter applied
+
+            # number of filter's coefficients:
+            numtaps = 29
+            
+            # use firwin to create a lowpass FIR filter:
+            fir_coeff = firwin(numtaps, cutoff)
+            
+            # downsize x according to user-defined spectral window:
+            skip = int(sampling_rate / spec_width)
+            timesignal.x = timesignal.x[::skip]
+            
+            for i in range(2):
+                # apply the filter to ith channel:
+                timesignal.y[i] = lfilter(fir_coeff, 1.0, timesignal.y[i])
+
+                # zeroize first N-1 "corrupted" samples:
+                timesignal.y[i][:numtaps-1] = 0.0
+
+                # circular left shift of y:
+                timesignal.y[i] = roll(timesignal.y[i], -(numtaps-1))    
+
+                # downsize y to user-defined number of samples (SI):
+                timesignal.y[i] = timesignal.y[i][::skip]
+
+            # update the sampling_rate attribute of the signal's:
+            timesignal.set_sampling_rate(spec_width)
+
+        # ----------------------------------------------------
+
+        # phase timesignal according to current rec_phase:
+        timesignal.phase(pars['rec_phase'])
+
+        # provide timesignal to the display tab:
+        data['Current scan'] = timesignal
+
+        # accumulate...
+        if not locals().get('accu'):
+            accu = Accumulation()
+
+        # skip dummy scans, if any:
+        if pars['run'] < 0: continue
+
+        # add up:
+        accu += timesignal
+
+        # provide accumulation to the display tab:
+        data['Accumulation'] = accu
+
+        # check how many scans are done:
+        if accu.n == pars['NS']: # accumulatioin is complete
+
+            # make a copy:
+            fid = accu + 0
+
+            # compute the signal's phase:
+            phi0 = arctan2(fid.y[1][0], fid.y[0][0]) * 180 / pi
+            if not 'ref' in locals(): ref = phi0   
+            print 'phi0 = ', phi0
+
+            # rotate the signal to maximize Re (optional):     
+            #fid.phase(-phi0) 
+
+            # do FFT: 
+            fid.exp_window(line_broadening=10)
+            spectrum = fid.fft(samples=2*pars['SI'])
+
+            # try zero-order phase correction:
+            spectrum.phase(-phi0)
+            if abs(abs(phi0)-abs(ref)) > 90:    
+                spectrum.phase(180)
+
+            # provide spectrum to the display tab:
+            data['Spectrum'] = spectrum
+
+            # check whether it is an arrayed experiment: 
+            var_key = pars.get('VAR_PAR')
+            if var_key:
+                # get variable parameter's value:
+                var_value = pars.get(var_key)
+                
+                # provide signal recorded with this var_value to the display tab:
+                data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu
+                
+                # measure signal intensity vs. var_value:
+                fid_phased = (accu + 0).phase(-ref)
+                if measurement_ranging == True:
+                   [start, stop] = accu.get_sampling_rate() * array(measurement_range)
+                   measurement[var_value] = sum(fid_phased.y[0][int(start):int(stop)]) 
+			
+                else:
+                   measurement[var_value] = sum(fid_phased.y[0][0:31])
+
+                # provide measurements to the display tab:
+                data[measurement.get_title()] = measurement
+
+                # update the file name suffix:
+                suffix = '_' + str(counter)
+                counter += 1
+
+            # save accu if required:   
+            outfile = pars.get('OUTFILE')
+            if outfile: 
+                datadir = pars.get('DATADIR')
+
+                # write raw data in Simpson format:
+                filename = datadir+outfile+suffix+'.dat'
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.dat')
+                accu.write_to_simpson(filename)
+
+                # write raw data in Tecmag format:
+                # filename = datadir+outfile+'.tnt'
+                # accu.write_to_tecmag(filename, nrecords=20)
+
+               # write parameters in a text file:
+                filename = datadir+outfile+suffix+'.par' 
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.par')
+
+                fileobject = open(filename, 'w')
+                for key in sorted(pars.iterkeys()):
+                    if key=='run': continue
+                    if key=='rec_phase': continue
+                    fileobject.write('%s%s%s'%(key,'=', pars[key]))
+                    fileobject.write('\n')
+                fileobject.close()
+
+            # reset accumulation:
+            del accu
+
+pass

Scripts/FID_with_Background_Suppression/op_zgbs_exp.py

@@ -0,0 +1,147 @@
+# -*- coding: iso-8859-1 -*-
+
+TXEnableDelay = 2e-6
+TXEnableValue = 0b0001		# TTL line blanking RF amplifier (bit 0)
+TXPulseValue  = 0b0010		# TTL line triggering RF pulses  (bit 1)
+ADCSensitivity = 2		# voltage span for ADC
+
+def experiment(): # fid aquisition with background signal suppresion (the Bruker's zgbs)
+
+  # set up acqusition parameters:      
+    pars = {}
+    pars['P90'] = 3.5e-6	# 90-degree pulse length (s)
+    pars['SF'] = 360.0e6	# spectrometer frequency (Hz)
+    pars['O1'] = -140e3  	# offset from SF (Hz)
+    pars['SW'] = 10e6		# spectrum width (Hz)
+    pars['SI'] = 4*1024		# number of acquisition points
+    pars['NS'] = 64             # number of scans
+    pars['DS'] = 0              # number of dummy scans
+    pars['RD'] = 1              # delay between scans (s)
+    pars['DEAD1'] = 5e-6	# receiver dead time (s)  
+    pars['PHA'] = 0             # receiver phase (degree)
+    pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+    pars['OUTFILE'] = None                     	# output file name
+
+  # specify a variable parameter (optional):
+    pars['VAR_PAR'] = None	# variable parameter name (a string)   
+    start = 1e-6           	# starting value
+    stop = 5e-6            	# end value
+    steps = 5              	# number of values
+    log_scale = False		# log scale flag 
+    stag_range = False		# staggered range flag
+
+  # check parameters for safety:     
+    if pars['PHA'] < 0:
+        pars['PHA'] = 360 + pars['PHA']
+
+    if pars['P90'] > 20e-6:
+        raise Exception("Pulse too long!!!")
+                    
+    # check whether a variable parameter is named: 
+    var_key = pars.get('VAR_PAR')
+    if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+        raise Exception("Pulse too long!!!")
+
+    if pars['NS']%16 != 0:
+        pars['NS'] = int(round(pars['NS'] / 16) + 1) * 16 
+        print 'Number of scans changed to ', pars['NS'], ' due to phase cycling'
+
+    # start the experiment:
+    if var_key: 
+        # this is an arrayed experiment:
+        if log_scale:
+            array = log_range(start,stop,steps)
+        else:
+            array = lin_range(start,stop,steps) 
+
+        if stag_range:
+            array = staggered_range(array, size = 2)
+
+        # estimate the experiment time:
+        if var_key == 'RD':
+            seconds = sum(array) * (pars['NS'] + pars['DS'])
+        else:
+            seconds = pars['RD'] * steps * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+ 
+        # loop for a variable parameter:
+        for index, pars[var_key] in enumerate(array):
+            print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\
+                             ' out of '+str(array.size)+': value = '+str(pars[var_key])
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield zgbs_experiment(pars, run) 
+            synchronize()
+    else:
+        # estimate the experiment time:
+        seconds = pars['RD'] * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for accumulation:
+        for run in xrange(pars['NS']+pars['DS']):
+            yield zgbs_experiment(pars, run)
+
+
+# the pulse program: 
+
+def zgbs_experiment(pars, run):
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+
+    pars['PROG'] = 'zgbs_experiment'
+
+    # phase lists [from Bruker's zgbs]:
+    pars['PH1'] = [0, 0, 0, 0, 90, 90, 90, 90, 180, 180, 180, 180, 270, 270, 270, 270]	# 90-degree pulse
+    pars['PH3'] = [0, 90, 180, 270]							# 180-degree pulse
+    pars['PH4'] = [0, 0, 0, 0, 180, 180, 180, 180, 270, 270, 270, 270, 90, 90, 90, 90]	# 180-degree pulse
+    pars['PH2'] = [0, 180, 0, 180, 90, 270, 90, 270]					# receiver
+
+    # read in variables:
+    P90 = pars['P90'] 
+    P180 = 2*pars['P90']
+    SF = pars['SF']
+    O1 = pars['O1']   
+    RD = pars['RD']
+    DEAD1 = pars['DEAD1'] 
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH4 = pars['PH4'][run%len(pars['PH4'])]
+    PHA = pars['PHA']
+
+    # set sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+
+    # run the pulse sequence: 
+    e.wait(RD)                                          # delay between scans
+    e.set_frequency(SF+O1, phase=PH1)                   # set frequency and phase for RF pulse
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# apply 90-degree pulse
+    e.set_phase(PH3) 
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P180, value=TXEnableValue|TXPulseValue)	# apply 180-degree pulse
+    e.set_phase(PH4)    
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P180, value=TXEnableValue|TXPulseValue)	# apply 180-degree pulse
+    e.set_phase(PHA)					# set phase for receiver 
+    e.wait(DEAD1)					# wait for coil ringdown
+    e.record(SI, SW, sensitivity=ADCSensitivity)	# acquire signal
+
+    # write experiment parameters:  
+    for key in pars.keys():
+        e.set_description(key, pars[key])		# acqusition parameters
+    e.set_description('run', run)			# current scan
+    e.set_description('rec_phase', -PH2)		# current receiver phase
+
+    return e

Scripts/FID_with_Background_Suppression/op_zgbs_res.py

@@ -0,0 +1,167 @@
+# -*- coding: iso-8859-1 -*-
+
+from numpy import *
+from scipy.signal import *
+from os import path, rename
+
+def result():
+
+    measurement = MeasurementResult('Magnetization')
+    
+    measurement_range = [0.0, 10e-6]
+    measurement_ranging = False
+
+    suffix = '' 	# output file name's suffix and...
+    counter = 1		# counter for arrayed experiments
+
+    # loop over the incoming results: 
+    for timesignal in results:
+        if not isinstance(timesignal,ADC_Result): 
+            continue     
+
+        # read experiment parameters:
+        pars = timesignal.get_description_dictionary()
+
+        # ---------------- digital filter ------------------ 
+             
+        # get actual sampling rate of timesignal:
+        sampling_rate = timesignal.get_sampling_rate()
+
+        # get user-defined spectrum width:
+        spec_width = pars['SW']
+        
+        # specify cutoff frequency, in relative units:
+        cutoff = spec_width / sampling_rate
+                
+        if cutoff < 1: # otherwise no filter applied
+
+            # number of filter's coefficients:
+            numtaps = 29
+            
+            # use firwin to create a lowpass FIR filter:
+            fir_coeff = firwin(numtaps, cutoff)
+            
+            # downsize x according to user-defined spectral window:
+            skip = int(sampling_rate / spec_width)
+            timesignal.x = timesignal.x[::skip]
+            
+            for i in range(2):
+                # apply the filter to ith channel:
+                timesignal.y[i] = lfilter(fir_coeff, 1.0, timesignal.y[i])
+
+                # zeroize first N-1 "corrupted" samples:
+                timesignal.y[i][:numtaps-1] = 0.0
+
+                # circular left shift of y:
+                timesignal.y[i] = roll(timesignal.y[i], -(numtaps-1))    
+
+                # downsize y to user-defined number of samples (SI):
+                timesignal.y[i] = timesignal.y[i][::skip]
+
+            # update the sampling_rate attribute of the signal's:
+            timesignal.set_sampling_rate(spec_width)
+
+        # ----------------------------------------------------
+
+        # phase timesignal according to current rec_phase:
+        timesignal.phase(pars['rec_phase'])
+
+        # provide timesignal to the display tab:
+        data['Current scan'] = timesignal
+
+        # accumulate...
+        if not locals().get('accu'):
+            accu = Accumulation()
+
+        # skip dummy scans, if any:
+        if pars['run'] < 0: continue
+
+        # add up:
+        accu += timesignal
+
+        # provide accumulation to the display tab:
+        data['Accumulation'] = accu
+
+        # check how many scans are done:
+        if accu.n == pars['NS']: # accumulatioin is complete
+
+            # make a copy:
+            fid = accu + 0
+
+            # compute the signal's phase:
+            phi0 = arctan2(fid.y[1][0], fid.y[0][0]) * 180 / pi
+            if not 'ref' in locals(): ref = phi0   
+            print 'phi0 = ', phi0
+
+            # rotate the signal to maximize Re (optional):     
+            #fid.phase(-phi0) 
+
+            # do FFT: 
+            fid.exp_window(line_broadening=10)
+            spectrum = fid.fft(samples=2*pars['SI'])
+
+            # try zero-order phase correction:
+            spectrum.phase(-phi0)
+            if abs(abs(phi0)-abs(ref)) > 90:    
+                spectrum.phase(180)
+
+            # provide spectrum to the display tab:
+            data['Spectrum'] = spectrum
+
+            # check whether it is an arrayed experiment: 
+            var_key = pars.get('VAR_PAR')
+            if var_key:
+                # get variable parameter's value:
+                var_value = pars.get(var_key)
+                
+                # provide signal recorded with this var_value to the display tab:
+                data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu
+                
+                # measure signal intensity vs. var_value:
+                fid_phased = (accu + 0).phase(-ref)
+                if measurement_ranging == True:
+                   [start, stop] = accu.get_sampling_rate() * array(measurement_range)
+                   measurement[var_value] = sum(fid_phased.y[0][int(start):int(stop)]) 
+			
+                else:
+                   measurement[var_value] = sum(fid_phased.y[0][0:31])
+
+                # provide measurements to the display tab:
+                data[measurement.get_title()] = measurement
+
+                # update the file name suffix:
+                suffix = '_' + str(counter)
+                counter += 1
+
+            # save accu if required:   
+            outfile = pars.get('OUTFILE')
+            if outfile: 
+                datadir = pars.get('DATADIR')
+
+                # write raw data in Simpson format:
+                filename = datadir+outfile+suffix+'.dat'
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.dat')
+                accu.write_to_simpson(filename)
+
+                # write raw data in Tecmag format:
+                # filename = datadir+outfile+'.tnt'
+                # accu.write_to_tecmag(filename, nrecords=20)
+
+               # write parameters in a text file:
+                filename = datadir+outfile+suffix+'.par' 
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.par')
+
+                fileobject = open(filename, 'w')
+                for key in sorted(pars.iterkeys()):
+                    if key=='run': continue
+                    if key=='rec_phase': continue
+                    fileobject.write('%s%s%s'%(key,'=', pars[key]))
+                    fileobject.write('\n')
+                fileobject.close()
+
+            # reset accumulation:
+            del accu
+
+pass

Scripts/Hahn_Echo/op_hahn_exp.py

@@ -0,0 +1,147 @@
+# -*- coding: iso-8859-1 -*-
+
+TXEnableDelay = 2e-6
+TXEnableValue = 0b0001	# TTL line blanking RF amplifier (bit 0)
+TXPulseValue  = 0b0010	# TTL line triggering RF pulses  (bit 1)
+ADCSensitivity = 2	# voltage span for ADC
+
+def experiment(): # solid echo (quadrupolar echo) experiment
+
+  # set up acquisition parameters:
+    pars = {}
+    pars['P90'] = 1.7e-6	# 90-degree pulse length (s)
+    pars['SF'] = 338.7e6	# spectrometer frequency (Hz)
+    pars['O1'] = -60e3	 	# offset from SF (Hz)
+    pars['SW'] = 200e3		# spectral window (Hz)
+    pars['SI'] = 1*256		# number of acquisition points
+    pars['NS'] = 8		# number of scans
+    pars['DS'] = 0		# number of dummy scans
+    pars['RD'] = 3		# delay between scans (s)
+    pars['TAU'] = 13e-6		# echo delay (s)
+    pars['D4'] = 2e-6		# echo pre-acquisition delay (s)
+    pars['PHA'] = 0		# receiver phase (degree)
+    pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+    pars['OUTFILE'] = None			# output file name
+
+  # specify a variable parameter (optional):
+    pars['VAR_PAR'] = None	# variable parameter name (a string)   
+    start = 20e-6		# starting value
+    stop = 100e-6		# end value
+    steps = 12			# number of values
+    log_scale = False		# log scale flag 
+    stag_range = False		# staggered range flag
+
+  # check parameters for safety: 
+    if pars['PHA'] < 0:
+        pars['PHA'] = 360 + pars['PHA']
+
+    if pars['P90'] > 20e-6:
+        raise Exception("Pulse too long!!!")
+
+    # check whether a variable parameter is named: 
+    var_key = pars.get('VAR_PAR')
+    if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+        raise Exception("Pulse too long!!!")
+
+    if pars['NS']%8 != 0:
+        pars['NS'] = int(round(pars['NS'] / 8) + 1) * 8 
+        print 'Number of scans changed to ',pars['NS'],' due to phase cycling' 
+
+  # start the experiment:
+    if var_key: 
+        # this is an arrayed experiment:
+        if log_scale:
+            array = log_range(start,stop,steps)
+        else:
+            array = lin_range(start,stop,steps)
+
+        if stag_range:
+            array = staggered_range(array, size = 2)
+
+        # estimate the experiment time:
+        if var_key == 'TAU':
+            seconds = (sum(array)*2 + pars['RD'] * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'RD':
+            seconds = (sum(array) + pars['TAU']*2 * steps) * (pars['NS'] + pars['DS'])
+        else:
+            seconds = (pars['TAU']*2 + pars['RD']) * steps * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for a variable parameter:
+        for index, pars[var_key] in enumerate(array):
+            print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\
+                             ' out of '+str(array.size)+': value = '+str(pars[var_key])
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield hahn_experiment(pars, run)
+            synchronize()
+    else:
+       # estimate the experiment time:
+        seconds = (pars['TAU']*2 + pars['RD']) * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for accumulation:
+        for run in xrange(pars['NS']+pars['DS']):
+            yield hahn_experiment(pars, run)       
+
+
+# the pulse program:
+
+def hahn_experiment(pars, run): 
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+
+    pars['PROG'] = 'hahn_experiment'
+
+    # phase lists [from Tecmag's pulse sequence]:
+    pars['PH1'] = [ 0, 180,   0, 180,  90, 270,  90, 270]	# 90-degree pulse
+    pars['PH3'] = [ 0,   0, 180, 180, 270, 270,  90,  90]	# 180-degree pulse
+    pars['PH2'] = [ 0, 180,   0, 180,  90, 270,  90, 270]	# receiver
+
+    # read in variables:
+    P90 = pars['P90']
+    P180 = pars['P90']*2
+    SF = pars['SF']
+    O1 = pars['O1'] 
+    RD =  pars['RD']
+    TAU = pars['TAU']
+    D4 = pars['D4']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']
+
+    # set sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+
+    # run the pulse sequence: 
+    e.wait(RD)						# delay between scans
+    e.set_frequency(SF+O1, phase=PH1)			# set frequency and phase for 1st RF pulse
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# apply 1st 90-degree pulse
+    e.wait(TAU-P90/2-TXEnableDelay)			# wait for TAU
+    e.set_phase(PH3)					# set phase for 2nd 90-degree pulse
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enalble RF amplifier
+    e.ttl_pulse(P180, value=TXEnableValue|TXPulseValue)	# apply 2nd 90-degree pulse
+    e.set_phase(PHA)					# set phase for receiver
+    e.wait(TAU-P180/2+D4)				# wait for TAU       
+    e.record(SI, SW, sensitivity=ADCSensitivity)	# acquire echo points
+
+    # write the experiment parameters:  
+    for key in pars.keys():
+        e.set_description(key, pars[key])	# pulse sequence parameters
+    e.set_description('run', run)		# current scan
+    e.set_description('rec_phase', -PH2)	# current receiver phase
+
+    return e

Scripts/Hahn_Echo/op_hahn_res.py

@@ -0,0 +1,167 @@
+# -*- coding: iso-8859-1 -*-
+
+from numpy import *
+from scipy.signal import *
+from os import path, rename
+
+def result():
+
+    measurement = MeasurementResult('Magnetization')
+
+    measurement_range = [0.0, 10e-6]
+    measurement_ranging = False
+
+    suffix = '' 	# output file name's suffix and...
+    counter = 1		# counter for arrayed experiments
+
+    # loop over the incoming results: 
+    for timesignal in results:
+        if not isinstance(timesignal,ADC_Result): 
+            continue   
+
+        # read experiment parameters:
+        pars = timesignal.get_description_dictionary()
+
+        # ---------------- digital filter ------------------ 
+
+        # get actual sampling rate of timesignal:
+        sampling_rate = timesignal.get_sampling_rate()
+
+        # get user-defined spectrum width:
+        spec_width = pars['SW']
+
+        # specify cutoff frequency, in relative units:
+        cutoff = spec_width / sampling_rate
+
+        if cutoff < 1: # no filter applied otherwise
+                
+            # number of filter's coefficients:
+            numtaps = 29 
+
+            # use firwin to create a lowpass FIR filter:
+            fir_coeff = firwin(numtaps, cutoff)
+                
+            # downsize x according to user-defined spectral window:
+            skip = int(sampling_rate / spec_width)
+            timesignal.x = timesignal.x[::skip]
+
+            for i in range(2):
+                # apply the filter to ith channel:
+                timesignal.y[i] = lfilter(fir_coeff, 1.0, timesignal.y[i])
+
+                # zeroize first N-1 "corrupted" samples:
+                timesignal.y[i][:numtaps-1] = 0.0
+
+                # circular left shift of y:
+                timesignal.y[i] = roll(timesignal.y[i], -(numtaps-1))    
+
+                # downsize y to user-defined number of samples (SI):
+                timesignal.y[i] = timesignal.y[i][::skip]
+
+            # replace the sampling_rate attribute of the signal:
+            timesignal.set_sampling_rate(spec_width)
+
+        # ----------------------------------------------------
+
+        # rotate timesignal by current receiver's phase:
+        timesignal.phase(pars['rec_phase'])
+
+        # provide timesignal to the display tab:
+        data['Current scan'] = timesignal
+
+        # accumulate...
+        if not locals().get('accu'):
+            accu = Accumulation()
+
+        # skip dummy scans, if any:
+        if pars['run'] < 0: continue
+
+        # add up:
+        accu += timesignal
+
+        # provide accumulation to the display tab:
+        data['Accumulation'] = accu
+
+        # check how many scans are done:
+        if accu.n == pars['NS']: # accumulation is complete
+
+            # make a copy of signal:
+            echo = accu + 0
+
+            # compute the signal's phase:
+            phi0 = arctan2(echo.y[1][0], echo.y[0][0]) * 180 / pi
+            if not locals().get('ref'):  ref = phi0 
+            print 'phi0 = ', phi0 
+
+            # rotate to maximize Re (optional):     
+            #echo.phase(-phi0) 
+
+            # do FFT: 
+            echo.exp_window(line_broadening=10)
+            spectrum = echo.fft(samples=2*pars['SI'])
+
+            # try zero-order phase correction:
+            spectrum.phase(-phi0)
+            
+            # try baseline correction:
+            spectrum.baseline()
+
+            # provide spectrum to the display tab:
+            data['Spectrum'] = spectrum
+
+            # check whether it is an arrayed experiment: 
+            var_key = pars.get('VAR_PAR')
+            if var_key:
+                # get variable parameter's value:
+                var_value = pars.get(var_key)
+                
+                # provide signal recorded with this var_value to the display tab:
+                data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu                
+
+                # measure signal intensity vs. var_value:
+                if measurement_ranging == True:
+                   [start, stop] = accu.get_sampling_rate() * array(measurement_range)
+                   measurement[var_value] = sum(accu.y[0][int(start):int(stop)]) 
+                   
+                else:
+                   measurement[var_value] = sum(accu.y[0][0:31])
+
+                # provide measurement to the display tab:
+                data[measurement.get_title()] = measurement
+
+                # update the file name suffix:
+                suffix = '_' + str(counter)
+                counter += 1
+
+           # save accu if required:   
+            outfile = pars.get('OUTFILE')
+            if outfile: 
+                datadir = pars.get('DATADIR')
+
+                # write raw data in Simpson format:
+                filename = datadir+outfile+suffix+'.dat'
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.dat')
+                accu.write_to_simpson(filename)
+
+                # write raw data in Tecmag format:
+                # filename = datadir+outfile+'.tnt'
+                # accu.write_to_tecmag(filename, nrecords=20)
+
+               # write parameters in a text file:
+                filename = datadir+outfile+suffix+'.par' 
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.par')
+                    
+                fileobject = open(filename, 'w')
+                for key in sorted(pars.iterkeys()):
+                    if key=='run': continue
+                    if key=='rec_phase': continue
+                    fileobject.write('%s%s%s'%(key,'=', pars[key]))
+                    fileobject.write('\n')
+                fileobject.close()
+
+            # reset accumulation:
+            del accu
+
+pass

Scripts/Miscellaneous/op_gs_exp.py

@@ -0,0 +1,39 @@
+# -*- coding: iso-8859-1 -*-
+
+def experiment(): # 'go setup' routine
+    
+    pars = {}
+    pars['P90'] = 2.5e-6		# 90-degree pulse length (s)
+    pars['RD'] = 1              	# delay between scans (s)
+    pars['SW'] = 100e3			# spectrum window (Hz)
+    pars['SI'] = 1*1024			# number of acquisition points
+    pars['DEAD1'] = 10e-6		# receiver dead time (s)  
+
+    while True:
+        yield gs_experiment(pars)       
+
+
+def gs_experiment(pars):
+    e=Experiment()
+
+    # read in variables:
+    P90 = pars['P90']      
+    RD =  pars['RD']
+    SI =  pars['SI']
+    SW =  pars['SW']
+    DEAD1 =  pars['DEAD1'] 
+
+    if P90 > 20e-6:
+        raise Exception("Pulse too long!!!")
+
+    e.ttl_pulse(2e-6, value=1)		# unblank RF amplifier
+    e.ttl_pulse(P90, value=3)		# apply 90-degree pulse
+    e.wait(DEAD1)			# wait for receiver dead time
+    e.record(SI, SW, sensitivity=2)	# acquire signal
+    e.wait(RD)				# wait for next scan
+
+    # write experiment parameters:  
+    for key in pars.keys():
+        e.set_description(key, pars[key])
+
+    return e

Scripts/Miscellaneous/op_gs_res.py

@@ -0,0 +1,54 @@
+# -*- coding: iso-8859-1 -*-
+
+from numpy import *
+
+def result():
+        
+    # loop over the incoming results: 
+    for timesignal in results:
+        if not isinstance(timesignal,ADC_Result): 
+            continue   
+
+        # get list of parameters:
+        pars = timesignal.get_description_dictionary()
+
+        # make a copy of timesignal:
+        fid = timesignal + 0
+
+        # correct for the baseline offset:
+        fid.baseline()
+        
+        # compute the initial phase:
+        phi0 = arctan2(fid.y[1][0], fid.y[0][0]) * 180 / pi
+
+        # do FFT: 
+        fid.exp_window(line_broadening=10)
+        spectrum = fid.fft(samples=2*pars['SI'])
+        spectrum.baseline()
+        spectrum.phase(-phi0)
+
+        # provide timesignal and spectrum to the display tab:
+        data['Timesignal'] = timesignal
+        data["Spectrum"] = spectrum
+
+        # ------------- optional measurements: --------------
+
+        # measure signal intensity:
+        signal = (timesignal +0).phase(-phi0)
+        signal_intensity = sum(signal.y[0:10])      
+
+        # measure spectrum integral:
+        #spectrum_integral = cumsum(spectrum.y[0])
+        spectrum_magnitude = (spectrum+0).magnitude()
+        spectrum_integral =  cumsum(spectrum_magnitude.y[0]) 
+
+        # find the centre of gravity of the spectrum:        
+        cg = argwhere(spectrum_integral > 0.5*spectrum_integral[-1])[0] 
+
+        print ''
+        print '%s%.3e'%('FID intensity = ', signal_intensity)
+        print '%s%.3e'%('Spectrum integral = ', spectrum_integral[-1])
+        print '%s%g%s'%("Spectrum middle frequency = ", spectrum.x[cg], ' Hz')
+        print ''
+
+pass

Scripts/Saturation_Recovery/op_satrec_exp.py

@@ -0,0 +1,169 @@
+# -*- coding: iso-8859-1 -*-
+
+TXEnableDelay = 2e-6
+TXEnableValue = 0b0001	# TTL line blanking RF amplifier (bit 0)
+TXPulseValue  = 0b0010	# TTL line triggering RF pulses  (bit 1)
+ADCSensitivity = 2	# voltage span for ADC
+
+def experiment(): # saturation-recovery experiment
+
+  # set up acquisition parameters:
+    pars = {}
+    pars['P90'] = 1.7e-6	# 90-degree pulse length (s)
+    pars['SF'] = 338.7e6	# spectrometer frequency (Hz)
+    pars['O1'] = -60e3		# offset from SF (Hz)
+    pars['SW'] = 200e3		# spectral window (Hz)
+    pars['SI'] = 1*256		# number of acquisition points
+    pars['NS'] = 8		# number of scans
+    pars['DS'] = 0		# number of dummy scans
+    pars['TAU'] = 1		# delay for recovery (s)
+    pars['DEAD1'] = 5e-6	# receiver dead time (s)
+    pars['PHA'] = 100		# receiver phase (degree) 
+      # -*- these aren't variable: -*-
+    pars['NECH'] = 40		# number of saturation pulses
+    pars['D1'] = 100e-3		# starting interval in saturation sequence (s) 
+    pars['D2'] = 1e-4		# end interval in saturation sequence (s) 
+    pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+    pars['OUTFILE'] = None			# output file name
+
+  # specify a variable parameter (optional):
+    pars['VAR_PAR'] = 'TAU'	# variable parameter name (a string)   
+    start = 1e-3		# starting value
+    stop = 5e-0			# end value
+    steps = 10			# number of values
+    log_scale = True		# log scale flag 
+    stag_range = False		# staggered range flag
+
+  # check parameters for safety: 
+    if pars['PHA'] < 0:
+        pars['PHA'] = 360 + pars['PHA']
+
+    if pars['P90'] > 20e-6:
+        raise Exception("Pulse too long!!!")
+
+    # check whether a variable parameter is named: 
+    var_key = pars.get('VAR_PAR')
+    if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+        raise Exception("Pulse too long!!!")
+
+    if pars['NS']%4 != 0:
+        pars['NS'] = int(round(pars['NS'] / 4) + 1) * 4 
+        print 'Number of scans changed to ',pars['NS'],' due to phase cycling' 
+
+    if pars['D1'] < pars['D2']:
+        raise Exception("D1 must be greater than D2")
+
+    sat_length = sum(log_range(pars['D1'],pars['D2'],pars['NECH']))
+    if sat_length > 1.:
+        raise Exception("Saturation sequence too long!!!")   
+
+  # start the experiment:
+    if var_key: 
+        # this is an arrayed experiment:
+        if log_scale:
+            array = log_range(start,stop,steps)
+        else:
+            array = lin_range(start,stop,steps)
+
+        if stag_range:
+            array = staggered_range(array, size = 2)
+
+        # estimate the experiment time:
+        if var_key == 'TAU':
+            seconds = (sat_length * steps + sum(array)) * (pars['NS'] + pars['DS'])
+        else:
+            seconds = (sat_length + pars['TAU']) * steps * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for a variable parameter:
+        for index, pars[var_key] in enumerate(array):
+            print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\
+                             ' out of '+str(array.size)+': value = '+str(pars[var_key])
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield satrec_experiment(pars, run)
+            synchronize()
+
+    else:
+       # estimate the experiment time:
+        seconds = (sat_length + pars['TAU']) * (pars['NS'] + pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for accumulation:
+        for run in xrange(pars['NS']+pars['DS']):
+            yield satrec_experiment(pars, run) 
+
+
+# the pulse program:
+
+def satrec_experiment(pars, run):
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+        
+    pars['PROG'] = 'satrec_experiment'
+
+    # phase lists:
+    pars['PH1'] = [0]	 		# saturation pulses
+    pars['PH3'] = [0,180,90,270]	# measuring pulse 
+    pars['PH2'] = [0,180,90,270]	# receiver
+
+    # read in variables:
+    P90 = pars['P90']
+    SF =  pars['SF']
+    O1 =  pars['O1']
+    DEAD1 = pars['DEAD1'] 
+    NECH = pars['NECH']
+    D1 =  pars['D1']
+    D2 =  pars['D2']
+    TAU = pars['TAU']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']
+
+    # set sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+
+    # set variable delay list for saturation pulses: 
+    vdlist = log_range(D2, D1, NECH-1)
+
+    # run the pulse sequence:
+
+    # saturation:
+    e.set_frequency(SF+O1, phase=PH1)				# set frequency and phase for saturation pulses
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)		# enable RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)		# apply 90-degree pulse
+    for delay in vdlist[::-1]:
+        e.wait(delay-P90-TXEnableDelay)				# wait for next pulse
+        e.ttl_pulse(TXEnableDelay, value=TXEnableValue)		# enable RF amplifier
+        e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# apply 90-degree pulse
+
+    # recovery:
+    e.wait(TAU)							# recovery time
+    e.set_phase(PH3)						# set phase for measuring pulse
+
+    # detection:
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)		# enable RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)  	# apply 90-degree pulse  
+    e.set_phase(PHA)						# set phase for receiver
+    e.wait(DEAD1)                                               # wait for coil ringdown  
+    e.record(SI, SW, sensitivity=ADCSensitivity)		# acquire signal
+
+    # write experiment parameters:  
+    for key in pars.keys():
+        e.set_description(key, pars[key])			# acquisition parameters
+    e.set_description('run', run)				# current scan
+    e.set_description('rec_phase', -PH2)			# current receiver phase
+
+    return e

Scripts/Saturation_Recovery/op_satrec_res.py

@@ -0,0 +1,208 @@
+# -*- coding: iso-8859-1 -*-
+
+from numpy import *
+from scipy.signal import *
+from scipy.optimize import *
+from os import path, rename
+
+def result():
+
+    measurement = MeasurementResult('Magnetization')
+    
+    measurement_range = [0.0, 10e-6]
+    measurement_ranging = False    
+
+    suffix = ''		# output file name's suffix and...
+    counter = 1		# counter for arrayed experiments
+    var_key = ''	# variable parameter name
+
+    # loop over the incoming results: 
+    for timesignal in results:
+        if not isinstance(timesignal,ADC_Result): 
+            continue   
+
+        # read experiment parameters:
+        pars = timesignal.get_description_dictionary()
+
+        # ---------------- digital filter ------------------ 
+             
+        # get actual sampling rate of timesignal:
+        sampling_rate = timesignal.get_sampling_rate()
+
+        # get user-defined spectrum width:
+        spec_width = pars['SW']
+        
+        # specify cutoff frequency, in relative units:
+        cutoff = spec_width / sampling_rate
+                
+        if cutoff < 1: # no filter applied otherwise
+
+            # number of filter's coefficients:
+            numtaps = 29
+            
+            # use firwin to create a lowpass FIR filter:
+            fir_coeff = firwin(numtaps, cutoff)
+            
+            # downsize x according to user-defined spectral window:
+            skip = int(sampling_rate / spec_width)
+            timesignal.x = timesignal.x[::skip]
+            
+            for i in range(2):
+                # apply the filter to ith channel:
+                timesignal.y[i] = lfilter(fir_coeff, 1.0, timesignal.y[i])
+
+                # zeroize first N-1 "corrupted" samples:
+                timesignal.y[i][:numtaps-1] = 0.0
+
+                # circular left shift of y:
+                timesignal.y[i] = roll(timesignal.y[i], -(numtaps-1))    
+
+                # downsize y to user-defined number of samples (SI):
+                timesignal.y[i] = timesignal.y[i][::skip]
+
+            # update the sampling_rate attribute of the signal's:
+            timesignal.set_sampling_rate(spec_width)
+
+        # ----------------------------------------------------
+
+        # rotate timesignal according to current receiver's phase:
+        timesignal.phase(pars['rec_phase'])
+
+        # provide timesignal to the display tab:
+        data['Current Scan'] = timesignal
+
+        # accumulate...
+        if not locals().get('accu'):
+            accu = Accumulation()
+
+        # skip dummy scans, if any:
+        if pars['run'] < 0: continue
+
+        # add up:
+        accu += timesignal
+
+        # provide accumulation to the display tab:
+        data['Accumulation'] = accu
+
+        # check how many scans are done:
+        if accu.n == pars['NS']: # accumulation is complete
+
+            # make a copy:
+            fid = accu + 0
+
+            # compute the initial phase of FID:
+            phi0 = arctan2(fid.y[1][0], fid.y[0][0]) * 180 / pi
+            if not locals().get('ref'): ref = phi0 
+            print 'phi0 = ', phi0  
+
+            # rotate FID to maximize Re (optional):     
+            #fid.phase(-phi0) 
+
+            # do FFT: 
+            fid.exp_window(line_broadening=10)
+            spectrum = fid.fft(samples=2*pars['SI'])
+
+            # try zero-order phase correction:
+            spectrum.phase(-phi0)
+
+            # provide spectrum to the display tab:
+            data['Spectrum'] = spectrum
+
+            # check whether it is an arrayed experiment: 
+            var_key = pars.get('VAR_PAR')
+            if var_key:
+                # get variable parameter's value:
+                var_value = pars.get(var_key)
+
+                # provide signal recorded with this var_value to the display tab:
+                data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu
+
+                # measure signal intensity vs. var_value:
+                if measurement_ranging == True:
+                   [start, stop] = accu.get_sampling_rate() * array(measurement_range)
+                   measurement[var_value] = sum(accu.y[0][int(start):int(stop)]) 
+                        
+                else:
+                   measurement[var_value] = sum(accu.y[0][0:31])
+
+                # provide measurement to the display tab:
+                data[measurement.get_title()] = measurement
+
+                # update the file name suffix:
+                suffix = '_' + str(counter)
+                counter += 1
+
+            # save accu if required:   
+            outfile = pars.get('OUTFILE')
+            if outfile: 
+                datadir = pars.get('DATADIR')
+
+                # write raw data in Simpson format:
+                filename = datadir+outfile+suffix+'.dat'
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.dat')
+                accu.write_to_simpson(filename)
+
+                # write raw data in Tecmag format:
+                # filename = datadir+outfile+'.tnt'
+                # accu.write_to_tecmag(filename, nrecords=20)
+
+               # write parameters in a text file:
+                filename = datadir+outfile+suffix+'.par' 
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.par')
+
+                fileobject = open(filename, 'w')
+                for key in sorted(pars.iterkeys()):
+                    if key=='run': continue
+                    if key=='rec_phase': continue
+                    fileobject.write('%s%s%s'%(key,'=', pars[key]))
+                    fileobject.write('\n')
+                fileobject.close()
+
+            # reset accumulation:
+            del accu
+
+    if var_key == 'TAU':
+        # mono-exponential recovery fit:
+        xdata = measurement.get_xdata()
+        ydata = measurement.get_ydata()
+        [amplitude, rate, offset] = fitfunc(xdata, ydata)
+        print '%s%02g' % ('Amplitude = ', amplitude)
+        print '%s%02g' % ('T1 [s] = ', 1./rate)
+
+        # update display for the fit:
+        measurement.y = func([amplitude, rate, offset], xdata)
+        data[measurement.get_title()] = measurement
+
+# the fitting procedure:
+def fitfunc(xdata, ydata):
+
+    # initialize variable parameters:
+    try: 
+        # solve Az = b:
+        A = array((ones(xdata.size/2), xdata[0:xdata.size/2]))
+        b = log(abs(mean(ydata[-2:]) - ydata[0:xdata.size/2]))     
+        z = linalg.lstsq(transpose(A), b)
+        amplitude = exp(z[0][0])
+        rate = -z[0][1]
+    except:
+        amplitude = abs(ydata[-1] - ydata[0])
+        rate = 1./(xdata[-1] - xdata[0])
+    offset = min(ydata)
+    p0 = [amplitude, rate, offset]
+
+    # run least-squares optimization:
+    plsq = leastsq(residuals, p0, args=(xdata, ydata))
+
+    return plsq[0] # best-fit parameters
+
+def residuals(p, xdata, ydata):
+    return ydata - func(p, xdata)
+
+# here is the function to fit:
+def func(p, xdata):
+    return p[0]*(1-exp(-p[1]*xdata)) + p[2]
+
+
+pass

Scripts/Saturation_Recovery_with_Solid_Echo_Detection/op_satrec2_exp.py

@@ -0,0 +1,178 @@
+# -*- coding: iso-8859-1 -*-
+
+TXEnableDelay = 2e-6
+TXEnableValue = 0b0001	# TTL line blanking RF amplifier (bit 0)
+TXPulseValue  = 0b0010	# TTL line triggering RF pulses  (bit 1)
+ADCSensitivity = 2	# voltage span for ADC
+
+def experiment(): # saturation-recovery with soild-echo detection
+
+  # set up acquisition parameters:
+    pars = {}
+    pars['P90'] = 2.3e-6	# 90-degree pulse length (s)
+    pars['SF'] = 46.7e6		# spectrometer frequency (Hz)
+    pars['O1'] = 5.6e3		# offset from SF (Hz)
+    pars['SW'] = 10e6		# spectral window (Hz)
+    pars['SI'] = 1*1024		# number of acquisition points
+    pars['NS'] = 16		# number of scans
+    pars['DS'] = 0		# number of dummy scans
+    pars['TAU'] = 1		# delay for recovery (s)
+    pars['D3'] = 20e-6		# echo delay (s)
+    pars['D4'] = 0			# echo pre-aquisition delay (s)
+    pars['PHA'] = -30		# receiver phase (degree)
+     # -*- these ain't variable: -*-
+    pars['NECH'] = 40		# number of saturation pulses
+    pars['D1'] = 100e-3		# starting interval in saturation sequence (s) 
+    pars['D2'] = 1e-4		# end interval in saturation sequence (s) 
+    pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+    pars['OUTFILE'] = None			# output file name
+
+  # specify a variable parameter (optional):
+    pars['VAR_PAR'] = 'TAU'	# variable parameter name (a string)   
+    start = 1e-4		# starting value
+    stop = 1e-0			# end value
+    steps = 10			# number of values
+    log_scale = True		# log scale flag 
+    stag_range = True		# staggered range flag
+
+  # check parameters for safety: 
+    if pars['PHA'] < 0:
+        pars['PHA'] = 360 + pars['PHA']
+
+    if pars['P90'] > 20e-6:
+        raise Exception("pulse too long!!!")
+
+    # check whether a variable parameter is named: 
+    var_key = pars.get('VAR_PAR')
+    if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+        raise Exception("Pulse too long!!!")
+
+    if pars['NS']%8 != 0:
+        pars['NS'] = int(round(pars['NS'] / 8) + 1) * 8 
+        print 'Number of scans changed to ',pars['NS'],' due to phase cycling'
+
+    if pars['D1'] < pars['D2']:
+        raise Exception("D1 must be greater than D2!")
+
+    sat_length = sum(log_range(pars['D1'],pars['D2'],pars['NECH']))
+    if sat_length > 1.:
+        raise Exception("saturation sequence too long!!!")   
+
+    # start the experiment:
+    if var_key: 
+        # this is an arrayed experiment:
+        if log_scale:
+            array = log_range(start, stop, steps)
+        else:
+            array = lin_range(start, stop, steps)  
+
+        if stag_range:
+            array = staggered_range(array, size = 2)
+
+        # estimate the experiment time:
+        if var_key == 'TAU':
+            seconds = ((sat_length + pars['D3']*2) * steps + sum(array)) * (pars['NS'] + pars['DS'])
+        elif var_key == 'D3':
+            seconds = ((sat_length + pars['TAU']) * steps + sum(array)*2) * (pars['NS'] + pars['DS'])
+        else:
+            seconds = (sat_length + pars['TAU'] + pars['D3']*2) * steps * (pars['NS']+ pars['DS']) 
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for a variable parameter:
+        for index, pars[var_key] in enumerate(array):
+            print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\
+                             ' out of '+str(array.size)+': value = '+str(pars[var_key])
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield satrec2_experiment(pars, run) 
+            synchronize()
+    else:
+       # estimate the experiment time:
+        seconds = (sat_length + pars['TAU'] + pars['D3']*2) * (pars['NS']+ pars['DS']) 
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for accumulation:
+        for run in xrange(pars['NS']+pars['DS']):
+            yield satrec2_experiment(pars, run)
+
+
+# the pulse program: 
+
+def satrec2_experiment(pars, run):
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+
+    pars['PROG'] = 'satrec2_experiment'
+
+    # phase lists:
+    pars['PH1'] = [  0] 					# saturation pulses
+    pars['PH3'] = [ 0, 180,   0, 180, 90, 270,  90, 270]	# 1st 90-degree pulse
+    pars['PH4'] = [90,  90, 270, 270,  0,   0, 180, 180]	# 2nd 90-degree pulse
+    pars['PH2'] = [ 0, 180,   0, 180, 90, 270,  90, 270]	# receiver
+
+    # read in variables:
+    P90 = pars['P90']
+    SF =  pars['SF']
+    O1 =  pars['O1']
+    NECH = pars['NECH']
+    D1 =  pars['D1']
+    D2 =  pars['D2']
+    D3 =  pars['D3']
+    D4 =  pars['D4']
+    TAU = pars['TAU']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH4 = pars['PH4'][run%len(pars['PH4'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']
+
+    # set variable delay list for saturation pulses: 
+    vdlist = log_range(D2, D1, NECH-1)
+
+    # set sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+
+    # the pulse sequence: 
+
+    # saturation:
+    e.set_frequency(SF+O1, phase=PH1)			# set frequency and phase for saturation pulses
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# apply 90-degree pulse
+    for delay in vdlist[::-1]:
+        e.wait(delay-P90-TXEnableDelay)			# wait for next saturation pulse
+        e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+        e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue) # apply 90-degree pulse
+
+    # recovery:
+    e.wait(TAU)						# wait for tau
+    e.set_phase(PH3)					# set phase for next pulse
+     
+    # echo detection:
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# apply 90-degree pulse
+    e.wait(D3-P90/2-TXEnableDelay)			# echo delay
+    e.set_phase(PH4)					# set phase for next pulse
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# apply 90-degree pulse 
+    e.set_phase(PHA)					# set phase for receiver
+    e.wait(D3-P90/2+D4)					# echo delay    
+    e.record(SI, SW, sensitivity=ADCSensitivity)	# acquisition
+
+    # write experiment attributes:  
+    for key in pars.keys():
+        e.set_description(key, pars[key])	# pulse sequence parameters
+    e.set_description('run', run)		# current scan
+    e.set_description('rec_phase', -PH2)	# current receiver phase
+
+    return e

Scripts/Saturation_Recovery_with_Solid_Echo_Detection/op_satrec2_res.py

@@ -0,0 +1,208 @@
+# -*- coding: iso-8859-1 -*-
+
+from numpy import *
+from scipy.signal import *
+from scipy.optimize import *
+from os import path, rename
+
+def result():
+
+    measurement = MeasurementResult('Magnetization')
+    
+    measurement_range = [0.0, 10e-6]
+    measurement_ranging = False    
+
+    suffix = ''		# output file name's suffix and...
+    counter = 1		# counter for arrayed experiments
+    var_key = ''	# variable parameter name
+
+    # loop over the incoming results: 
+    for timesignal in results:
+        if not isinstance(timesignal,ADC_Result): 
+            continue   
+
+        # read experiment parameters:
+        pars = timesignal.get_description_dictionary()
+
+        # ---------------- digital filter ------------------ 
+             
+        # get actual sampling rate of timesignal:
+        sampling_rate = timesignal.get_sampling_rate()
+
+        # get user-defined spectrum width:
+        spec_width = pars['SW']
+        
+        # specify cutoff frequency, in relative units:
+        cutoff = spec_width / sampling_rate
+                
+        if cutoff < 1: # no filter applied otherwise
+
+            # number of filter's coefficients:
+            numtaps = 29
+            
+            # use firwin to create a lowpass FIR filter:
+            fir_coeff = firwin(numtaps, cutoff)
+            
+            # downsize x according to user-defined spectral window:
+            skip = int(sampling_rate / spec_width)
+            timesignal.x = timesignal.x[::skip]
+            
+            for i in range(2):
+                # apply the filter to ith channel:
+                timesignal.y[i] = lfilter(fir_coeff, 1.0, timesignal.y[i])
+
+                # zeroize first N-1 "corrupted" samples:
+                timesignal.y[i][:numtaps-1] = 0.0
+
+                # circular left shift of y:
+                timesignal.y[i] = roll(timesignal.y[i], -(numtaps-1))    
+
+                # downsize y to user-defined number of samples (SI):
+                timesignal.y[i] = timesignal.y[i][::skip]
+
+            # update the sampling_rate attribute of the signal's:
+            timesignal.set_sampling_rate(spec_width)
+
+        # ----------------------------------------------------
+
+        # rotate timesignal according to current receiver's phase:
+        timesignal.phase(pars['rec_phase'])
+
+        # provide timesignal to the display tab:
+        data['Current Scan'] = timesignal
+
+        # accumulate...
+        if not locals().get('accu'):
+            accu = Accumulation()
+
+        # skip dummy scans, if any:
+        if pars['run'] < 0: continue
+
+        # add up:
+        accu += timesignal
+
+        # provide accumulation to the display tab:
+        data['Accumulation'] = accu
+
+        # check how many scans are done:
+        if accu.n == pars['NS']: # accumulation is complete
+
+            # make a copy:
+            fid = accu + 0
+
+            # compute the initial phase of FID:
+            phi0 = arctan2(fid.y[1][0], fid.y[0][0]) * 180 / pi
+            if not locals().get('ref'): ref = phi0 
+            print 'phi0 = ', phi0  
+
+            # rotate FID to maximize Re (optional):     
+            #fid.phase(-phi0) 
+
+            # do FFT: 
+            fid.exp_window(line_broadening=10)
+            spectrum = fid.fft(samples=2*pars['SI'])
+
+            # try zero-order phase correction:
+            spectrum.phase(-phi0)
+
+            # provide spectrum to the display tab:
+            data['Spectrum'] = spectrum
+
+            # check whether it is an arrayed experiment: 
+            var_key = pars.get('VAR_PAR')
+            if var_key:
+                # get variable parameter's value:
+                var_value = pars.get(var_key)
+
+                # provide signal recorded with this var_value to the display tab:
+                data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu
+
+                # measure signal intensity vs. var_value:
+                if measurement_ranging == True:
+                   [start, stop] = accu.get_sampling_rate() * array(measurement_range)
+                   measurement[var_value] = sum(accu.y[0][int(start):int(stop)]) 
+                        
+                else:
+                   measurement[var_value] = sum(accu.y[0][0:31])
+
+                # provide measurement to the display tab:
+                data[measurement.get_title()] = measurement
+
+                # update the file name suffix:
+                suffix = '_' + str(counter)
+                counter += 1
+
+            # save accu if required:   
+            outfile = pars.get('OUTFILE')
+            if outfile: 
+                datadir = pars.get('DATADIR')
+
+                # write raw data in Simpson format:
+                filename = datadir+outfile+suffix+'.dat'
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.dat')
+                accu.write_to_simpson(filename)
+
+                # write raw data in Tecmag format:
+                # filename = datadir+outfile+'.tnt'
+                # accu.write_to_tecmag(filename, nrecords=20)
+
+               # write parameters in a text file:
+                filename = datadir+outfile+suffix+'.par' 
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.par')
+
+                fileobject = open(filename, 'w')
+                for key in sorted(pars.iterkeys()):
+                    if key=='run': continue
+                    if key=='rec_phase': continue
+                    fileobject.write('%s%s%s'%(key,'=', pars[key]))
+                    fileobject.write('\n')
+                fileobject.close()
+
+            # reset accumulation:
+            del accu
+
+    if var_key == 'TAU':
+        # mono-exponential recovery fit:
+        xdata = measurement.get_xdata()
+        ydata = measurement.get_ydata()
+        [amplitude, rate, offset] = fitfunc(xdata, ydata)
+        print '%s%02g' % ('Amplitude = ', amplitude)
+        print '%s%02g' % ('T1 [s] = ', 1./rate)
+
+        # update display for the fit:
+        measurement.y = func([amplitude, rate, offset], xdata)
+        data[measurement.get_title()] = measurement
+
+# the fitting procedure:
+def fitfunc(xdata, ydata):
+
+    # initialize variable parameters:
+    try: 
+        # solve Az = b:
+        A = array((ones(xdata.size/2), xdata[0:xdata.size/2]))
+        b = log(abs(mean(ydata[-2:]) - ydata[0:xdata.size/2]))     
+        z = linalg.lstsq(transpose(A), b)
+        amplitude = exp(z[0][0])
+        rate = -z[0][1]
+    except:
+        amplitude = abs(ydata[-1] - ydata[0])
+        rate = 1./(xdata[-1] - xdata[0])
+    offset = min(ydata)
+    p0 = [amplitude, rate, offset]
+
+    # run least-squares optimization:
+    plsq = leastsq(residuals, p0, args=(xdata, ydata))
+
+    return plsq[0] # best-fit parameters
+
+def residuals(p, xdata, ydata):
+    return ydata - func(p, xdata)
+
+# here is the function to fit:
+def func(p, xdata):
+    return p[0]*(1-exp(-p[1]*xdata)) + p[2]
+
+
+pass

Scripts/Solid_Echo/op_solidecho_exp.py

@@ -0,0 +1,147 @@
+# -*- coding: iso-8859-1 -*-
+
+TXEnableDelay = 2e-6
+TXEnableValue = 0b0001	# TTL line blanking RF amplifier (bit 0)
+TXPulseValue  = 0b0010	# TTL line triggering RF pulses  (bit 1)
+ADCSensitivity = 2	# voltage span for ADC
+
+def experiment(): # solid echo (quadrupolar echo) experiment
+
+  # set up acquisition parameters:
+    pars = {}
+    pars['P90'] = 1.7e-6	# 90-degree pulse length (s)
+    pars['SF'] = 338.7e6		# spectrometer frequency (Hz)
+    pars['O1'] = -60e3	 	# offset from SF (Hz)
+    pars['SW'] = 500e3		# spectral window (Hz)
+    pars['SI'] = 1*1024		# number of acquisition points
+    pars['NS'] = 8		# number of scans
+    pars['DS'] = 0		# number of dummy scans
+    pars['RD'] = 3		# delay between scans (s)
+    pars['TAU'] = 20e-6		# echo delay (s)
+    pars['D4'] = 0e-6		# echo pre-acquisition delay (s)
+    pars['PHA'] = -30		# receiver phase (degree)
+    pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+    pars['OUTFILE'] = None			# output file name
+
+  # specify a variable parameter (optional):
+    pars['VAR_PAR'] = None	# variable parameter name (a string)   
+    start = 1e-6		# starting value
+    stop = 3e-6			# end value
+    steps = 3			# number of values
+    log_scale = False		# log scale flag 
+    stag_range = False		# staggered range flag
+
+  # check parameters for safety: 
+    if pars['PHA'] < 0:
+        pars['PHA'] = 360 + pars['PHA']
+
+    if pars['P90'] > 20e-6:
+        raise Exception("Pulse too long!!!")
+
+    # check whether a variable parameter is named: 
+    var_key = pars.get('VAR_PAR')
+    if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+        raise Exception("Pulse too long!!!")
+
+    if pars['NS']%8 != 0:
+        pars['NS'] = int(round(pars['NS'] / 8) + 1) * 8 
+        print 'Number of scans changed to ',pars['NS'],' due to phase cycling' 
+
+  # start the experiment:
+    if var_key: 
+        # this is an arrayed experiment:
+        if log_scale:
+            array = log_range(start,stop,steps)
+        else:
+            array = lin_range(start,stop,steps)
+
+        if stag_range:
+            array = staggered_range(array, size = 2)
+
+        # estimate the experiment time:
+        if var_key == 'TAU':
+            seconds = (sum(array)*2 + pars['RD'] * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'RD':
+            seconds = (sum(array) + pars['TAU']*2 * steps) * (pars['NS'] + pars['DS'])
+        else:
+            seconds = (pars['TAU']*2 + pars['RD']) * steps * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for a variable parameter:
+        for index, pars[var_key] in enumerate(array):
+            print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\
+                             ' out of '+str(array.size)+': value = '+str(pars[var_key])
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield solidecho_experiment(pars, run)
+            synchronize()
+    else:
+       # estimate the experiment time:
+        seconds = (pars['TAU']*2 + pars['RD']) * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for accumulation:
+        for run in xrange(pars['NS']+pars['DS']):
+            yield solidecho_experiment(pars, run)       
+
+
+# the pulse program:
+
+def solidecho_experiment(pars, run): 
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+
+    pars['PROG'] = 'solidecho_experiment'
+
+    # phase lists [from Tecmag's pulse sequence]:
+    pars['PH1'] = [ 0, 180,   0, 180, 90, 270,  90, 270]	# 1st 90-degree pulse
+    pars['PH3'] = [90,  90, 270, 270,  0,   0, 180, 180]	# 2nd 90-degree pulse
+    pars['PH2'] = [ 0, 180,   0, 180, 90, 270,  90, 270]	# receiver
+
+
+    # read in variables:
+    P90 = pars['P90']
+    SF = pars['SF']
+    O1 = pars['O1'] 
+    RD =  pars['RD']
+    TAU = pars['TAU']
+    D4 = pars['D4']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']
+
+    # set sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+
+    # run the pulse sequence: 
+    e.wait(RD)						# delay between scans
+    e.set_frequency(SF+O1, phase=PH1)			# set frequency and phase for 1st RF pulse
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# apply 1st 90-degree pulse
+    e.wait(TAU-P90/2-TXEnableDelay)			# wait for TAU
+    e.set_phase(PH3)					# set phase for 2nd 90-degree pulse
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enalble RF amplifier
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# apply 2nd 90-degree pulse
+    e.set_phase(PHA)					# set phase for receiver
+    e.wait(TAU-P90/2+D4)				# wait for TAU       
+    e.record(SI, SW, sensitivity=ADCSensitivity)	# acquire echo points
+
+    # write the experiment parameters:  
+    for key in pars.keys():
+        e.set_description(key, pars[key])	# pulse sequence parameters
+    e.set_description('run', run)		# current scan
+    e.set_description('rec_phase', -PH2)	# current receiver phase
+
+    return e

Scripts/Solid_Echo/op_solidecho_res.py

@@ -0,0 +1,167 @@
+# -*- coding: iso-8859-1 -*-
+
+from numpy import *
+from scipy.signal import *
+from os import path, rename
+
+def result():
+
+    measurement = MeasurementResult('Magnetization')
+
+    measurement_range = [0.0, 10e-6]
+    measurement_ranging = False
+
+    suffix = '' 	# output file name's suffix and...
+    counter = 1		# counter for arrayed experiments
+
+    # loop over the incoming results: 
+    for timesignal in results:
+        if not isinstance(timesignal,ADC_Result): 
+            continue   
+
+        # read experiment parameters:
+        pars = timesignal.get_description_dictionary()
+
+        # ---------------- digital filter ------------------ 
+
+        # get actual sampling rate of timesignal:
+        sampling_rate = timesignal.get_sampling_rate()
+
+        # get user-defined spectrum width:
+        spec_width = pars['SW']
+
+        # specify cutoff frequency, in relative units:
+        cutoff = spec_width / sampling_rate
+
+        if cutoff < 1: # no filter applied otherwise
+                
+            # number of filter's coefficients:
+            numtaps = 29 
+
+            # use firwin to create a lowpass FIR filter:
+            fir_coeff = firwin(numtaps, cutoff)
+                
+            # downsize x according to user-defined spectral window:
+            skip = int(sampling_rate / spec_width)
+            timesignal.x = timesignal.x[::skip]
+
+            for i in range(2):
+                # apply the filter to ith channel:
+                timesignal.y[i] = lfilter(fir_coeff, 1.0, timesignal.y[i])
+
+                # zeroize first N-1 "corrupted" samples:
+                timesignal.y[i][:numtaps-1] = 0.0
+
+                # circular left shift of y:
+                timesignal.y[i] = roll(timesignal.y[i], -(numtaps-1))    
+
+                # downsize y to user-defined number of samples (SI):
+                timesignal.y[i] = timesignal.y[i][::skip]
+
+            # replace the sampling_rate attribute of the signal:
+            timesignal.set_sampling_rate(spec_width)
+
+        # ----------------------------------------------------
+
+        # rotate timesignal by current receiver's phase:
+        timesignal.phase(pars['rec_phase'])
+
+        # provide timesignal to the display tab:
+        data['Current scan'] = timesignal
+
+        # accumulate...
+        if not locals().get('accu'):
+            accu = Accumulation()
+
+        # skip dummy scans, if any:
+        if pars['run'] < 0: continue
+
+        # add up:
+        accu += timesignal
+
+        # provide accumulation to the display tab:
+        data['Accumulation'] = accu
+
+        # check how many scans are done:
+        if accu.n == pars['NS']: # accumulation is complete
+
+            # make a copy of signal:
+            echo = accu + 0
+
+            # compute the signal's phase:
+            phi0 = arctan2(echo.y[1][0], echo.y[0][0]) * 180 / pi
+            if not locals().get('ref'):  ref = phi0 
+            print 'phi0 = ', phi0 
+
+            # rotate to maximize Re (optional):     
+            #echo.phase(-phi0) 
+
+            # do FFT: 
+            echo.exp_window(line_broadening=10)
+            spectrum = echo.fft(samples=2*pars['SI'])
+
+            # try zero-order phase correction:
+            spectrum.phase(-phi0)
+            
+            # try baseline correction:
+            spectrum.baseline()
+
+            # provide spectrum to the display tab:
+            data['Spectrum'] = spectrum
+
+            # check whether it is an arrayed experiment: 
+            var_key = pars.get('VAR_PAR')
+            if var_key:
+                # get variable parameter's value:
+                var_value = pars.get(var_key)
+
+                # provide signal recorded with this var_value to the display tab:
+                data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu                
+
+                # measure signal intensity vs. var_value:
+                if measurement_ranging == True:
+                   [start, stop] = accu.get_sampling_rate() * array(measurement_range)
+                   measurement[var_value] = sum(accu.y[0][int(start):int(stop)]) 
+                   
+                else:
+                   measurement[var_value] = sum(accu.y[0][0:31])
+
+                # provide measurement to the display tab:
+                data[measurement.get_title()] = measurement
+
+                # update the file name suffix:
+                suffix = '_' + str(counter)
+                counter += 1
+
+           # save accu if required:   
+            outfile = pars.get('OUTFILE')
+            if outfile: 
+                datadir = pars.get('DATADIR')
+
+                # write raw data in Simpson format:
+                filename = datadir+outfile+suffix+'.dat'
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.dat')
+                accu.write_to_simpson(filename)
+
+                # write raw data in Tecmag format:
+                # filename = datadir+outfile+'.tnt'
+                # accu.write_to_tecmag(filename, nrecords=20)
+
+               # write parameters in a text file:
+                filename = datadir+outfile+suffix+'.par' 
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.par')
+                    
+                fileobject = open(filename, 'w')
+                for key in sorted(pars.iterkeys()):
+                    if key=='run': continue
+                    if key=='rec_phase': continue
+                    fileobject.write('%s%s%s'%(key,'=', pars[key]))
+                    fileobject.write('\n')
+                fileobject.close()
+
+            # reset accumulation:
+            del accu
+
+pass

Scripts/Spin_Alignment/op_spinal_exp.py

@@ -0,0 +1,163 @@
+# -*- coding: iso-8859-1 -*-
+
+TXEnableDelay = 2e-6
+TXEnableValue = 0b0001  # TTL line blanking RF amplifier (bit 0)
+TXPulseValue  = 0b0010  # TTL line triggering RF pulses  (bit 1)
+ADCSensitivity = 2	# voltage span for ADC
+
+def experiment(): # Jeener-Broekaert echo sequence (a.k.a. spin-alignment)
+
+  # set up acquisition parameters:         
+    pars = {}
+    pars['P90'] = 2.3e-6	# 90-degree pulse length (s)
+    pars['SF'] = 46.7e6		# spectrometer frequency (Hz)
+    pars['O1'] = -30e3	 	# offset from SF (Hz)
+    pars['SW'] = 1e6    	# spectral window (Hz)
+    pars['SI'] = 1*512  	# number of acquisition points
+    pars['NS'] = 8		# number of scans
+    pars['DS'] = 0		# number of dummy scans
+    pars['RD'] = 3		# delay between scans (s)
+    pars['D1'] = 30e-6		# delay after first pulse, or tp (s)
+    pars['D2'] = 100e-6		# delay after second pulse, or tm (s)
+    pars['PHA'] = -36		# receiver phase (degree)
+    pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+    pars['OUTFILE'] = None			# output file name
+    
+  # specify a variable parameter (optional):
+    pars['VAR_PAR'] = 'D2'	# variable parameter name (a string)   
+    start = 30e-6		# starting value
+    stop = 2e-0			# end value
+    steps = 24			# number of values
+    log_scale = True		# log scale flag 
+    stag_range = False		# staggered range flag
+
+  # check parameters for safety:
+    if pars['PHA'] < 0:
+        pars['PHA'] = 360 + pars['PHA']
+        
+    if pars['P90'] > 20e-6:
+        raise Exception("Pulse too long!!!")
+                    
+    # check whether a variable parameter is named: 
+    var_key = pars.get('VAR_PAR')
+    if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+        raise Exception("Pulse too long!!!")
+
+    if pars['NS']%8 != 0:
+        pars['NS'] = int(round(pars['NS'] / 8) + 1) * 8 
+        print 'Number of scans changed to ', pars['NS'], ' due to phase cycling'
+
+    # start the experiment:
+    if var_key: 
+        # this is an arrayed experiment:
+        if log_scale:
+            array = log_range(start,stop,steps)
+        else:
+            array = lin_range(start,stop,steps) 
+        
+        if stag_range:
+            array = staggered_range(array, size = 2)
+
+       # estimate the experiment time:
+        if var_key == 'D1':
+            seconds = (sum(array)*2 + (pars['D2'] + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'D2':
+            seconds = (sum(array) + (pars['D1']*2 + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'RD':
+            seconds = (sum(array) + (pars['D1']*2 + pars['D2']) * steps) * (pars['NS'] + pars['DS'])
+        else:
+            seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * steps * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for a variable parameter:
+        for index, pars[var_key] in enumerate(array):
+            print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\
+                             ' out of '+str(array.size)+': value = '+str(pars[var_key])
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield spinal_experiment(pars, run)
+            synchronize()
+
+    else:
+       # estimate the experiment time:
+        seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+       # loop for accumulation:
+        for run in xrange(pars['NS']+pars['DS']):
+            yield spinal_experiment(pars, run)
+
+
+# the pulse program: 
+
+def spinal_experiment(pars, run):
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+                       
+    pars['PROG'] = 'spinal_experiment'                       
+                       
+   # 8-step phase cycle (1-14 modifided to deal with T1-recovery and extended for Re/Im imbalance)
+    pars['PH1'] = [0, 270,   0, 270,     90,  90, 180, 180 ] # 1st (90-degree) pulse
+    pars['PH3'] = [90,180,  90, 180,    180, 180,  90,  90 ] # 2nd (45-degree) pulse
+    pars['PH4'] = [90, 90, 270, 270,    180,   0,   0, 180 ] # 3rd (45-degree) pulse
+    pars['PH2'] = [0, 180, 180,   0,     90, 270,  90, 270 ] # receiver
+        
+   # read in variables:
+    P90 = pars['P90']
+    P45 = pars['P90']*0.5
+    P1  = pars['P90']*0.5
+    SF =  pars['SF']
+    O1 =  pars['O1']    
+    RD =  pars['RD']
+    D1 =  pars['D1']
+    D2 =  pars['D2']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH4 = pars['PH4'][run%len(pars['PH4'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']
+
+   # set sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+
+   # run the pulse sequence:
+    e.wait(RD)						# relaxation delay between scans
+    e.set_frequency(SF+O1, phase=PH1)
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse 
+
+    e.wait(D1-P90/2-TXEnableDelay)			# 'short tau'
+    e.set_phase(PH3)
+
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)
+    e.ttl_pulse(P45, value=TXEnableValue|TXPulseValue)	# 45-degree pulse 
+
+    e.wait(D2-P45/2-TXEnableDelay)			# 'long tau'
+    e.set_phase(PH4)
+
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)
+    e.ttl_pulse(P1, value=TXEnableValue|TXPulseValue)	# 45-degree pulse
+
+    e.wait(TXEnableDelay)
+    e.set_phase(PHA)
+    e.wait(5e-6)#D1-P45/2-TXEnableDelay)    	                # 'short tau'
+    e.record(SI, SW, sensitivity=ADCSensitivity)        # acquisition                                     
+
+   # write experiment parameters:
+    for key in pars.keys():
+        e.set_description(key, pars[key])	# acquisition parameters
+    e.set_description('run', run)		# current scan
+    e.set_description('rec_phase', -PH2)	# current receiver phase
+
+    return e

Scripts/Spin_Alignment/op_spinal_res.py

@@ -0,0 +1,210 @@
+# -*- coding: iso-8859-1 -*-
+
+from numpy import *
+from scipy.signal import *
+from scipy.optimize import *
+from os import path, rename
+
+def result():
+
+    measurement = MeasurementResult('Magnetization')
+
+    measurement_range = [0.0, 10e-6]
+    measurement_ranging = False
+
+    suffix = '' 	# output file name's suffix and...
+    counter = 1		# counter for arrayed experiments
+    var_key = ''	# variable parameter name
+
+    # loop over the incoming results: 
+    for timesignal in results:
+        if not isinstance(timesignal,ADC_Result): 
+            continue
+
+        # read experiment parameters:
+        pars = timesignal.get_description_dictionary()
+
+        # ---------------- digital filter ------------------ 
+             
+        # get actual sampling rate of timesignal:
+        sampling_rate = timesignal.get_sampling_rate()
+
+        # get user-defined spectrum width:
+        spec_width = pars['SW']
+        
+        # specify cutoff frequency, in relative units:
+        cutoff = spec_width / sampling_rate
+                
+        if cutoff < 1: # no filter applied otherwise
+
+            # number of filter's coefficients:
+            numtaps = 29
+            
+            # use firwin to create a lowpass FIR filter:
+            fir_coeff = firwin(numtaps, cutoff)
+            
+            # downsize x according to user-defined spectral window:
+            skip = int(sampling_rate / spec_width)
+            timesignal.x = timesignal.x[::skip]
+            
+            for i in range(2):
+                # apply the filter to ith channel:
+                timesignal.y[i] = lfilter(fir_coeff, 1.0, timesignal.y[i])
+
+                # zeroize first N-1 "corrupted" samples:
+                timesignal.y[i][:numtaps-1] = 0.0
+
+                # circular left shift of y:
+                timesignal.y[i] = roll(timesignal.y[i], -(numtaps-1))    
+
+                # downsize y to user-defined number of samples (SI):
+                timesignal.y[i] = timesignal.y[i][::skip]
+
+            # update the sampling_rate attribute of the signal's:
+            timesignal.set_sampling_rate(spec_width)
+
+        # ----------------------------------------------------
+
+        # rotate timesignal according to current receiver's phase:
+        timesignal.phase(pars['rec_phase'])
+
+        # provide timesignal to the display tab:
+        data['Current scan'] = timesignal
+
+        # accumulate...
+        if not locals().get('accu'):
+            accu = Accumulation()
+
+        # skip dummy scans, if any:
+        if pars['run'] < 0: continue
+
+        # add up:
+        accu += timesignal
+
+        # provide accumulation to the display tab:
+        data['Accumulation'] = accu
+
+        # check how many scans are done:
+        if accu.n == pars['NS']: # accumulation is complete
+
+            # make a copy:
+            echo = accu + 0
+
+            # compute the initial phase of FID:
+            phi0 = arctan2(echo.y[1][0], echo.y[0][0]) * 180 / pi
+            if not 'ref' in locals(): ref = phi0   
+            print 'phi0 = ', phi0
+
+            # rotate FID to maximize y[0][0]:     
+            #echo.phase(-phi0) 
+
+            # do FFT:
+            echo.exp_window(line_broadening=10)
+            spectrum = echo.fft(samples=2*pars['SI'])
+
+            # try zero-order phase correction: 
+            spectrum.phase(-phi0)
+
+            # provide spectrum to the display tab:
+            data['Spectrum'] = spectrum
+
+            # check whether it is an arrayed experiment: 
+            var_key = pars.get('VAR_PAR')
+            if var_key:
+                # get variable parameter's value:
+                var_value = pars.get(var_key)
+
+                # provide signal recorded with this var_value to the display tab:
+                data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu                
+
+                # measure signal intensity vs. var_value:
+                if measurement_ranging == True:
+                   [start, stop] = accu.get_sampling_rate() * array(measurement_range)
+                   measurement[var_value] = sum(accu.y[0][int(start):int(stop)]) 
+                   
+                else:
+                   measurement[var_value] = sum(accu.y[0][0:31])
+
+                # provide measurement to the display tab:
+                data[measurement.get_title()] = measurement
+
+                # update the file name suffix:
+                suffix = '_' + str(counter)
+                counter += 1
+
+            # save accu if required:
+            outfile = pars.get('OUTFILE')
+            if outfile: 
+                datadir = pars.get('DATADIR')
+
+                # write raw data in Simpson format:
+                filename = datadir+outfile+suffix+'.dat'
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.dat')
+                accu.write_to_simpson(filename)
+
+                # write raw data in Tecmag format:
+                # filename = datadir+outfile+'.tnt'
+                # accu.write_to_tecmag(filename, nrecords=20)
+
+               # write parameters in a text file:
+                filename = datadir+outfile+suffix+'.par' 
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.par')
+
+                fileobject = open(filename, 'w')
+                for key in sorted(pars.iterkeys()):
+                    if key=='run': continue
+                    if key=='rec_phase': continue
+                    fileobject.write('%s%s%s'%(key,'=', pars[key]))
+                    fileobject.write('\n')
+                fileobject.close()
+
+            # reset accumulation:
+            del accu
+
+
+    if var_key == 'D2':
+        # KWW fit:
+        xdata = measurement.get_xdata()
+        ydata = measurement.get_ydata()
+        [amplitude, rate, beta] = fitfunc(xdata, ydata)
+        print '%s%02g' % ('Amplitude = ', amplitude)
+        print '%s%02g' % ('T2 [s] = ', 1./rate)
+        print '%s%02g' % ('Beta = ', beta)
+
+        # update display for the fit:
+        measurement.y = func([amplitude, rate, beta], xdata)
+        data[measurement.get_title()] = measurement
+
+# the fitting procedure:
+def fitfunc(xdata, ydata): 
+
+    # initialize variable parameters:
+    try: 
+        # solve Az = b:
+        A = array((ones(xdata.size/2), xdata[0:xdata.size/2]))    
+        b = log(abs(ydata[0:xdata.size/2]))     
+        z = linalg.lstsq(transpose(A), b)
+        amplitude = exp(z[0][0])
+        rate = -z[0][1]
+    except:
+        amplitude = abs(ydata[0])
+        rate = 1./(xdata[-1] - xdata[0])
+    beta = 1
+    p0 = [amplitude, rate, beta]
+
+    # run least-squares optimization:
+    plsq = leastsq(residuals, p0, args=(xdata, ydata))
+
+    return plsq[0] # best-fit parameters
+
+def residuals(p, xdata, ydata):
+    return ydata - func(p, xdata)
+
+# here is the function to fit:
+def func(p, xdata):
+    return p[0]*exp(-(p[1]*xdata)**p[2])
+
+
+pass

Scripts/Spin_Alignment_Spin32/op_spinal32_exp.py

@@ -0,0 +1,163 @@
+# -*- coding: iso-8859-1 -*-
+
+TXEnableDelay = 2e-6
+TXEnableValue = 0b0001  # TTL line blanking RF amplifier (bit 0)
+TXPulseValue  = 0b0010  # TTL line triggering RF pulses  (bit 1)
+ADCSensitivity = 2	# voltage span for ADC
+
+def experiment(): # Jeener-Broekaert echoes (a.k.a. spin-alignment) for spins-3/2
+    
+  # set up acquisition parameters:         
+    pars = {}
+    pars['P90'] = 2.0e-6	# 90-degree pulse length (s)
+    pars['SF'] = 139.9e6 	# spectrometer frequency (Hz)
+    pars['O1'] = -31e3    	# offset from SF (Hz)
+    pars['SW'] = 10e6  	 	# spectral window (Hz)
+    pars['SI'] = 4*1024  	# number of acquisition points
+    pars['NS'] = 32     	# number of scans
+    pars['DS'] = 0          	# number of dummy scans
+    pars['RD'] = 25            	# delay between scans (s)
+    pars['D1'] = 30e-6  	# delay after first pulse, or short tau (s)
+    pars['D2'] = 10e-6          # delay after second pulse, or long tau (s) 
+    pars['PHA'] = 65	        # receiver phase (degree)   
+    pars['DATADIR'] = '/home/fprak/Desktop/'	# data directory
+    pars['OUTFILE'] = None      		# output file name
+    
+  # specify a variable parameter (optional):
+    pars['VAR_PAR'] = 'D2'	# variable parameter name (a string)   
+    start = 10e-6            	# starting value
+    stop = 1            	# end value
+    steps = 24               	# number of values
+    log_scale = True		# log scale flag 
+    stag_range = False		# staggered range flag
+
+  # check parameters for safety:
+    if pars['PHA'] < 0:
+        pars['PHA'] = 360 + pars['PHA']
+        
+    if pars['P90'] > 20e-6:
+        raise Exception("Pulse too long!!!")
+                    
+    # check whether a variable parameter is named: 
+    var_key = pars.get('VAR_PAR')
+    if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+        raise Exception("Pulse too long!!!")
+
+    if pars['NS']%8 != 0:
+        pars['NS'] = int(round(pars['NS'] / 8) + 1) * 8 
+        print 'Number of scans changed to ', pars['NS'], ' due to phase cycling'
+
+    # start the experiment:
+    if var_key: 
+        # this is an arrayed experiment:
+        if log_scale:
+            array = log_range(start,stop,steps)
+        else:
+            array = lin_range(start,stop,steps) 
+        
+        if stag_range:
+            array = staggered_range(array, size = 2)
+
+       # estimate the experiment time:
+        if var_key == 'D1':
+            seconds = (sum(array)*2 + (pars['D2'] + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'D2':
+            seconds = (sum(array) + (pars['D1']*2 + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'RD':
+            seconds = (sum(array) + (pars['D1']*2 + pars['D2']) * steps) * (pars['NS'] + pars['DS'])
+        else:
+            seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * steps * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for a variable parameter:
+        for index, pars[var_key] in enumerate(array):
+            print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\
+                             ' out of '+str(array.size)+': value = '+str(pars[var_key])
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield spinal32_experiment(pars, run)
+            synchronize()
+
+    else:
+       # estimate the experiment time:
+        seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+       # loop for accumulation:
+        for run in xrange(pars['NS']+pars['DS']):
+            yield spinal32_experiment(pars, run)
+
+
+# the pulse program: 
+
+def spinal32_experiment(pars, run):
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+        
+    pars['PROG'] = 'spinal32_experiment'
+	
+        #  phase cycle by F. Qi et al. [JMR 169 (2004) 225-239] with 3rd-phase invertion
+    pars['PH1'] = [0, 180,   0, 180,    90, 270,  90, 270,         90, 270, 90, 270,    180,   0, 180,   0]	
+    pars['PH3'] = [90, 90, 270, 270,     0,   0, 180, 180,        180, 180,  0,   0,     90,  90, 270, 270]	
+    pars['PH4'] = [0,   0,   0,   0,   180, 180, 180, 180,         90,  90, 90,  90,    270, 270, 270, 270]	
+    pars['PH2'] = [180, 0,   0, 180,   180,   0,   0, 180,        270,  90, 90, 270,    270,  90,  90, 270]	
+
+
+   # read in variables:
+    P90 = pars['P90']
+    P45 = pars['P90']*0.5
+    SF =  pars['SF']
+    O1 =  pars['O1']    
+    RD =  pars['RD']
+    D1 =  pars['D1']
+    D2 =  pars['D2']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH4 = pars['PH4'][run%len(pars['PH4'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']
+    
+   # set sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+    
+   # run the pulse program: 
+    e.wait(RD)   				        # relaxation delay between scans
+    e.set_frequency(SF+O1, phase=PH1)           				 
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)		
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse 
+
+    e.wait(D1-P90/2-TXEnableDelay)			# 'short tau'
+    e.set_phase(PH3)
+    				
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)		
+    e.ttl_pulse(P45, value=TXEnableValue|TXPulseValue)	# 45-degree pulse 
+   
+    e.wait(D2-P45/2-TXEnableDelay)			# 'long tau'
+    e.set_phase(PH4)
+
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)		
+    e.ttl_pulse(P45, value=TXEnableValue|TXPulseValue)	# 45-degree pulse
+    
+    e.wait(TXEnableDelay)
+    e.set_phase(PHA)
+    e.wait(D1-P45/2-TXEnableDelay)    	                # 'short tau'
+    e.record(SI, SW, sensitivity=ADCSensitivity)	# acquisition                                     
+  
+   # write experiment parameters:  
+    for key in pars.keys():
+        e.set_description(key, pars[key])	# acqusition parameters
+    e.set_description('run', run)		# current scan
+    e.set_description('rec_phase', -PH2)        # current receiver phase
+
+    return e

Scripts/Spin_Alignment_Spin32/op_spinal32_res.py

@@ -0,0 +1,210 @@
+# -*- coding: iso-8859-1 -*-
+
+from numpy import *
+from scipy.signal import *
+from scipy.optimize import *
+from os import path, rename
+
+def result():
+
+    measurement = MeasurementResult('Magnetization')
+
+    measurement_range = [0.0, 10e-6]
+    measurement_ranging = False
+
+    suffix = '' 	# output file name's suffix and...
+    counter = 1		# counter for arrayed experiments
+    var_key = ''	# variable parameter name
+
+    # loop over the incoming results: 
+    for timesignal in results:
+        if not isinstance(timesignal,ADC_Result): 
+            continue
+
+        # read experiment parameters:
+        pars = timesignal.get_description_dictionary()
+
+        # ---------------- digital filter ------------------ 
+             
+        # get actual sampling rate of timesignal:
+        sampling_rate = timesignal.get_sampling_rate()
+
+        # get user-defined spectrum width:
+        spec_width = pars['SW']
+        
+        # specify cutoff frequency, in relative units:
+        cutoff = spec_width / sampling_rate
+                
+        if cutoff < 1: # no filter applied otherwise
+
+            # number of filter's coefficients:
+            numtaps = 29
+            
+            # use firwin to create a lowpass FIR filter:
+            fir_coeff = firwin(numtaps, cutoff)
+            
+            # downsize x according to user-defined spectral window:
+            skip = int(sampling_rate / spec_width)
+            timesignal.x = timesignal.x[::skip]
+            
+            for i in range(2):
+                # apply the filter to ith channel:
+                timesignal.y[i] = lfilter(fir_coeff, 1.0, timesignal.y[i])
+
+                # zeroize first N-1 "corrupted" samples:
+                timesignal.y[i][:numtaps-1] = 0.0
+
+                # circular left shift of y:
+                timesignal.y[i] = roll(timesignal.y[i], -(numtaps-1))    
+
+                # downsize y to user-defined number of samples (SI):
+                timesignal.y[i] = timesignal.y[i][::skip]
+
+            # update the sampling_rate attribute of the signal's:
+            timesignal.set_sampling_rate(spec_width)
+
+        # ----------------------------------------------------
+
+        # rotate timesignal according to current receiver's phase:
+        timesignal.phase(pars['rec_phase'])
+
+        # provide timesignal to the display tab:
+        data['Current scan'] = timesignal
+
+        # accumulate...
+        if not locals().get('accu'):
+            accu = Accumulation()
+
+        # skip dummy scans, if any:
+        if pars['run'] < 0: continue
+
+        # add up:
+        accu += timesignal
+
+        # provide accumulation to the display tab:
+        data['Accumulation'] = accu
+
+        # check how many scans are done:
+        if accu.n == pars['NS']: # accumulation is complete
+
+            # make a copy:
+            echo = accu + 0
+
+            # compute the initial phase of FID:
+            phi0 = arctan2(echo.y[1][0], echo.y[0][0]) * 180 / pi
+            if not 'ref' in locals(): ref = phi0   
+            print 'phi0 = ', phi0
+
+            # rotate FID to maximize y[0][0]:     
+            #echo.phase(-phi0) 
+
+            # do FFT:
+            echo.exp_window(line_broadening=10)
+            spectrum = echo.fft(samples=2*pars['SI'])
+
+            # try zero-order phase correction: 
+            spectrum.phase(-phi0)
+
+            # provide spectrum to the display tab:
+            data['Spectrum'] = spectrum
+
+            # check whether it is an arrayed experiment: 
+            var_key = pars.get('VAR_PAR')
+            if var_key:
+                # get variable parameter's value:
+                var_value = pars.get(var_key)
+
+                # provide signal recorded with this var_value to the display tab:
+                data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu                
+
+                # measure signal intensity vs. var_value:
+                if measurement_ranging == True:
+                   [start, stop] = accu.get_sampling_rate() * array(measurement_range)
+                   measurement[var_value] = sum(accu.y[0][int(start):int(stop)]) 
+                   
+                else:
+                   measurement[var_value] = sum(accu.y[0][0:31])
+
+                # provide measurement to the display tab:
+                data[measurement.get_title()] = measurement
+
+                # update the file name suffix:
+                suffix = '_' + str(counter)
+                counter += 1
+
+            # save accu if required:
+            outfile = pars.get('OUTFILE')
+            if outfile: 
+                datadir = pars.get('DATADIR')
+
+                # write raw data in Simpson format:
+                filename = datadir+outfile+suffix+'.dat'
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.dat')
+                accu.write_to_simpson(filename)
+
+                # write raw data in Tecmag format:
+                # filename = datadir+outfile+'.tnt'
+                # accu.write_to_tecmag(filename, nrecords=20)
+
+               # write parameters in a text file:
+                filename = datadir+outfile+suffix+'.par' 
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.par')
+
+                fileobject = open(filename, 'w')
+                for key in sorted(pars.iterkeys()):
+                    if key=='run': continue
+                    if key=='rec_phase': continue
+                    fileobject.write('%s%s%s'%(key,'=', pars[key]))
+                    fileobject.write('\n')
+                fileobject.close()
+
+            # reset accumulation:
+            del accu
+
+
+    if var_key == 'D2':
+        # KWW fit:
+        xdata = measurement.get_xdata()
+        ydata = measurement.get_ydata()
+        [amplitude, rate, beta] = fitfunc(xdata, ydata)
+        print '%s%02g' % ('Amplitude = ', amplitude)
+        print '%s%02g' % ('T2 [s] = ', 1./rate)
+        print '%s%02g' % ('Beta = ', beta)
+
+        # update display for the fit:
+        measurement.y = func([amplitude, rate, beta], xdata)
+        data[measurement.get_title()] = measurement
+
+# the fitting procedure:
+def fitfunc(xdata, ydata): 
+
+    # initialize variable parameters:
+    try: 
+        # solve Az = b:
+        A = array((ones(xdata.size/2), xdata[0:xdata.size/2]))    
+        b = log(abs(ydata[0:xdata.size/2]))     
+        z = linalg.lstsq(transpose(A), b)
+        amplitude = exp(z[0][0])
+        rate = -z[0][1]
+    except:
+        amplitude = abs(ydata[0])
+        rate = 1./(xdata[-1] - xdata[0])
+    beta = 1
+    p0 = [amplitude, rate, beta]
+
+    # run least-squares optimization:
+    plsq = leastsq(residuals, p0, args=(xdata, ydata))
+
+    return plsq[0] # best-fit parameters
+
+def residuals(p, xdata, ydata):
+    return ydata - func(p, xdata)
+
+# here is the function to fit:
+def func(p, xdata):
+    return p[0]*exp(-(p[1]*xdata)**p[2])
+
+
+pass

Scripts/Steady_Gradient_Spin_Echo/op_sgse_exp.py

@@ -0,0 +1,160 @@
+# -*- coding: iso-8859-1 -*-
+
+TXEnableDelay = 2e-6
+TXEnableValue = 0b0001  # TTL line blanking RF amplifier (bit 0)
+TXPulseValue  = 0b0010  # TTL line triggering RF pulses  (bit 1)
+ADCSensitivity = 2	# voltage span for ADC
+
+def experiment(): # the diffusion editing sequence with stimulated echo
+    
+  # set up acquisition parameters:
+    pars = {}
+    pars['P90'] = 1.7e-6	# 90-degree pulse length (s)
+    pars['SF'] = 338.7e6	# spectrometer frequency (Hz)
+    pars['O1'] = -60e3		# offset from SF (Hz)
+    pars['SW'] = 10e6		# spectral window (Hz)
+    pars['SI'] = 1*1024		# number of acquisition points
+    pars['NS'] = 8		# number of scans
+    pars['DS'] = 0		# number of dummy scans
+    pars['RD'] = 3		# delay between scans (s)
+    pars['D1'] = 20e-6		# delay after first pulse (short tau) (s)
+    pars['D2'] = 100e-6		# delay after second pulse (long tau) (s) 
+    pars['PHA'] = -150		# receiver phase (degree)    
+    pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+    pars['OUTFILE'] = None  			# output file name
+
+  # specify a variable parameter (optional):
+    pars['VAR_PAR'] = 'D2'	# variable parameter name (a string)   
+    start = 20e-6		# starting value
+    stop = 4e-0			# end value
+    steps = 24			# number of values
+    log_scale = True		# log scale flag 
+    stag_range = False		# staggered range flag
+
+  # check parameters for safety: 
+    if pars['PHA'] < 0:
+        pars['PHA'] = 360 + pars['PHA']
+
+    if pars['P90'] > 20e-6:
+        raise Exception("Pulse too long!!!")
+
+    # check whether a variable parameter is named: 
+    var_key = pars.get('VAR_PAR')
+    if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+        raise Exception("Pulse too long!!!")
+
+    if pars['NS']%16 != 0:
+        pars['NS'] = int(round(pars['NS'] / 16) + 1) * 16 
+        print 'Number of scans changed to ',pars['NS'],' due to phase cycling'
+
+    # start the experiment:
+    if var_key: 
+        # this is an arrayed experiment:
+        if log_scale:
+            array = log_range(start,stop,steps)
+        else:
+            array = lin_range(start,stop,steps) 
+
+        if stag_range:
+            array = staggered_range(array, size = 2)
+
+        # estimate the experiment time:
+        if var_key == 'D1':
+            seconds = (sum(array)*2 + (pars['D2'] + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'D2':
+            seconds = (sum(array) + (pars['D1']*2 + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'RD':
+            seconds = (sum(array) + (pars['D1']*2 + pars['D2']) * steps) * (pars['NS'] + pars['DS'])
+        else:
+            seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * steps * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for a variable parameter:
+        for index, pars[var_key] in enumerate(array):
+            print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\
+                             ' out of '+str(array.size)+': value = '+str(pars[var_key])
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield ste_experiment(pars, run) 
+            synchronize()
+    else:
+       # estimate the experiment time:
+        seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for accumulation:
+        for run in xrange(pars['NS']+pars['DS']):
+            yield ste_experiment(pars, run)       
+
+
+# the pulse program:
+
+def ste_experiment(pars, run):
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+        
+    pars['PROG'] = 'ste_experiment'
+
+    # phase lists (16-phase cycle from JMR 157, 31 (2002)):
+    pars['PH1'] = [0, 180,   0, 180,   0, 180,   0, 180,    90, 270,  90, 270,  90, 270,  90, 270] # 1st 90-degree pulse
+    pars['PH3'] = [0,   0, 180, 180,   0,   0, 180, 180,     0,   0, 180, 180,   0,   0, 180, 180] # 2nd 90-degree pulse
+    pars['PH4'] = [0,   0,   0,   0, 180, 180, 180, 180,     0,   0,   0,   0, 180, 180, 180, 180] # 3nd 90-degree pulse
+    pars['PH2'] = [0, 180, 180,   0, 180,   0,   0, 180,   270,  90,  90, 270,  90, 270, 270,  90] # receiver
+
+    # read in variables:
+    P90 = pars['P90']
+    SF =  pars['SF']
+    O1 =  pars['O1']    
+    RD =  pars['RD']
+    D1 =  pars['D1']
+    D2 =  pars['D2']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH4 = pars['PH4'][run%len(pars['PH4'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']
+
+    # set sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+
+   # run the pulse sequence: 
+    e.wait(RD)						# delay between scans
+    e.set_frequency(SF+O1, phase=PH1)
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse 
+
+    e.wait(D1-P90/2-TXEnableDelay)			# 'short tau'
+    e.set_phase(PH3)
+
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)		
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse 
+
+    e.wait(D2-P90/2-TXEnableDelay)			# 'long tau'
+    e.set_phase(PH4)
+
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)		
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse
+
+    e.set_phase(PHA)
+    e.wait(D1-P90/2-TXEnableDelay)			# 'short tau'
+    e.record(SI, SW, sensitivity=ADCSensitivity)	# acquisition
+
+
+    # write experiment parameters:  
+    for key in pars.keys():
+        e.set_description(key, pars[key])		# acquisition parameters
+    e.set_description('run', run)			# current scan
+    e.set_description('rec_phase', -PH2)		# current receiver phase
+
+    return e

Scripts/Steady_Gradient_Spin_Echo/op_sgse_res.py

@@ -0,0 +1,201 @@
+# -*- coding: iso-8859-1 -*-
+
+from numpy import *
+from scipy.signal import *
+from scipy.optimize import *
+from os import path, rename
+
+def result():
+
+    measurement = MeasurementResult('Magnetization')
+
+    measurement_range = [0.0, 10e-6]
+    measurement_ranging = True
+  
+    suffix = '' 	# output file name's suffix and...
+    counter = 1		# counter for arrayed experiments
+    var_key = ''	# variable parameter name
+
+    # loop over the incoming results: 
+    for timesignal in results:
+        if not isinstance(timesignal,ADC_Result): 
+            continue     
+
+        # read experiment parameters:
+        pars = timesignal.get_description_dictionary()
+
+        # ---------------- digital filter ------------------ 
+             
+        # get actual sampling rate of timesignal:
+        sampling_rate = timesignal.get_sampling_rate()
+
+        # get user-defined spectrum width:
+        spec_width = pars['SW']
+        
+        # specify cutoff frequency, in relative units:
+        cutoff = spec_width / sampling_rate
+                
+        if cutoff < 1: # otherwise no filter applied
+
+            # number of filter's coefficients:
+            numtaps = 29
+            
+            # use firwin to create a lowpass FIR filter:
+            fir_coeff = firwin(numtaps, cutoff)
+            
+            # downsize x according to user-defined spectral window:
+            skip = int(sampling_rate / spec_width)
+            timesignal.x = timesignal.x[::skip]
+            
+            for i in range(2):
+                # apply the filter to ith channel:
+                timesignal.y[i] = lfilter(fir_coeff, 1.0, timesignal.y[i])
+
+                # zeroize first N-1 "corrupted" samples:
+                timesignal.y[i][:numtaps-1] = 0.0
+
+                # circular left shift of y:
+                timesignal.y[i] = roll(timesignal.y[i], -(numtaps-1))    
+
+                # downsize y to user-defined number of samples (SI):
+                timesignal.y[i] = timesignal.y[i][::skip]
+
+            # update the sampling_rate attribute of the signal's:
+            timesignal.set_sampling_rate(spec_width)
+
+        # ----------------------------------------------------
+
+        # rotate timesignal according to current receiver's phase:
+        timesignal.phase(pars['rec_phase'])
+
+        # provide timesignal to the display tab:
+        data['Current scan'] = timesignal
+
+        # accumulate...
+        if not locals().get('accu'):
+            accu = Accumulation()
+
+        # skip dummy scans, if any:
+        if pars['run'] < 0: continue
+
+        # add up:
+        accu += timesignal
+
+        # provide accumulation to the display tab:
+        data['Accumulation'] = accu
+
+        # check how many scans are done:
+        if accu.n == pars['NS']: # accumulation is complete
+
+            # make a copy:
+            echo = accu + 0
+
+            # compute the signal's phase:
+            phi0 = arctan2(echo.y[1][0], echo.y[0][0]) * 180 / pi
+            if not locals().get('ref'): ref = phi0   
+            print 'phi0 = ', phi0
+
+            # rotate the signal to maximize Re (optional):
+            #echo.phase(-phi0) 
+            
+            # specify noise level: 
+            if not locals().get('noise'):
+                noise = 0.1*max(abs(echo.y[0]))
+                samples = echo.y[0] > noise
+
+            # check whether it is an arrayed experiment: 
+            var_key = pars.get('VAR_PAR')
+            if var_key:
+                # get variable parameter's value:
+                var_value = pars.get(var_key)
+                
+                # provide signal recorded with this var_value to the display tab:
+                data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu                
+
+                # measure signal intensity vs. var_value:
+                if measurement_ranging == True:
+                   [start, stop] = echo.get_sampling_rate() * array(measurement_range)
+                   measurement[var_value] = sum(echo.y[0][int(start):int(stop)]) 
+                   
+                else:
+                   measurement[var_value] = sum(echo.y[0][samples])
+
+                # provide the measurement result to the display tab:
+                data[measurement.get_title()] = measurement
+
+                # update the file name suffix:
+                suffix = '_' + str(counter)
+                counter += 1
+
+            # save accu if required:   
+            outfile = pars.get('OUTFILE')
+            if outfile: 
+                datadir = pars.get('DATADIR')
+
+                # write raw data in Simpson format:
+                filename = datadir+outfile+suffix+'.dat'
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.dat')
+                accu.write_to_simpson(filename)
+
+                # write raw data in Tecmag format:
+                # filename = datadir+outfile+'.tnt'
+                # accu.write_to_tecmag(filename, nrecords=20)
+
+               # write parameters in a text file:
+                filename = datadir+outfile+suffix+'.par' 
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.par')
+
+                fileobject = open(filename, 'w')
+                for key in sorted(pars.iterkeys()):
+                    if key=='run': continue
+                    if key=='rec_phase': continue
+                    fileobject.write('%s%s%s'%(key,'=', pars[key]))
+                    fileobject.write('\n')
+                fileobject.close()
+
+            # reset accumulation:
+            del accu
+
+    if var_key == 'D2':
+        # mono-exponential decay fit:
+        xdata = measurement.get_xdata()
+        ydata = measurement.get_ydata()
+        [amplitude, rate] = fitfunc(xdata, ydata)
+        print '%s%02g' % ('Amplitude = ', amplitude)
+        print '%s%02g' % ('T1 [s] = ', 1./rate)
+
+        # update display for the fit:
+        measurement.y = func([amplitude, rate], xdata)
+        data[measurement.get_title()] = measurement
+
+# the fitting procedure:
+def fitfunc(xdata, ydata): 
+
+    # initialize variable parameters:
+    try: 
+        # solve Az = b:
+        A = array((ones(xdata.size/2), xdata[0:xdata.size/2]))
+        b = log(abs(ydata[0:xdata.size/2]))
+        z = linalg.lstsq(transpose(A), b)
+        amplitude = exp(z[0][0])
+        rate = -z[0][1]
+    except:
+        amplitude = abs(ydata[0])
+        rate = 1./(xdata[-1] - xdata[0])
+    p0 = [amplitude, rate]
+
+    # run least-squares optimization:
+    plsq = leastsq(residuals, p0, args=(xdata, ydata))
+
+    return plsq[0] # best-fit parameters
+
+def residuals(p, xdata, ydata):
+    return ydata - func(p, xdata)
+
+# here is the function to fit:
+def func(p, xdata):
+    return p[0]*exp(-p[1]*xdata)
+
+pass

Scripts/Steady_Gradient_Spin_Echo_with_CPMG_Detection/op_sgse2_exp.py

@@ -0,0 +1,177 @@
+# -*- coding: iso-8859-1 -*-
+
+TXEnableDelay = 2e-6
+TXEnableValue = 0b0001  # TTL line blanking RF amplifier (bit 0)
+TXPulseValue  = 0b0010  # TTL line triggering RF pulses  (bit 1)
+ADCSensitivity = 2	# voltage span for ADC
+
+def experiment(): # the diffusion editing sequence with stimulated echo and CPMG detection
+
+  # set up acqusition parameters:
+    pars = {}
+    pars['P90'] = 1.7e-6	# 90-degree pulse length (s)
+    pars['SF'] = 338.7e6  	# spectrometer frequency (Hz)
+    pars['O1'] = -60e3		# offset from SF (Hz)
+    pars['NS'] = 16		# number of scans
+    pars['DS'] = 0		# number of dummy scans
+    pars['RD'] = 3		# delay between scans (s)
+    pars['D1'] = 20e-6		# delay after first STE pulse (s)
+    pars['D2'] = 100e-6		# delay after second STE pulse (s) 
+    pars['NECH'] = 16   	# number of 180-degree pulses in the CPMG train
+    pars['TAU'] = 40e-6         # half pulse period in the CPMG train (s)
+    pars['PHA'] = -127   	# receiver phase (degree) 
+    pars['DATADIR'] = '/home/fprak/Desktop/test/'	# data directory
+    pars['OUTFILE'] = None				# output file name
+
+  # specify a variable parameter (optional):
+    pars['VAR_PAR'] = 'D2'	# variable parameter name (a string)   
+    start = 20e-6		# starting value
+    stop = 4e-3			# end value
+    steps = 3			# number of values
+    log_scale = True		# log scale flag 
+    stag_range = False		# staggered range flag
+
+   # check parameters for safety: 
+    if pars['PHA'] < 0:
+        pars['PHA'] = 360 + pars['PHA']
+
+    if pars['P90'] > 20e-6:
+        raise Exception("Pulse too long!!!")
+
+    # check whether a variable parameter is named: 
+    var_key = pars.get('VAR_PAR')
+    if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+        raise Exception("Pulse too long!!!")
+
+    if pars['NS']%16 != 0:
+        pars['NS'] = int(round(pars['NS'] / 16) + 1) * 16 
+        print 'Number of scans changed to ',pars['NS'],' due to phase cycling'
+
+    # start the experiment:
+    if var_key: 
+        # this is an arrayed experiment:
+        if log_scale:
+            array = log_range(start,stop,steps)
+        else:
+            array = lin_range(start,stop,steps) 
+
+        if stag_range:
+            array = staggered_range(array, size = 2) 
+
+        # estimate the experiment time:
+        if var_key == 'D1':
+            seconds = (sum(array)*2 + (pars['D2'] + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'D2':
+            seconds = (sum(array) + (pars['D1']*2 + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'RD':
+            seconds = (sum(array) + (pars['D1']*2 + pars['D2']) * steps) * (pars['NS'] + pars['DS'])
+        else:
+            seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * steps * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+ 
+        # loop for a variable parameter:
+        for index, pars[var_key] in enumerate(array):
+            print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\
+                             ' out of '+str(array.size)+': value = '+str(pars[var_key])
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield ste2_experiment(pars, run)
+            synchronize() 
+    else:
+       # estimate the experiment time:
+        seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for accumulation:
+        for run in xrange(pars['NS']+pars['DS']):
+            yield ste2_experiment(pars, run)  
+
+
+# the pulse program: 
+
+def ste2_experiment(pars, run):
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+
+    pars['PROG'] = 'ste2_experiment'
+
+    # phase lists [from J. Magn. Reson. 157, 31 (2002)]: 
+    pars['PH1'] = [0, 180,   0, 180,   0, 180,   0, 180,  90, 270,  90, 270,  90, 270,  90, 270] # 1st 90-degree pulse
+    pars['PH3'] = [0,   0, 180, 180,   0,   0, 180, 180,   0,   0, 180, 180,   0,   0, 180, 180] # 2nd 90-degree pulse
+    pars['PH4'] = [0,   0,   0,   0, 180, 180, 180, 180,   0,   0,   0,   0, 180, 180, 180, 180] # 3nd 90-degree pulse
+    pars['PH5'] = [90, 90,  90,  90,  90,  90,  90,  90,   0,   0,   0,   0,   0,   0,   0,   0] # 180-degree pulses
+    pars['PH2'] = [0, 180, 180,   0, 180,   0,   0, 180, 270,  90,  90, 270,  90, 270, 270,  90] # receiver
+
+    # read in variables:
+    P90 = pars['P90']
+    P180 = pars['P90']*2
+    SF =  pars['SF']
+    O1 =  pars['O1']    
+    RD =  pars['RD']
+    D1 =  pars['D1']
+    D2 =  pars['D2']
+    TAU = pars['TAU']
+    NECH = pars['NECH']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH4 = pars['PH4'][run%len(pars['PH4'])]
+    PH5 = pars['PH5'][run%len(pars['PH5'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']  
+
+    # set sampling parameters:
+    SI = 128    	# number of echo samples
+    SW = 20e6   	# sample rate
+    AQ = (SI+6)/SW	# acquisition window
+	
+    if TAU < (P90+P180)/2+TXEnableDelay or TAU < (P180+TXEnableDelay+AQ)/2:
+        raise Exception('pulse period is too short!')
+
+    if 2*TAU < P180+TXEnableDelay+SI/SW:
+        raise Exception('pulse period too short!')
+    
+   # run the pulse sequence: 
+    e.wait(RD)						# delay between scans
+    e.set_frequency(SF+O1, phase=PH1)
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 1st 90-degree pulse 
+
+    e.wait(D1-P90/2-TXEnableDelay)			# short delay
+    e.set_phase(PH3)
+
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)		
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 2nd 90-degree pulse 
+
+    e.wait(D2-P90/2-TXEnableDelay)			# long delay
+    e.set_phase(PH4)
+
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)		
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 3rd 90-degree pulse
+    
+    e.wait(D1+TAU-P90/2-TXEnableDelay)			# wait for first echo and tau
+    e.set_phase(PH5)
+
+    e.loop_start(NECH)                          	# ----- loop for CPMG pulse train: -----
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)	# enable RF amplifier
+    e.ttl_pulse(P180, value=TXEnableValue|TXPulseValue)	# apply a 180-degree pulse
+    e.set_phase(PHA)                            	# set phase for receiver     
+    e.wait(TAU-(P180+TXEnableDelay+AQ)/2)       	# pre-acquisition delay
+    e.record(SI, SW, timelength=AQ, sensitivity=ADCSensitivity)	# acquire echo samples
+    e.wait(TAU-(P180+TXEnableDelay+AQ)/2)       	# post-acquisition delay
+    e.set_phase(PH5)                            	# set phase for theta-degree pulse
+    e.loop_end()                                	# --------------------------------------
+
+    # write experiment attributes:  
+    for key in pars.keys():
+        e.set_description(key, pars[key])		# acqusition parameters
+    e.set_description('run', run)			# current scan
+    e.set_description('rec_phase', -PH2)		# current receiver phase
+
+    return e

Scripts/Steady_Gradient_Spin_Echo_with_CPMG_Detection/op_sgse2_res.py

@@ -0,0 +1,185 @@
+# -*- coding: iso-8859-1 -*-
+
+from numpy import *
+from scipy.signal import *
+from scipy.optimize import *
+from os import path, rename
+
+def result():
+
+    measurement = MeasurementResult('Magnetization')
+
+    suffix = '' 	# output file name's suffix and...
+    counter = 1		# counter for arrayed experiments
+    var_key = ''	# variable parameter name
+
+    # loop over the incoming results: 
+    for timesignal in results:
+        if not isinstance(timesignal,ADC_Result): 
+            continue   
+
+        # read experiment parameters:
+        pars = timesignal.get_description_dictionary()
+
+        # rotate timesignal acording to current receiver's phase:
+        timesignal.phase(pars['rec_phase'])
+
+        # provide timesignal to the display tab:
+        data['Current scan'] = timesignal
+
+        # accumulate...
+        if not locals().get('accu'):
+            accu = Accumulation()
+
+        # skip dummy scans, if any:
+        if pars['run'] < 0: continue
+        
+        # add up:
+        accu += timesignal
+                
+        # provide accumulation to the display tab:
+        data['Accumulation'] = accu
+
+        # check how many scans are done:
+        if accu.n == pars['NS']: # accumulation is complete
+
+            # get number of echoes:
+            num_echoes = pars['NECH'] 
+
+            # downsize accu to one point per echo:
+            echodecay = accu + 0
+            echodecay.x = resize(echodecay.x, int(num_echoes))
+            echodecay.y[0] = resize(echodecay.y[0], int(num_echoes))
+            echodecay.y[1] = resize(echodecay.y[1], int(num_echoes))
+            
+            # find where first echo exceeds noise:
+#            if not locals().get('noise'):
+#                last_echo = accu.get_accu_by_index(num_echoes-1)
+#                first_echo = accu.get_accu_by_index(0)
+#                noise = 3*std(last_echo.y[1][:])
+#                samples = first_echo.y[0] > noise
+                
+            if not locals().get('noise'):
+                echo = accu.get_accu_by_index(0)
+                noise = 0.1*max(abs(echo.y[0]))
+                samples = abs(echo.y[0]) > noise
+
+            # set echo times and intensities:
+            for i in range(num_echoes):
+                # get ith echo:
+                echo = accu.get_accu_by_index(i)
+                # set echo time:
+                echodecay.x[i] = i*2*pars['TAU']
+                # set echo value:
+                echodecay.y[0][i] = sum(echo.y[0][samples])	# the sum of echo points that exceed noise...
+                echodecay.y[1][i] = sum(echo.y[1][samples]) 
+                #echodecay.y[0][i] = sum(echo.y[0])             # or the sum of all echo points...
+                #echodecay.y[1][i] = sum(echo.y[1]) 
+                #echodecay.y[0][i] = echo.y[0][echo.x.size/2]	# or a middle echo point
+                #echodecay.y[1][i] = echo.y[1][echo.x.size/2]
+                
+            # compute a signal's phase:
+            phi0 = arctan2(echodecay.y[1][0], echodecay.y[0][0]) * 180 / pi
+            if not locals().get('ref'):  ref = phi0 
+            print 'phi0 = ', phi0 
+                     
+            # rotate signal to maximize Re (optional): 
+            #echodecay.phase(-phi0)
+
+            # provide echo decay to the display tab:
+            data['Echo Decay'] = echodecay
+
+            # check whether it is an arrayed experiment: 
+            var_key = pars.get('VAR_PAR')
+            if var_key:
+                # get variable parameter's value:
+                var_value = pars.get(var_key)
+
+                # provide signal recorded with different var_value's to the display tab:
+                data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu       
+                data['Echo Decay'+"/"+var_key+"=%e"%(var_value)] = echodecay
+
+                # measure signal intensity vs. var_value:
+#                [amplitude, rate] = fitfunc(echodecay.x, echodecay.y[0])
+#                print '%s%02g' % ('Amplitude = ', amplitude)
+#                print '%s%02g' % ('T2 [s] = ', 1./rate)
+
+#                measurand[var_value] = amplitude
+                measurement[var_value] = sum(echodecay.y[0][:])
+
+                # provide measurement to the display tab:
+                data[measurement.get_title()] = measurement
+
+            # save accu if required:   
+            outfile = pars.get('OUTFILE')
+            if outfile: 
+                datadir = pars.get('DATADIR')
+
+                # write raw data in Simpson format:
+                filename = datadir+outfile+suffix+'.dat'
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.dat')
+                accu.write_to_simpson(filename)
+
+                # write raw data in Tecmag format:
+                # filename = datadir+outfile+'.tnt'
+                # accu.write_to_tecmag(filename, nrecords=20)
+
+               # write parameters in a text file:
+                filename = datadir+outfile+suffix+'.par' 
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.par')
+
+                fileobject = open(filename, 'w')
+                for key in sorted(pars.iterkeys()):
+                    if key=='run': continue
+                    if key=='rec_phase': continue
+                    fileobject.write('%s%s%s'%(key,'=', pars[key]))
+                    fileobject.write('\n')
+                fileobject.close()
+
+            # reset accumulation:
+            del accu
+
+    if var_key == 'D2':
+        # mono-exponential decay fit:
+        xdata = measurement.get_xdata()
+        ydata = measurement.get_ydata()
+        [amplitude, rate] = fitfunc(xdata, ydata)
+        print '%s%02g' % ('Amplitude = ', amplitude)
+        print '%s%02g' % ('T2 [s] = ', 1./rate)
+
+        # update display for the fit:
+        measurement.y = func([amplitude, rate], xdata)
+        data[measurement.get_title()] = measurement
+
+# the fitting procedure:
+def fitfunc(xdata, ydata):
+
+    # initialize variable parameters:
+    try: 
+        # solve Az = b:
+        A = array((ones(xdata.size/2), xdata[0:xdata.size/2]))
+        b = log(abs(ydata[0:xdata.size/2]))
+        z = linalg.lstsq(transpose(A), b)
+        amplitude = exp(z[0][0])
+        rate = -z[0][1]
+        print amplitude, 1./rate
+    except:
+        amplitude = abs(ydata[0])
+        rate = 1./(xdata[-1] - xdata[0])
+    p0 = [amplitude, rate]
+
+    # run least-squares optimization:
+    plsq = leastsq(residuals, p0, args=(xdata, ydata))
+
+    return plsq[0] # best-fit parameters
+
+def residuals(p, xdata, ydata):
+    return ydata - func(p, xdata)
+
+# here is the function to fit:
+def func(p, xdata):
+    return p[0]*exp(-p[1]*xdata)
+
+pass

Scripts/Stimulated_Echo/op_ste_exp.py

@@ -0,0 +1,165 @@
+# -*- coding: iso-8859-1 -*-
+
+TXEnableDelay = 2e-6
+TXEnableValue = 0b0001	# TTL line enabling RF amplifier (bit 0)
+TXPulseValue  = 0b0010	# TTL line triggering RF pulses  (bit 1)
+ADCSensitivity = 2	# voltage span for ADC
+
+def experiment(): # stimulated echo experiment
+
+  # set up acquisition parameters:
+    pars = {}
+    pars['P90'] = 1.7e-6	# 90-degree pulse length (s)
+    pars['SF'] = 338.7e6	# spectrometer frequency (Hz)
+    pars['O1'] = -87e3		# offset from SF (Hz)
+    pars['SW'] = 500e3		# spectrum width (Hz)
+    pars['SI'] = 1*1024		# number of acquisition points
+    pars['NS'] = 16		# number of scans
+    pars['DS'] = 0		# number of dummy scans
+    pars['RD'] = 1		# delay between scans (s)
+    pars['D1'] = 100e-6		# delay after first pulse (short tau) (s)
+    pars['D2'] = 100e-6		# delay after second pulse (long tau) (s) 
+    pars['D4'] = 0e-6		# echo pre-acquisition delay (s)
+    pars['PHA'] = 30		# receiver phase (degree) 
+    pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+    pars['OUTFILE'] = None			# output file name
+
+  # specify a variable parameter (optional):
+    pars['VAR_PAR'] = None	# variable parameter name (a string)   
+    start = 50e-6		# starting value
+    stop = 3e-3			# end value
+    steps = 10			# number of values
+    log_scale = False		# log scale flag 
+    stag_range = False		# staggered range flag
+
+  # check parameters for safety: 
+    if pars['PHA'] < 0:
+        pars['PHA'] = 360 + pars['PHA']
+
+    if pars['P90'] > 20e-6:
+        raise Exception("Pulse too long!!!")
+
+    # check whether a variable parameter is named: 
+    var_key = pars.get('VAR_PAR')
+    if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+        raise Exception("Pulse too long!!!")
+
+    if pars['NS']%16 != 0:
+        pars['NS'] = int(round(pars['NS'] / 16) + 1) * 16 
+        print 'Number of scans changed to ',pars['NS'],' due to phase cycling'
+
+   # start the experiment:
+
+    # check if a variable parameter is named: 
+    var_key = pars.get('VAR_PAR')
+    if var_key: 
+        # this is an arrayed experiment:
+        if log_scale:
+            array = log_range(start,stop,steps)
+        else:
+            array = lin_range(start,stop,steps)
+
+        if stag_range:
+            array = staggered_range(array, size = 2)
+
+        # estimate the experiment time:
+        if var_key == 'D1':
+            seconds = (sum(array)*2 + (pars['D2'] + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'D2':
+            seconds = (sum(array) + (pars['D1']*2 + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'RD':
+            seconds = (sum(array) + (pars['D1']*2 + pars['D2']) * steps) * (pars['NS'] + pars['DS'])
+        else:
+            seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * steps * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for a variable parameter:
+        for index, pars[var_key] in enumerate(array):
+            print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\
+                             ' out of '+str(array.size)+': value = '+str(pars[var_key])
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield ste_experiment(pars, run) 
+            synchronize()
+    else:
+        # estimate the experiment time:
+        seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for accumulation:
+        for run in xrange(pars['NS']+pars['DS']):
+            yield ste_experiment(pars, run)
+
+
+# the pulse program:
+
+def ste_experiment(pars, run):
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+
+    pars['PROG'] = 'ste_experiment'
+
+    # phase lists [16-phase cycle from JMR 157, 31 (2002)]:
+    pars['PH1'] = [0, 180,   0, 180,   0, 180,   0, 180,  90, 270,  90, 270,  90, 270,  90, 270] # 1st 90-degree pulse
+    pars['PH3'] = [0,   0, 180, 180,   0,   0, 180, 180,   0,   0, 180, 180,   0,   0, 180, 180] # 2nd 90-degree pulse
+    pars['PH4'] = [0,   0,   0,   0, 180, 180, 180, 180,   0,   0,   0,   0, 180, 180, 180, 180] # 3nd 90-degree pulse
+    pars['PH2'] = [0, 180, 180,   0, 180,   0,   0, 180, 270,  90,  90, 270,  90, 270, 270,  90] # receiver
+
+    # read in variables:
+    P90 = pars['P90']
+    SF =  pars['SF']
+    O1 =  pars['O1']
+    RD =  pars['RD']
+    D1 =  pars['D1']
+    D2 =  pars['D2']
+    D4 =  pars['D4']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH4 = pars['PH4'][run%len(pars['PH4'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']  
+
+    # set sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+
+   # run the pulse sequence: 
+    e.wait(RD)						# delay between scans
+    e.set_frequency(SF+O1, phase=PH1)
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse 
+
+    e.wait(D1-P90/2-TXEnableDelay)			# 'short tau'
+    e.set_phase(PH3)
+
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)		
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse 
+
+    e.wait(D2-P90/2-TXEnableDelay)			# 'long tau'
+    e.set_phase(PH4)
+
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)		
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse
+
+    e.set_phase(PHA)
+    e.wait(D1-P90/2-TXEnableDelay+D4)      		# 'short tau'
+    e.record(SI, SW, sensitivity=ADCSensitivity)	# acquisition
+
+
+    # write experiment parameters:  
+    for key in pars.keys():
+        e.set_description(key, pars[key])		# acquisition parameters
+    e.set_description('run', run)			# current scan
+    e.set_description('rec_phase', -PH2)		# current receiver phase
+
+    return e

Scripts/Stimulated_Echo/op_ste_res.py

@@ -0,0 +1,209 @@
+# -*- coding: iso-8859-1 -*-
+
+from numpy import *
+from scipy.signal import *
+from scipy.optimize import *
+from os import path, rename
+
+def result():
+
+    measurement = MeasurementResult('Magnetization')
+
+    measurement_range = [0.0, 10e-6]
+    measurement_ranging = False
+
+    suffix = '' 	# output file name's suffix and...
+    counter = 1		# counter for arrayed experiments
+    var_key = ''	# variable parameter name
+
+    # loop over the incoming results: 
+    for timesignal in results:
+        if not isinstance(timesignal,ADC_Result): 
+            continue
+
+        # read experiment parameters:
+        pars = timesignal.get_description_dictionary()
+
+        # ---------------- digital filter ------------------ 
+             
+        # get actual sampling rate of timesignal:
+        sampling_rate = timesignal.get_sampling_rate()
+
+        # get user-defined spectrum width:
+        spec_width = pars['SW']
+        
+        # specify cutoff frequency, in relative units:
+        cutoff = spec_width / sampling_rate
+                
+        if cutoff < 1: # otherwise no filter applied
+
+            # number of filter's coefficients:
+            numtaps = 29
+            
+            # use firwin to create a lowpass FIR filter:
+            fir_coeff = firwin(numtaps, cutoff)
+            
+            # downsize x according to user-defined spectral window:
+            skip = int(sampling_rate / spec_width)
+            timesignal.x = timesignal.x[::skip]
+            
+            for i in range(2):
+                # apply the filter to ith channel:
+                timesignal.y[i] = lfilter(fir_coeff, 1.0, timesignal.y[i])
+
+                # zeroize first N-1 "corrupted" samples:
+                timesignal.y[i][:numtaps-1] = 0.0
+
+                # circular left shift of y:
+                timesignal.y[i] = roll(timesignal.y[i], -(numtaps-1))    
+
+                # downsize y to user-defined number of samples (SI):
+                timesignal.y[i] = timesignal.y[i][::skip]
+
+            # update the sampling_rate attribute of the signal's:
+            timesignal.set_sampling_rate(spec_width)
+
+        # ----------------------------------------------------
+
+        # phase timesignal according to current rec_phase:
+        timesignal.phase(pars['rec_phase'])
+
+        # provide timesignal to the display tab:
+        data['Current scan'] = timesignal
+
+        # accumulate...
+        if not locals().get('accu'):
+            accu = Accumulation()
+
+        # skip dummy scans, if any:
+        if pars['run'] < 0: continue
+
+        # add up:
+        accu += timesignal
+
+        # provide accumulation to the display tab:
+        data['Accumulation'] = accu
+
+        # check how many scans are done:
+        if accu.n == pars['NS']: # accumulation is complete
+
+            # make a copy:
+            echo = accu + 0
+
+            # compute the initial phase of the signal:
+            phi0 = arctan2(accu.y[1][0], accu.y[0][0]) * 180 / pi
+            if not locals().get('ref'): ref = phi0   
+            print 'phi0 = ', phi0
+
+            # rotate the signal to maximize Re (optional):     
+            #echo.phase(-phi0) 
+
+            # do FFT: 
+            echo.exp_window(line_broadening=10)
+            spectrum = echo.fft(samples=2*pars['SI'])
+
+            # try zero-order phase correction:
+            spectrum.phase(-phi0)
+
+            # provide spectrum to the display tab:
+            data['Spectrum'] = spectrum
+
+            # check whether it is an arrayed experiment:
+            var_key = pars.get('VAR_PAR')
+            if var_key:
+                # get variable parameter's value:
+                var_value = pars.get(var_key)
+
+                # provide signal recorded with var_value to the display tab:
+                data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu                
+
+                # measure signal intensity vs. var_value:
+                if measurement_ranging == True:
+                   [start, stop] = accu.get_sampling_rate() * array(measurement_range)
+                   measurement[var_value] = sum(accu.y[0][int(start):int(stop)]) 
+                   
+                else:
+                   measurement[var_value] = sum(accu.y[0][0:31])
+
+                # provide measurement to the display tab:
+                data[measurement.get_title()] = measurement
+
+                # update the file name suffix:
+                suffix = '_' + str(counter)
+                counter += 1
+
+            # save accu if required:
+            outfile = pars.get('OUTFILE')
+            if outfile:
+                datadir = pars.get('DATADIR')
+
+                # write raw data in Simpson format:
+                filename = datadir+outfile+suffix+'.dat'
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.dat')
+                accu.write_to_simpson(filename)
+
+                # write raw data in Tecmag format:
+                # filename = datadir+outfile+'.tnt'
+                # accu.write_to_tecmag(filename, nrecords=20)
+
+                # write parameters in a text file:
+                filename = datadir+outfile+suffix+'.par' 
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.par')
+
+                fileobject = open(filename, 'w')
+                for key in sorted(pars.iterkeys()):
+                    if key=='run': continue
+                    if key=='rec_phase': continue
+                    fileobject.write('%s%s%s'%(key,'=', pars[key]))
+                    fileobject.write('\n')
+                fileobject.close()
+
+            # reset accumulation:
+            del accu
+
+    if var_key == 'D2':
+        # mono-exponential decay fit:
+        xdata = measurement.get_xdata()
+        ydata = measurement.get_ydata()
+        [amplitude, rate, offset] = fitfunc(xdata, ydata)
+        print '%s%02g' % ('Amplitude = ', amplitude)
+        print '%s%02g' % ('T1 [s] = ', 1./rate)
+        print '%s%02g' % ('Offset = ', offset)
+
+        # update display for the fit:
+        measurement.y = func([amplitude, rate, offset], xdata)
+        data[measurement.get_title()] = measurement
+
+# the fitting procedure:
+def fitfunc(xdata, ydata):
+
+    # initialize variable parameters:
+    try: 
+        # solve Az = b:
+        A = array((ones(xdata.size/2), xdata[0:xdata.size/2]))
+        b = log(abs(ydata[0:xdata.size/2]))
+        z = linalg.lstsq(transpose(A), b)
+        amplitude = exp(z[0][0])
+        rate = -z[0][1]
+    except:
+        amplitude = abs(ydata[0])
+        rate = 1./(xdata[-1] - xdata[0])
+    offset = min(ydata)
+    p0 = [amplitude, rate, offset]
+
+    # run least-squares optimization:
+    plsq = leastsq(residuals, p0, args=(xdata, ydata))
+
+    return plsq[0] # best-fit parameters
+
+def residuals(p, xdata, ydata):
+    return ydata - func(p, xdata)
+
+# here is the function to fit:
+def func(p, xdata):
+    return p[0]*exp(-p[1]*xdata) + p[2]
+
+
+pass

Scripts/T1Q/op_t1q_exp.py

@@ -0,0 +1,169 @@
+# -*- coding: iso-8859-1 -*-
+
+TXEnableDelay = 2e-6
+TXEnableValue = 0b0001  # TTL line blanking RF amplifier (bit 0)
+TXPulseValue  = 0b0010  # TTL line triggering RF pulses  (bit 1)
+ADCSensitivity = 2	# voltage span for ADC
+
+def experiment(): # Jeener-Broekaert sequence with solid-echo detection to measure T1Q [JMR 43, 213 (1981)].
+
+  # set up acquisition parameters:
+    pars = {}
+    pars['P90'] = 2.3e-6	# 90-degree pulse length (s)
+    pars['SF'] = 46.7e6 	# spectrometer frequency (Hz)
+    pars['O1'] = 5.6e3		# offset from SF (Hz)
+    pars['SW'] = 200e3		# spectrum width (Hz)
+    pars['SI'] = 1*256		# number of acquisition points
+    pars['NS'] = 32		# number of scans
+    pars['DS'] = 0		# number of dummy scans
+    pars['RD'] = .5		# delay between scans (s)
+    pars['D1'] = 30e-6		# delay after first pulse or 'short tau' (s)
+    pars['D2'] = 30e-6		# delay after second pulse or 'long tau' (s) 
+    pars['PHA'] = -27		# receiver phase (degree)
+    pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+    pars['OUTFILE'] = None			# output file name
+
+  # specify a variable parameter (optional):
+    pars['VAR_PAR'] = 'D2'	# variable parameter name (a string)   
+    start = 30e-6		# starting value
+    stop = 2e-0			# end value
+    steps = 16			# number of values
+    log_scale = True		# log scale flag 
+    stag_range = False		# staggered range flag
+
+  # check parameters for safety: 
+    if pars['PHA'] < 0:
+        pars['PHA'] = 360 + pars['PHA']
+
+    if pars['P90'] > 20e-6:
+        raise Exception("Pulse too long!!!")
+
+    # check whether a variable parameter is named: 
+    var_key = pars.get('VAR_PAR')
+    if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+        raise Exception("Pulse too long!!!")
+
+    if pars['NS']%8 != 0:
+        pars['NS'] = int(round(pars['NS'] / 8) + 1) * 8 
+        print 'Number of scans changed to ',pars['NS'],' due to phase cycling'
+
+    # start the experiment:
+    if var_key:
+        # this is an arrayed experiment:
+        if log_scale:
+            array = log_range(start,stop,steps)
+        else:
+            array = lin_range(start,stop,steps)
+
+        if stag_range:
+            array = staggered_range(array, size = 2)
+
+        # estimate the experiment time:
+        if var_key == 'D1':
+            seconds = (sum(array)*2 + (pars['D2'] + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'D2':
+            seconds = (sum(array) + (pars['D1']*2 + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'RD':
+            seconds = (sum(array) + (pars['D1']*2 + pars['D2']) * steps) * (pars['NS'] + pars['DS'])
+        else:
+            seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * steps * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for a variable parameter:
+        for index, pars[var_key] in enumerate(array):
+            print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\
+                             ' out of '+str(array.size)+': value = '+str(pars[var_key])
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield t1q_experiment(pars, run) 
+            synchronize()
+    else:
+       # estimate the experiment time:
+        seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for accumulation:
+        for run in xrange(pars['NS']+pars['DS']):
+            yield t1q_experiment(pars, run)
+
+
+# the pulse program:
+
+def t1q_experiment(pars, run):
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+
+    pars['PROG'] = 't1q_experiment'
+
+    # phase cycle from J. Magn. Reson. 43, 213-223 (1981) extended for 90-deg refocusing pulse which also eliminates echo(es) after 3rd 45-deg pulse: 
+    pars['PH1'] = [0,   0,  90,  90,     180, 180, 270, 270] # 90-deg pulse
+    pars['PH3'] = [90, 90, 180, 180,     270, 270,   0,   0] # 1st 45-deg pulse
+    pars['PH4'] = [0, 180,   0, 180,     180,   0, 180,   0] # 2nd 45-deg pulse
+    pars['PH5'] = [0,   0, 180, 180,     270,  90,  90, 270] # refocucing 90-deg pulse
+    pars['PH2'] = [0, 180,   0, 180,     180,   0, 180,   0] # receiver
+
+
+    # read in variables:
+    P90 = pars['P90']
+    P45 = pars['P90']*0.5
+    SF =  pars['SF']
+    O1 =  pars['O1']
+    RD =  pars['RD']
+    D1 =  pars['D1']
+    D2 =  pars['D2']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH4 = pars['PH4'][run%len(pars['PH4'])]
+    PH5 = pars['PH5'][run%len(pars['PH5'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']
+
+    # set sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+
+   # run the pulse sequence: 
+    e.wait(RD)   					# delay between scans
+    e.set_frequency(SF+O1, phase=PH1)           				 
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)		
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse 
+
+    e.wait(D1-P90/2-TXEnableDelay)			# 'short tau'
+    e.set_phase(PH3)
+
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)		
+    e.ttl_pulse(P45, value=TXEnableValue|TXPulseValue)	# 45-degree pulse 
+
+    e.wait(D2-P45/2-TXEnableDelay)			# 'long tau'
+    e.set_phase(PH4)
+
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)		
+    e.ttl_pulse(P45, value=TXEnableValue|TXPulseValue)	# 45-degree pulse
+
+    e.wait(10e-6-TXEnableDelay)
+    e.set_phase(PH5)
+
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)		
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse
+
+    e.set_phase(PHA)
+    e.wait(13e-6)       	    	               		
+    e.record(SI, SW, sensitivity=ADCSensitivity)	# acquisition                                     
+
+   # write experiment attributes:  
+    for key in pars.keys():
+        e.set_description(key, pars[key])		# acqusition parameters
+    e.set_description('run', run)			# current scan
+    e.set_description('rec_phase', -PH2)		# current receiver phase
+
+    return e

Scripts/T1Q/op_t1q_res.py

@@ -0,0 +1,209 @@
+# -*- coding: iso-8859-1 -*-
+
+from numpy import *
+from scipy.signal import *
+from scipy.optimize import *
+from scipy.fftpack import rfft
+from os import path, rename
+
+def result(): 
+
+    measurement = MeasurementResult('Magnetization')
+
+    suffix = '' 	# output file name's suffix and...
+    counter = 1		# counter for arrayed experiments
+    var_key = ''	# variable parameter name
+
+    # loop over the incoming results: 
+    for timesignal in results:
+        if not isinstance(timesignal,ADC_Result): 
+            continue
+
+        # read experiment parameters:
+        pars = timesignal.get_description_dictionary()
+
+        # ---------------- digital filter ------------------ 
+             
+        # get actual sampling rate of timesignal:
+        sampling_rate = timesignal.get_sampling_rate()
+
+        # get user-defined spectrum width:
+        spec_width = pars['SW']
+        
+        # specify cutoff frequency, in relative units:
+        cutoff = spec_width / sampling_rate
+                
+        if cutoff < 1: # otherwise no filter applied
+
+            # number of filter's coefficients:
+            numtaps = 29
+            
+            # use firwin to create a lowpass FIR filter:
+            fir_coeff = firwin(numtaps, cutoff)
+            
+            # downsize x according to user-defined spectral window:
+            skip = int(sampling_rate / spec_width)
+            timesignal.x = timesignal.x[::skip]
+            
+            for i in range(2):
+                # apply the filter to ith channel:
+                timesignal.y[i] = lfilter(fir_coeff, 1.0, timesignal.y[i])
+
+                # zeroize first N-1 "corrupted" samples:
+                timesignal.y[i][:numtaps-1] = 0.0
+
+                # circular left shift of y:
+                timesignal.y[i] = roll(timesignal.y[i], -(numtaps-1))    
+
+                # downsize y to user-defined number of samples (SI):
+                timesignal.y[i] = timesignal.y[i][::skip]
+
+            # update the sampling_rate attribute of the signal's:
+            timesignal.set_sampling_rate(spec_width)
+
+        # ----------------------------------------------------
+
+        # phase timesignal according to current rec_phase:
+        timesignal.phase(pars['rec_phase'])
+
+        # provide timesignal to the display tab:
+        data['Current scan'] = timesignal
+
+        # accumulate...
+        if not locals().get('accu'):
+            accu = Accumulation()
+
+        # skip dummy scans, if any:
+        if pars['run'] < 0: continue
+
+        # add up:
+        accu += timesignal
+
+        # provide accumulation to the display tab:
+        data['Accumulation'] = accu
+
+        # check how many scans are done:
+        if accu.n == pars['NS']: # accumulation is complete
+
+            # make a copy:
+            fid = accu + 0
+
+            # compute the signal's phase:
+            phi0 = arctan2(fid.y[1][0], fid.y[0][0]) * 180 / pi
+            if not 'ref' in locals(): ref = phi0
+            print 'phi0 = ', phi0
+
+            # do FFT: 
+            fid.exp_window(line_broadening=10)
+            spectrum = fid.fft(samples=2*pars['SI'])
+            spectrum.baseline()
+
+            # provide spectrum to the display tab:
+            data['Spectrum'] = spectrum
+
+            # check whether it is an arrayed experiment:
+            var_key = pars.get	('VAR_PAR')
+            if var_key:
+                # get variable parameter's value:
+                var_value = pars.get(var_key)
+                
+                # provide signal recorded with the var_value to the display tab:
+                data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu 
+
+                # measure signal intensity vs. var_value:
+                signal = (accu + 0).y[1]
+                
+                # -*- discrete cosine transform of Im -*-
+                N = len(signal)
+                y = empty(2*N, float)
+                y[:N] = signal[:]
+                y[N:] = signal[::-1]
+                c = rfft(y)
+                phi = exp(-1j*pi*arange(N)/(2*N))
+                dct = real(phi*c[:N])
+                # ---------------------------------------
+                measurement[var_value] = sum(dct[0:9])
+
+                # provide measurement to the display tab:
+                data[measurement.get_title()] = measurement
+
+                # update the file name suffix:
+                suffix = '_' + str(counter)
+                counter += 1
+
+            # save accu if required:
+            outfile = pars.get('OUTFILE')
+            if outfile: 
+                datadir = pars.get('DATADIR')
+
+                # write raw data in Simpson format:
+                filename = datadir+outfile+suffix+'.dat'
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.dat')
+                accu.write_to_simpson(filename)
+
+                # write raw data in Tecmag format:
+                # filename = datadir+outfile+'.tnt'
+                # accu.write_to_tecmag(filename, nrecords=20)
+
+               # write parameters in a text file:
+                filename = datadir+outfile+suffix+'.par' 
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.par')
+
+                fileobject = open(filename, 'w')
+                for key in sorted(pars.iterkeys()):
+                    if key=='run': continue
+                    if key=='rec_phase': continue
+                    fileobject.write('%s%s%s'%(key,'=', pars[key]))
+                    fileobject.write('\n')
+                fileobject.close()
+
+            # reset accumulation:
+            del accu
+
+    if var_key == 'D2':
+        # mono-exponential decay fit:
+        xdata = measurement.get_xdata()
+        ydata = measurement.get_ydata()
+
+        [amplitude, rate, offset] = fitfunc(xdata, ydata)
+        print '%s%02g' % ('Amplitude = ', amplitude)
+        print '%s%02g' % ('T1Q [s] = ', 1./rate)
+
+        # update display for the fit:
+        measurement.y = func([amplitude, rate, offset], xdata)
+        data[measurement.get_title()] = measurement
+
+
+# the fitting procedure:
+def fitfunc(xdata, ydata): 
+
+    # initialize variable parameters:
+    try: 
+        # solve Az = b:
+        A = array((ones(xdata.size/2), xdata[0:xdata.size/2]))    
+        b = log(abs(ydata[0:xdata.size/2]))     
+        z = linalg.lstsq(transpose(A), b)
+        amplitude = exp(z[0][0])
+        rate = -z[0][1]
+    except:
+        amplitude = abs(ydata[0])
+        rate = 1./(xdata[-1] - xdata[0])
+    offset = min(ydata)        
+    p0 = [amplitude, rate, offset]
+
+    # run least-squares optimization:
+    plsq = leastsq(residuals, p0, args=(xdata, ydata))
+
+    return plsq[0] # best-fit parameters
+
+def residuals(p, xdata, ydata):
+    return ydata - func(p, xdata)
+
+# here is the function to fit:
+def func(p, xdata):
+    return p[0]*exp(-p[1]*xdata)+p[2]
+
+
+pass

Scripts/ZZ_Exchange_with_T2_Selection/op_t2zz_exp.py

@@ -0,0 +1,172 @@
+# -*- coding: iso-8859-1 -*-
+
+TXEnableDelay = 2e-6
+TXEnableValue = 0b0001  # TTL line blanking RF amplifier (bit 0)
+TXPulseValue  = 0b0010  # TTL line triggering RF pulses  (bit 1)
+ADCSensitivity = 2	# voltage span for ADC
+
+def experiment(): # ZZ-exchange with T2-selection  
+    
+  # set up acqusition parameters:      
+    pars = {}
+    pars['P90'] = 1.7e-6	# 90-degree pulse length (s)
+    pars['SF'] = 338.7e6	# spectrometer frequency (Hz)
+    pars['O1'] = -60e3  	# offset from SF (Hz)
+    pars['SW'] = 200e3 		# spectral width (Hz)
+    pars['SI'] = 128	 	# acquisition points
+    pars['D1'] = 10e-6   	# echo delay (s)
+    pars['D2'] = 100e-6         # z-storage duration (s)
+    pars['NS'] = 16      	# number of scans
+    pars['DS'] = 0      	# number of dummy scans
+    pars['RD'] = 3      	# delay between scans (s)
+    pars['DEAD1'] = 5e-6	# receiver dead time (s)  
+    pars['PHA'] = 0     	# receiver phase (degree)
+    pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+    pars['OUTFILE'] = None                            	# output file name
+    
+  # specify a variable parameter:
+    pars['VAR_PAR'] = 'D2'      	# variable parameter name
+    start = 20e-6          	# first value
+    stop  = 1e-3      	 	# last value
+    steps = 20  	      	# number of values
+    log_scale = True		# log scale flag 
+    stag_range = False		# staggered range flag
+    
+  # check parameters for safety:     
+    if pars['PHA'] < 0:
+        pars['PHA'] = 360 + pars['PHA']
+        
+    if pars['P90'] > 20e-6:
+        raise Exception("Pulse too long!!!")                    
+
+    if pars['NS']%16 != 0:
+        pars['NS'] = (pars['NS'] / 16 + 1) * 16 
+        print 'Number of scans changed to', pars['NS'], 'due to phase cycling'
+
+    # check whether a variable parameter is named: 
+    var_key = pars.get('VAR_PAR')
+    if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+        raise Exception("Pulse too long!!!")
+
+    if pars['NS']%8 != 0:
+        pars['NS'] = int(round(pars['NS'] / 8) + 1) * 8 
+        print 'Number of scans changed to ',pars['NS'],' due to phase cycling'
+
+    # start the experiment:
+    if var_key:
+        # this is an arrayed experiment:
+        if log_scale:
+            array = log_range(start,stop,steps)
+        else:
+            array = lin_range(start,stop,steps)
+
+        if stag_range:
+            array = staggered_range(array, size = 2)
+
+        # estimate the experiment time:
+        if var_key == 'D1':
+            seconds = (sum(array)*2 + (pars['D2'] + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'D2':
+            seconds = (sum(array) + (pars['D1']*2 + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'RD':
+            seconds = (sum(array) + (pars['D1']*2 + pars['D2']) * steps) * (pars['NS'] + pars['DS'])
+        else:
+            seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * steps * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for a variable parameter:
+        for index, pars[var_key] in enumerate(array):
+            print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\
+                             ' out of '+str(array.size)+': value = '+str(pars[var_key])
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield t2zz_experiment(pars, run) 
+            synchronize()
+    else:
+       # estimate the experiment time:
+        seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for accumulation:
+        for run in xrange(pars['NS']+pars['DS']):
+            yield t2zz_experiment(pars, run)
+
+# the pulse program: 
+
+def t2zz_experiment(pars, run):
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+        
+    pars['PROG'] = 't2zz_experiment'
+  
+    # eliminates T1 recovery and contaminating echo signals:
+    pars['PH1'] = [0, 180, 90, 270, 180, 0, 270, 90,       180, 0, 270, 90, 0, 180, 90, 270]
+    pars['PH3'] = [90, 90, 180, 180, 270, 270, 0, 0,       270, 270, 0, 0, 90, 90, 180, 180]	
+    pars['PH4'] = [0, 0, 90, 90, 180, 180, 270, 270,       180, 180, 270, 270, 0, 0, 90, 90]
+    pars['PH5'] = [0, 0, 90, 90, 180, 180, 270, 270,       0, 0, 90, 90, 180, 180, 270, 270]
+    pars['PH2'] = [0, 180, 90, 270, 180, 0, 270, 90,       0, 180, 90, 270, 180, 0, 270, 90]    # data routing
+
+    # read in variables:
+    P90 = pars['P90'] 
+    SF = pars['SF']
+    O1 = pars['O1']   
+    RD = pars['RD']
+    DEAD1 = pars['DEAD1'] 
+    D1 = pars['D1']
+    D2 = pars['D2']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH4 = pars['PH4'][run%len(pars['PH4'])]
+    PH5 = pars['PH5'][run%len(pars['PH5'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']
+
+    # sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+    
+    # run the pulse sequence: 
+    e.wait(RD)                                  	# delay between scans
+    e.set_frequency(SF+O1, phase=PH1)
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 1st pulse
+
+    e.wait(D1-P90/2-TXEnableDelay)     	              	# echo delay
+    e.set_phase(PH3)
+
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)		
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 2nd pulse (90- or 180-degree)
+    
+    e.wait(D1-P90/2-TXEnableDelay)     	              	# echo delay
+    e.set_phase(PH4)
+    
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)		
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)  # 3rd pulse (z-storage pulse)
+
+    e.wait(D2-P90/2-TXEnableDelay)    	              	# mixing time
+    e.set_phase(PH5)
+
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)		
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 4th pulse (readout pulse)
+
+    e.set_phase(PHA)                                    # set phase for receiver 
+    e.wait(DEAD1)                                       # wait for coil ringdown
+    e.record(SI, SW, sensitivity=ADCSensitivity)       	# acquire signal
+
+    # write experiment attributes:  
+    for key in pars.keys():
+        e.set_description(key, pars[key])		# acqusition parameters
+    e.set_description('run', run)               	# current scan
+    e.set_description('rec_phase', -PH2)		# current receiver phase
+
+    return e

Scripts/ZZ_Exchange_with_T2_Selection/op_t2zz_res.py

@@ -0,0 +1,173 @@
+# -*- coding: iso-8859-1 -*-
+
+from numpy import *
+from scipy.signal import *
+from os import path, rename 
+
+def result():
+   
+    measurement = MeasurementResult('Magnetization')
+
+    measurement_range = [0.0, 10e-6]
+    measurement_ranging = False
+
+    suffix = '' 	# output file name's suffix and...
+    counter = 1		# counter for arrayed experiments
+    var_key = ''	# variable parameter name
+        
+
+    # loop over the incoming results: 
+    for timesignal in results:
+        if not isinstance(timesignal,ADC_Result): 
+            continue
+
+        # read experiment parameters:
+        pars = timesignal.get_description_dictionary()
+
+        # ---------------- digital filter ------------------ 
+             
+        # get actual sampling rate of timesignal:
+        sampling_rate = timesignal.get_sampling_rate()
+
+        # get user-defined spectrum width:
+        spec_width = pars['SW']
+        
+        # specify cutoff frequency, in relative units:
+        cutoff = spec_width / sampling_rate
+                
+        if cutoff < 1: # otherwise no filter applied
+
+            # number of filter's coefficients:
+            numtaps = 29
+            
+            # use firwin to create a lowpass FIR filter:
+            fir_coeff = firwin(numtaps, cutoff)
+            
+            # downsize x according to user-defined spectral window:
+            skip = int(sampling_rate / spec_width)
+            timesignal.x = timesignal.x[::skip]
+            
+            for i in range(2):
+                # apply the filter to ith channel:
+                timesignal.y[i] = lfilter(fir_coeff, 1.0, timesignal.y[i])
+
+                # zeroize first N-1 "corrupted" samples:
+                timesignal.y[i][:numtaps-1] = 0.0
+
+                # circular left shift of y:
+                timesignal.y[i] = roll(timesignal.y[i], -(numtaps-1))    
+
+                # downsize y to user-defined number of samples (SI):
+                timesignal.y[i] = timesignal.y[i][::skip]
+
+            # update the sampling_rate attribute of the signal's:
+            timesignal.set_sampling_rate(spec_width)
+
+        # ----------------------------------------------------
+        
+        # sort timesignal data (data routing):
+        timesignal.phase(pars['rec_phase'])
+        
+        # provide timesignal to the display tab:
+        data['Current scan'] = timesignal
+        
+        # accumulate...
+        if not locals().get('accu'):
+            accu = Accumulation()
+
+        # skip dummy scans, if any:
+        if pars['run'] < 0: continue
+        
+        # add up:
+        accu += timesignal
+       # print 'job id: ', accu.get_job_id() 
+                
+        # provide accumulation to the display tab:
+        data['Accumulation'] = accu#(accu+0).magnitude()
+
+        # check how many scans are done: 
+        if accu.n == 1 and pars.has_key('COMMENT'):
+            print pars['COMMENT']
+       
+        if accu.n == pars['NS']: # accumulatioin is complete
+        
+            # correct accu's baselines:   
+            accu.baseline(last_part=0.2) 
+                    
+            # make a copy to process:
+            fid = accu + 0
+
+            # compute the initial phase of FID:
+            phi0 = arctan2(fid.y[1][0], fid.y[0][0]) * 180 / pi
+            if not 'ref' in locals(): ref = phi0   
+            print 'phi0 = ', phi0
+
+            # rotate FID to maximize y[0][0]:     
+            #fid.phase(-phi0) 
+
+            # do FFT:
+            fid.exp_window(line_broadening=10)
+            spectrum = fid.fft(samples=2*pars['SI'])
+
+            # try zero-order phase correction: 
+            spectrum.phase(-phi0)
+
+            # provide spectrum to the display tab:
+            data['Spectrum'] = spectrum
+
+            # check whether it is an arrayed experiment: 
+            var_key = pars.get('VAR_PAR')
+            if var_key:
+                # get variable parameter's value:
+                var_value = pars.get(var_key)
+                
+                # provide signal recorded with this var_value to the display tab:
+                data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu                
+
+                # measure signal intensity vs. var_value:
+                if measurement_ranging == True:
+                   [start, stop] = accu.get_sampling_rate() * array(measurement_range)
+                   measurement[var_value] = sum(accu.y[0][int(start):int(stop)]) 
+                   
+                else:
+                   measurement[var_value] = sum(accu.y[0][0:31])
+
+                # provide measurement to the display tab:
+                data[measurement.get_title()] = measurement
+
+                # update the file name suffix:
+                suffix = '_' + str(counter)
+                counter += 1
+
+            # save accu if required:
+            outfile = pars.get('OUTFILE')
+            if outfile: 
+                datadir = pars.get('DATADIR')
+
+                # write raw data in Simpson format:
+                filename = datadir+outfile+suffix+'.dat'
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.dat')
+                accu.write_to_simpson(filename)
+
+                # write raw data in Tecmag format:
+                # filename = datadir+outfile+'.tnt'
+                # accu.write_to_tecmag(filename, nrecords=20)
+
+               # write parameters in a text file:
+                filename = datadir+outfile+suffix+'.par' 
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.par')
+
+                fileobject = open(filename, 'w')
+                for key in sorted(pars.iterkeys()):
+                    if key=='run': continue
+                    if key=='rec_phase': continue
+                    fileobject.write('%s%s%s'%(key,'=', pars[key]))
+                    fileobject.write('\n')
+                fileobject.close()
+
+            # reset accumulation:
+            del accu
+            
+pass

Scripts/Zeeman_Order/op_zeeman_exp.py

@@ -0,0 +1,159 @@
+# -*- coding: iso-8859-1 -*-
+
+TXEnableDelay = 2e-6
+TXEnableValue = 0b0001	# TTL line blanking RF amplifier (bit 0)
+TXPulseValue  = 0b0010	# TTL line triggering RF pulses  (bit 1)
+ADCSensitivity = 2	# voltage span for ADC
+
+def experiment(): # Three-pulse STE sequence (Zeeman order)
+
+  # set up acquisition parameters:         
+    pars = {}
+    pars['P90'] = 1.7e-6	# 90-degree pulse length (s)
+    pars['SF'] = 338.7e6 	# spectrometer frequency (Hz)
+    pars['O1'] = -60e3		# offset from SF (Hz)
+    pars['SW'] = 1e6		# spectral window (Hz)
+    pars['SI'] = 1*1024		# number of acquisition points
+    pars['NS'] = 8		# number of scans
+    pars['DS'] = 0		# number of dummy scans
+    pars['RD'] = .5		# delay between scans (s)
+    pars['D1'] = 20e-6		# delay after first pulse, or 'short tau' (s)
+    pars['D2'] = 100e-6		# delay after second pulse, or 'long tau' (s)
+    pars['PHA'] = 30		# receiver phase (degree)
+    pars['DATADIR'] = '/home/fprak/Students/'	# data directory
+    pars['OUTFILE'] = None			# output file name
+
+  # specify a variable parameter (optional):
+    pars['VAR_PAR'] = 'D2'	# variable parameter name (a string)
+    start = 20e-6		# starting value
+    stop = 5e-0			# end value
+    steps = 24			# number of values
+    log_scale = True		# log scale flag 
+    stag_range = False		# staggered range flag
+
+  # check parameters for safety: 
+    if pars['PHA'] < 0:
+        pars['PHA'] = 360 + pars['PHA']
+
+    if pars['P90'] > 20e-6:
+        raise Exception("Pulse too long!!!")
+
+    # check whether a variable parameter is named: 
+    var_key = pars.get('VAR_PAR')
+    if var_key == 'P90' and (start > 20e-6 or stop > 20e-6):
+        raise Exception("Pulse too long!!!")
+
+    if pars['NS']%8 != 0:
+        pars['NS'] = int(round(pars['NS'] / 8) + 1) * 8 
+        print 'Number of scans changed to ',pars['NS'],' due to phase cycling'
+
+    # start the experiment:
+    if var_key:
+        # this is an arrayed experiment:
+        if log_scale:
+            array = log_range(start,stop,steps)
+        else:
+            array = lin_range(start,stop,steps)
+
+        if stag_range:
+            array = staggered_range(array, size = 2)
+
+        # estimate the experiment time:
+        if var_key == 'D1':
+            seconds = (sum(array)*2 + (pars['D2'] + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'D2':
+            seconds = (sum(array) + (pars['D1']*2 + pars['RD']) * steps) * (pars['NS'] + pars['DS'])
+        elif var_key == 'RD':
+            seconds = (sum(array) + (pars['D1']*2 + pars['D2']) * steps) * (pars['NS'] + pars['DS'])
+        else:
+            seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * steps * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for a variable parameter:
+        for index, pars[var_key] in enumerate(array):
+            print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\
+                             ' out of '+str(array.size)+': value = '+str(pars[var_key])
+            # loop for accumulation:
+            for run in xrange(pars['NS']+pars['DS']):
+                yield zeeman_experiment(pars, run) 
+            synchronize()
+    else:
+       # estimate the experiment time:
+        seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * (pars['NS']+ pars['DS'])
+        m, s = divmod(seconds, 60)
+        h, m = divmod(m, 60)
+        print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s)
+
+        # loop for accumulation:
+        for run in xrange(pars['NS']+pars['DS']):
+            yield zeeman_experiment(pars, run)
+
+# the pulse program: 
+
+def zeeman_experiment(pars, run):
+    e=Experiment()
+
+    dummy_scans = pars.get('DS')
+    if dummy_scans:
+        run -= dummy_scans
+        
+    pars['PROG'] = 'zeeman_experiment'
+
+   # 8-step phase cycle (1-21 modifided to deal with T1-recovery and extended for Re/Im imbalance)
+    pars['PH1'] = [0, 270,   0, 270,   180,  90, 180, 90] # 1st (90-degree) pulse
+    pars['PH3'] = [0,  90,   0,  90,     0,  90,   0, 90] # 2nd (90-degree) pulse
+    pars['PH4'] = [0,   0, 180, 180,   270, 270,  90, 90]  # 3rd (90-degree) pulse
+    pars['PH2'] = [0, 180, 180,   0,    90, 270, 270, 90]  # receiver
+
+   # read in variables:
+    P90 = pars['P90']
+    SF =  pars['SF']
+    O1 =  pars['O1']    
+    RD =  pars['RD']
+    D1 =  pars['D1']
+    D2 =  pars['D2']
+    PH1 = pars['PH1'][run%len(pars['PH1'])]
+    PH3 = pars['PH3'][run%len(pars['PH3'])]
+    PH4 = pars['PH4'][run%len(pars['PH4'])]
+    PH2 = pars['PH2'][run%len(pars['PH2'])]
+    PHA = pars['PHA']
+
+   # set sampling parameters:
+    SI = pars['SI']
+    SW = pars['SW']
+    while SW <= 10e6 and SI < 256*1024:
+        SI *= 2
+        SW *= 2
+
+   # run the pulse sequence:
+    e.wait(RD)                                  	# relaxation delay between scans
+    e.set_frequency(SF+O1, phase=PH1)
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse
+
+    e.wait(D1-P90/2-TXEnableDelay)			# 'short tau'
+
+    e.set_phase(PH3)
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse 
+
+    e.wait(D2-P90/2-TXEnableDelay)			# 'long tau'
+
+    e.set_phase(PH4)
+    e.ttl_pulse(TXEnableDelay, value=TXEnableValue)
+    e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue)	# 90-degree pulse
+
+    e.wait(TXEnableDelay)
+    e.set_phase(PHA)
+    e.wait(D1-P90/2-TXEnableDelay)			# 'short tau'
+    e.record(SI, SW, sensitivity=ADCSensitivity)	# acquisition
+
+   # write experiment parameters:
+    for key in pars.keys():
+        e.set_description(key, pars[key])	# acquisition parameters
+    e.set_description('run', run)		# current scan
+    e.set_description('rec_phase', -PH2)	# current receiver phase
+
+    return e

Scripts/Zeeman_Order/op_zeeman_res.py

@@ -0,0 +1,210 @@
+# -*- coding: iso-8859-1 -*-
+
+from numpy import *
+from scipy.signal import *
+from scipy.optimize import *
+from os import path, rename
+
+def result():
+
+    measurement = MeasurementResult('Magnetization')
+
+    measurement_range = [0.0, 10e-6]
+    measurement_ranging = False
+
+    suffix = '' 	# output file name's suffix and...
+    counter = 1		# counter for arrayed experiments
+    var_key = ''	# variable parameter name
+
+    # loop over the incoming results: 
+    for timesignal in results:
+        if not isinstance(timesignal,ADC_Result): 
+            continue
+
+        # read experiment parameters:
+        pars = timesignal.get_description_dictionary()
+
+        # ---------------- digital filter ------------------ 
+             
+        # get actual sampling rate of timesignal:
+        sampling_rate = timesignal.get_sampling_rate()
+
+        # get user-defined spectrum width:
+        spec_width = pars['SW']
+        
+        # specify cutoff frequency, in relative units:
+        cutoff = spec_width / sampling_rate
+                
+        if cutoff < 1: # otherwise no filter applied
+
+            # number of filter's coefficients:
+            numtaps = 29
+            
+            # use firwin to create a lowpass FIR filter:
+            fir_coeff = firwin(numtaps, cutoff)
+            
+            # downsize x according to user-defined spectral window:
+            skip = int(sampling_rate / spec_width)
+            timesignal.x = timesignal.x[::skip]
+            
+            for i in range(2):
+                # apply the filter to ith channel:
+                timesignal.y[i] = lfilter(fir_coeff, 1.0, timesignal.y[i])
+
+                # zeroize first N-1 "corrupted" samples:
+                timesignal.y[i][:numtaps-1] = 0.0
+
+                # circular left shift of y:
+                timesignal.y[i] = roll(timesignal.y[i], -(numtaps-1))    
+
+                # downsize y to user-defined number of samples (SI):
+                timesignal.y[i] = timesignal.y[i][::skip]
+
+            # update the sampling_rate attribute of the signal's:
+            timesignal.set_sampling_rate(spec_width)
+
+        # ----------------------------------------------------
+
+        # rotate timesignal according to current receiver's phase:
+        timesignal.phase(pars['rec_phase'])
+
+        # provide timesignal to the display tab:
+        data['Current scan'] = timesignal
+
+        # accumulate...
+        if not locals().get('accu'):
+            accu = Accumulation()
+
+        # skip dummy scans, if any:
+        if pars['run'] < 0: continue
+
+        # add up:
+        accu += timesignal
+
+        # provide accumulation to the display tab:
+        data['Accumulation'] = accu
+
+        # check how many scans are done:
+        if accu.n == pars['NS']: # accumulation is complete
+
+            # make a copy:
+            echo = accu + 0
+
+            # compute the initial phase of FID:
+            phi0 = arctan2(echo.y[1][0], echo.y[0][0]) * 180 / pi
+            if not 'ref' in locals(): ref = phi0   
+            print 'phi0 = ', phi0
+
+            # rotate FID to maximize y[0][0]:     
+            #echo.phase(-phi0) 
+
+            # do FFT:
+            echo.exp_window(line_broadening=10)
+            spectrum = echo.fft(samples=2*pars['SI'])
+
+            # try zero-order phase correction: 
+            spectrum.phase(-phi0)
+
+            # provide spectrum to the display tab:
+            data['Spectrum'] = spectrum
+
+            # check whether it is an arrayed experiment: 
+            var_key = pars.get('VAR_PAR')
+            if var_key:
+                # get variable parameter's value:
+                var_value = pars.get(var_key)
+                
+                # provide signal recorded with this var_value to the display tab:
+                data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu                
+
+                # measure signal intensity vs. var_value:
+                if measurement_ranging == True:
+                   [start, stop] = accu.get_sampling_rate() * array(measurement_range)
+                   measurement[var_value] = sum(accu.y[0][int(start):int(stop)]) 
+                   
+                else:
+                   measurement[var_value] = sum(accu.y[0][0:31])
+
+                # provide measurement to the display tab:
+                data[measurement.get_title()] = measurement
+
+                # update the file name suffix:
+                suffix = '_' + str(counter)
+                counter += 1
+
+            # save accu if required:
+            outfile = pars.get('OUTFILE')
+            if outfile: 
+                datadir = pars.get('DATADIR')
+
+                # write raw data in Simpson format:
+                filename = datadir+outfile+suffix+'.dat'
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.dat')
+                accu.write_to_simpson(filename)
+
+                # write raw data in Tecmag format:
+                # filename = datadir+outfile+'.tnt'
+                # accu.write_to_tecmag(filename, nrecords=20)
+
+               # write parameters in a text file:
+                filename = datadir+outfile+suffix+'.par' 
+                if path.exists(filename):
+                    rename(filename, datadir+'~'+outfile+suffix+'.par')
+
+                fileobject = open(filename, 'w')
+                for key in sorted(pars.iterkeys()):
+                    if key=='run': continue
+                    if key=='rec_phase': continue
+                    fileobject.write('%s%s%s'%(key,'=', pars[key]))
+                    fileobject.write('\n')
+                fileobject.close()
+
+            # reset accumulation:
+            del accu
+
+
+    if var_key == 'D2':
+        # KWW fit:
+        xdata = measurement.get_xdata()
+        ydata = measurement.get_ydata()
+        [amplitude, rate, beta] = fitfunc(xdata, ydata)
+        print '%s%02g' % ('Amplitude = ', amplitude)
+        print '%s%02g' % ('T2 [s] = ', 1./rate)
+        print '%s%02g' % ('Beta = ', beta)
+
+        # update display for the fit:
+        measurement.y = func([amplitude, rate, beta], xdata)
+        data[measurement.get_title()] = measurement
+
+# the fitting procedure:
+def fitfunc(xdata, ydata): 
+
+    # initialize variable parameters:
+    try: 
+        # solve Az = b:
+        A = array((ones(xdata.size/2), xdata[0:xdata.size/2]))    
+        b = log(abs(ydata[0:xdata.size/2]))     
+        z = linalg.lstsq(transpose(A), b)
+        amplitude = exp(z[0][0])
+        rate = -z[0][1]
+    except:
+        amplitude = abs(ydata[0])
+        rate = 1./(xdata[-1] - xdata[0])
+    beta = 1
+    p0 = [amplitude, rate, beta]
+
+    # run least-squares optimization:
+    plsq = leastsq(residuals, p0, args=(xdata, ydata))
+
+    return plsq[0] # best-fit parameters
+
+def residuals(p, xdata, ydata):
+    return ydata - func(p, xdata)
+
+# here is the function to fit:
+def func(p, xdata):
+    return p[0]*exp(-(p[1]*xdata)**p[2])
+
+
+pass