commit d994875c0fb7edf17f862902cf90b41a6b738980 Author: Markus Rosenstihl Date: Fri Jun 26 12:17:24 2015 +0000 initial check-in diff --git a/AU_Programs/2H/op_2h_exp.py b/AU_Programs/2H/op_2h_exp.py new file mode 100644 index 0000000..53d89ed --- /dev/null +++ b/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 \ No newline at end of file diff --git a/AU_Programs/2H/op_2h_res.py b/AU_Programs/2H/op_2h_res.py new file mode 100644 index 0000000..4e9acdb --- /dev/null +++ b/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 \ No newline at end of file diff --git a/AU_Programs/Diffusiometry/op_diff_exp.py b/AU_Programs/Diffusiometry/op_diff_exp.py new file mode 100644 index 0000000..511626a --- /dev/null +++ b/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 \ No newline at end of file diff --git a/AU_Programs/Diffusiometry/op_diff_res.py b/AU_Programs/Diffusiometry/op_diff_res.py new file mode 100644 index 0000000..48858e4 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/CPMG/op_cpmg_exp.py b/Scripts/CPMG/op_cpmg_exp.py new file mode 100644 index 0000000..a8f94cc --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/CPMG/op_cpmg_res.py b/Scripts/CPMG/op_cpmg_res.py new file mode 100644 index 0000000..bb14b6e --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/FID/op_fid_exp.py b/Scripts/FID/op_fid_exp.py new file mode 100644 index 0000000..32a3a5b --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/FID/op_fid_res.py b/Scripts/FID/op_fid_res.py new file mode 100644 index 0000000..0d903a1 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/FID_with_Background_Suppression/op_zgbs_exp.py b/Scripts/FID_with_Background_Suppression/op_zgbs_exp.py new file mode 100644 index 0000000..c999682 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/FID_with_Background_Suppression/op_zgbs_res.py b/Scripts/FID_with_Background_Suppression/op_zgbs_res.py new file mode 100644 index 0000000..0d903a1 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Hahn_Echo/op_hahn_exp.py b/Scripts/Hahn_Echo/op_hahn_exp.py new file mode 100644 index 0000000..91d0fff --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Hahn_Echo/op_hahn_res.py b/Scripts/Hahn_Echo/op_hahn_res.py new file mode 100644 index 0000000..64dd8b4 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Miscellaneous/op_gs_exp.py b/Scripts/Miscellaneous/op_gs_exp.py new file mode 100644 index 0000000..7f9ec60 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Miscellaneous/op_gs_res.py b/Scripts/Miscellaneous/op_gs_res.py new file mode 100644 index 0000000..f4f1553 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Saturation_Recovery/op_satrec_exp.py b/Scripts/Saturation_Recovery/op_satrec_exp.py new file mode 100644 index 0000000..02515fd --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Saturation_Recovery/op_satrec_res.py b/Scripts/Saturation_Recovery/op_satrec_res.py new file mode 100644 index 0000000..32a96b4 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Saturation_Recovery_with_Solid_Echo_Detection/op_satrec2_exp.py b/Scripts/Saturation_Recovery_with_Solid_Echo_Detection/op_satrec2_exp.py new file mode 100644 index 0000000..28c1ab4 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Saturation_Recovery_with_Solid_Echo_Detection/op_satrec2_res.py b/Scripts/Saturation_Recovery_with_Solid_Echo_Detection/op_satrec2_res.py new file mode 100644 index 0000000..32a96b4 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Solid_Echo/op_solidecho_exp.py b/Scripts/Solid_Echo/op_solidecho_exp.py new file mode 100644 index 0000000..4cd6f79 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Solid_Echo/op_solidecho_res.py b/Scripts/Solid_Echo/op_solidecho_res.py new file mode 100644 index 0000000..8aec750 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Spin_Alignment/op_spinal_exp.py b/Scripts/Spin_Alignment/op_spinal_exp.py new file mode 100644 index 0000000..4d7d1fa --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Spin_Alignment/op_spinal_res.py b/Scripts/Spin_Alignment/op_spinal_res.py new file mode 100644 index 0000000..5dc8a3c --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Spin_Alignment_Spin32/op_spinal32_exp.py b/Scripts/Spin_Alignment_Spin32/op_spinal32_exp.py new file mode 100644 index 0000000..6e517b5 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Spin_Alignment_Spin32/op_spinal32_res.py b/Scripts/Spin_Alignment_Spin32/op_spinal32_res.py new file mode 100644 index 0000000..5dc8a3c --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Steady_Gradient_Spin_Echo/op_sgse_exp.py b/Scripts/Steady_Gradient_Spin_Echo/op_sgse_exp.py new file mode 100644 index 0000000..c70e533 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Steady_Gradient_Spin_Echo/op_sgse_res.py b/Scripts/Steady_Gradient_Spin_Echo/op_sgse_res.py new file mode 100644 index 0000000..d5144d2 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Steady_Gradient_Spin_Echo_with_CPMG_Detection/op_sgse2_exp.py b/Scripts/Steady_Gradient_Spin_Echo_with_CPMG_Detection/op_sgse2_exp.py new file mode 100644 index 0000000..81102bf --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Steady_Gradient_Spin_Echo_with_CPMG_Detection/op_sgse2_res.py b/Scripts/Steady_Gradient_Spin_Echo_with_CPMG_Detection/op_sgse2_res.py new file mode 100644 index 0000000..f8b7d2c --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Stimulated_Echo/op_ste_exp.py b/Scripts/Stimulated_Echo/op_ste_exp.py new file mode 100644 index 0000000..29d7e21 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Stimulated_Echo/op_ste_res.py b/Scripts/Stimulated_Echo/op_ste_res.py new file mode 100644 index 0000000..1007a1a --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/T1Q/op_t1q_exp.py b/Scripts/T1Q/op_t1q_exp.py new file mode 100644 index 0000000..3200415 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/T1Q/op_t1q_res.py b/Scripts/T1Q/op_t1q_res.py new file mode 100644 index 0000000..206f025 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/ZZ_Exchange_with_T2_Selection/op_t2zz_exp.py b/Scripts/ZZ_Exchange_with_T2_Selection/op_t2zz_exp.py new file mode 100644 index 0000000..033f32a --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/ZZ_Exchange_with_T2_Selection/op_t2zz_res.py b/Scripts/ZZ_Exchange_with_T2_Selection/op_t2zz_res.py new file mode 100644 index 0000000..ae90fd2 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Zeeman_Order/op_zeeman_exp.py b/Scripts/Zeeman_Order/op_zeeman_exp.py new file mode 100644 index 0000000..91db8e1 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Zeeman_Order/op_zeeman_res.py b/Scripts/Zeeman_Order/op_zeeman_res.py new file mode 100644 index 0000000..3f5cb26 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/The DAMARIS Script Library.pdf b/The DAMARIS Script Library.pdf new file mode 100644 index 0000000..cc842fa Binary files /dev/null and b/The DAMARIS Script Library.pdf differ