From 8cdd3a73efc75f7132116085ab84de1320680fa2 Mon Sep 17 00:00:00 2001 From: Markus Rosenstihl Date: Tue, 27 Nov 2018 17:53:06 +0100 Subject: [PATCH 1/2] Added phase control testing for RF sources --- Scripts/Miscellaneous/PhaseTest/README.md | 7 ++++++ Scripts/Miscellaneous/PhaseTest/phase_exp.py | 25 ++++++++++++++++++++ Scripts/Miscellaneous/PhaseTest/phase_res.py | 23 ++++++++++++++++++ 3 files changed, 55 insertions(+) create mode 100644 Scripts/Miscellaneous/PhaseTest/README.md create mode 100755 Scripts/Miscellaneous/PhaseTest/phase_exp.py create mode 100755 Scripts/Miscellaneous/PhaseTest/phase_res.py diff --git a/Scripts/Miscellaneous/PhaseTest/README.md b/Scripts/Miscellaneous/PhaseTest/README.md new file mode 100644 index 0000000..1cab6ee --- /dev/null +++ b/Scripts/Miscellaneous/PhaseTest/README.md @@ -0,0 +1,7 @@ +Phase Test +========== +This script helps in testing the correct function of the phase programming +of RF sources. It measures a first and second, albeit phase shifted, +interval and calculatets the dot product of these to intervals. +The result should be a smooth full period of a cosine. + diff --git a/Scripts/Miscellaneous/PhaseTest/phase_exp.py b/Scripts/Miscellaneous/PhaseTest/phase_exp.py new file mode 100755 index 0000000..c5b5e92 --- /dev/null +++ b/Scripts/Miscellaneous/PhaseTest/phase_exp.py @@ -0,0 +1,25 @@ +import numpy as N +# connect PTS RF source to ADC card, set to generate a few MHz +# result should be a fulle period of sine/cosine +def experiment(): + #step = 0.36 # for PTS 500 + step = 0.225 # for PTS 310 + for i in list(xrange(1)): # number of periods + samples = 2048 + print samples + for phase in N.arange(0,360,step): + #phase=ph + e=Experiment() + #e.set_frequency(1e6,0) + e.wait(20e-6) + e.set_description("phase", phase) + e.set_description("i", i) + e.set_phase(0.0) + #e.loop_start(1) + e.wait(20e-6) + e.record(samples, 20e6,sensitivity=10) + e.set_phase(phase) + e.wait(20e-6) + #e.loop_end() + e.record(samples, 20e6,sensitivity=10) + yield e diff --git a/Scripts/Miscellaneous/PhaseTest/phase_res.py b/Scripts/Miscellaneous/PhaseTest/phase_res.py new file mode 100755 index 0000000..960cd75 --- /dev/null +++ b/Scripts/Miscellaneous/PhaseTest/phase_res.py @@ -0,0 +1,23 @@ +import numpy + +def result(): + o=MeasurementResult("overview") + for r in results: + if r is None: continue + #print r + data["single scan"]=r + phase = r.get_description("phase") + i = int(r.get_description("i")) + r1=numpy.array(r.get_result_by_index(0).y[0], dtype="Float64") + r2=numpy.array(r.get_result_by_index(1).y[0], dtype="Float64") + #print r1 + r1-=r1.mean() + r1/=r1.std() + r2-=r2.mean() + r2/=r2.std() + + c = numpy.dot(r1, r2) + print c + o[phase+360*i]=AccumulatedValue(c) + + data["overview"]=o \ No newline at end of file From cae15e0f1a4483e116c6c90a2950c9921541b209 Mon Sep 17 00:00:00 2001 From: Markus Rosenstihl Date: Tue, 27 Nov 2018 18:29:29 +0100 Subject: [PATCH 2/2] PFG STE added --- README.md | 26 +++ Scripts/PFG/Stimulated_Echo/pfg_ste_exp.py | 174 +++++++++++++++++ Scripts/PFG/Stimulated_Echo/pfg_ste_res.py | 209 +++++++++++++++++++++ 3 files changed, 409 insertions(+) create mode 100644 Scripts/PFG/Stimulated_Echo/pfg_ste_exp.py create mode 100644 Scripts/PFG/Stimulated_Echo/pfg_ste_res.py diff --git a/README.md b/README.md index f20ac55..ac66edc 100644 --- a/README.md +++ b/README.md @@ -7,4 +7,30 @@ You can get a clone of this repository via: Details of the pulse programs can be found [[https://chaos3.fkp.physik.tu-darmstadt.de/diffusion/DSL/browse/master/The%20DAMARIS%20Script%20Library.pdf | here]]. + +## For people working with this and thy want to revert the local changes: + +(This is from [[https://stackoverflow.com/questions/1146973/how-do-i-revert-all-local-changes-in-git-managed-project-to-previous-state|from StackOverflow]]) +If you want to revert changes made to your working copy, do this: + + git checkout . + +If you want to revert changes made to the index (i.e., that you have added), do this. Warning this will reset all of your unpushed commits to master!: + + git reset + +If you want to revert a change that you have committed, do this: + + git revert + +If you want to remove untracked files (e.g., new files, generated files): + + git clean -f + +Or untracked directories (e.g., new or automatically generated directories): + + git clean -fd + + + Current maintainer of this library is @markusro. diff --git a/Scripts/PFG/Stimulated_Echo/pfg_ste_exp.py b/Scripts/PFG/Stimulated_Echo/pfg_ste_exp.py new file mode 100644 index 0000000..1aa51d1 --- /dev/null +++ b/Scripts/PFG/Stimulated_Echo/pfg_ste_exp.py @@ -0,0 +1,174 @@ +# -*- coding: iso-8859-1 -*- + +TXEnableDelay = 2e-6 # test +TXEnableValue = 0b0001 # TTL line enabling RF amplifier (bit 0) +TXPulseValue = 0b0010 # TTL line triggering RF pulses (bit 1) +ADCSensitivity = 2 # voltage span for ADC +DAC_conv = 6.32e-5 # T/dac_value + +def experiment(): # stimulated echo experiment + + # set up acquisition parameters: + pars = {} + pars['P90'] = 3.2e-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'] = 8*1024 # number of acquisition points + pars['NS'] = 16*2 # number of scans + pars['DS'] = 0 # number of dummy scans + pars['RD'] = 3*7.7 # delay between scans (s) + pars['D1'] = 1e-3 # delay after first pulse (short tau) (s) + pars['D2'] = 50e-3 # delay after second pulse (long tau) (s) + pars['D4'] = 0e-6 # echo pre-acquisition delay (s) + pars["DAC1"] = 1000 # DAC value (PFG) + pars["D5"] = 0.9e-3 # PFG pulse length + pars['PHA'] = 30+180 # receiver phase (degree) + pars['DATADIR'] = '/home/markusro/STE' # data directory + pars['OUTFILE'] = "" # output file name + + # specify a variable parameter (optional): + pars['VAR_PAR'] = "DAC1" # variable parameter name (a string) + start = 0 # starting value + stop = int(5/DAC_conv) #1.5e5 # end value + steps = 21 # 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'] + D5 = pars['D5'] + DAC1 = pars['DAC1'] + 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.set_phase(PH3) + + e.set_pfg(dac_value=DAC1, length=D5, shape=("rec",100e-6)) + e.wait(D1-P90/2-TXEnableDelay - D5) # 'short tau' + + 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.set_pfg(dac_value=DAC1, length=D5, shape=("rec",100e-6)) + e.wait(D1-P90/2-TXEnableDelay+D4-D5) # '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/PFG/Stimulated_Echo/pfg_ste_res.py b/Scripts/PFG/Stimulated_Echo/pfg_ste_res.py new file mode 100644 index 0000000..1007a1a --- /dev/null +++ b/Scripts/PFG/Stimulated_Echo/pfg_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