Thorlabs OSA and thorlabs rotation stage

StageAndThorUpdated4.py

# -*- coding: utf-8 -*-
"""
StageAndThor.py
By Grace Kerber and Dan Hicktein (danhickstein@gmail.com)

This program connects a Thorlabs Rotation Stage and a Thorlabs OSA (ours OSA 205) 
# to take power scans while recording the spectra
"""
# This program defines zero as the position when stage is plugged into computer, will home in first part of code

import ctypes as c
import numpy as np
import os, time
import sys
import matplotlib.pyplot as plt
import datetime

t=time.time()

#Variable input Section
#Determines what angles to collect data at
ExtinguishAngle = 102 # Angle at which minimum power occurs
MaxPower        = 102. #Maximum Power through half wave plate
MinPower        = 1.3  # Minimum Power through half wave plate (power measured at ExtinguishAngle)
NumOfPoints     = 100  # number of data points (different power) to be collected 

avenum = 10  #Number of Thor OSA Spectra to average


#Accesses the DLLs that correspond to the rotation stage and Thor OSA
dllnameStage = os.path.join(os.path.dirname(__file__), 'Thorlabs.MotionControl.IntegratedStepperMotors.dll')

if not os.path.exists(dllnameStage):
    print "ERROR: DLL not found"

#p = c.windll.LoadLibrary(dllnameStage)  #Alternate between dll loading method
p = c.CDLL(dllnameStage) #Alternate between dll loading method

dllnameOSA = os.path.join('C:\\Program Files\\Thorlabs\\Thorlabs OSA\\FTSLib.dll')

print(dllnameOSA)
if not os.path.exists(dllnameOSA):
    print("ERROR: DLL not found")
#o = c.windll.LoadLibrary(dllnameOSA)  #Alternate between dll loading method
o = c.CDLL(dllnameOSA) #Altnate between dll loading method

#p.function for rotation stage
#o.function for Thor OSA

#Function calculates power at each angle from input power information
def angleFromPower(power, minPower=MinPower, maxPower=MaxPower, extinguishingAngle=ExtinguishAngle):
    return -(np.arcsin( ((power-minPower)/(maxPower-minPower))**0.5))*(180/np.pi)/2.  + extinguishingAngle

powerarray = np.linspace(MaxPower, MinPower, NumOfPoints)
anglearray = angleFromPower(powerarray)


#Function collects information from the Rotation Stage
def getHardwareInfo(SN):
    
    modelNo =         c.c_buffer(255)
    sizeOfModelNo =   c.c_ulong(255)
    hardwareType =    c.c_ushort()
    numChannels =     c.c_short()
    notes =           c.c_buffer(255)
    sizeOfNotes =     c.c_ulong(255)
    firmwareVersion = c.c_ulong()
    hardwareVersion = c.c_ushort()
    modState        = c.c_ushort()
    # p.PCC_GetHardwareInfo(SN)

    p.ISC_GetHardwareInfo(SN, 
                          c.pointer(modelNo), 
                          c.pointer(sizeOfModelNo), 
                          c.pointer(hardwareType),
                          c.pointer(numChannels),
                          c.pointer(notes),
                          c.pointer(sizeOfNotes),
                          c.pointer(firmwareVersion),
                          c.pointer(hardwareVersion),
                          c.pointer(modState) )
    

    return [x.value for x in (modelNo, sizeOfModelNo, hardwareType, numChannels, notes, sizeOfNotes, firmwareVersion, hardwareVersion, modState)]


#Function collectes motor parameter from Rotation Stage
def getMotorParamsExt(SN):
	# this doesn't work for some reason...
	stepsPerRev =  c.c_double()
	gearBoxRatio = c.c_double()
	pitch =        c.c_double()

	p.ISC_GetMotorParamsExt(SN, c.pointer(stepsPerRev), 
							       c.pointer(gearBoxRatio), 
							       c.pointer(pitch))

	return stepsPerRev.value, gearBoxRatio.value, pitch.value


#Function collects all Thorlabs rotation stages connected to the computer
def getDeviceList():
    p.TLI_BuildDeviceList()
    receiveBuffer = c.c_buffer(200)
    sizeOfBuffer  = c.c_ulong(255)
    p.TLI_GetDeviceListExt(c.pointer(receiveBuffer), c.pointer(sizeOfBuffer))
    return [x for x in (receiveBuffer.value).split(',')[:-1]]


#Function moves Rotation Stage to specific angle, in terms of device units (needs to be converted from angle)
def MoveToPosition(SN, deviceUnits, timeout=20, queryDelay=0.01, tolerance=1):
    """
    Moves the rotation stage to a certain position (given by device units).
    This call blocks future action until the move is complete.
    The timeout is in seconds
    
    SN is a c_buffer of the serial number string
    deviceUnits shold be a int.
    tolerance is when the blocking should end (device units)
    """
    
    GetStatus(SN)
    p.ISC_MoveToPosition(SN, c.c_int(int(deviceUnits)))

    t = time.time()

    while time.time()<(t+timeout):
        GetStatus(SN)
        p.ISC_RequestStatus(SN) # order the stage to find out its location
        currentPosition = p.ISC_GetPosition(SN)
        error = currentPosition - deviceUnits
        if np.abs(error)<tolerance: 
            return
        else:
            time.sleep(queryDelay)
    raise ValueError('Oh no!!! We never got there!! Maybe you should make the timeout longer than %.3f seconds dude.'%timeout)


try:
    serialNumber = getDeviceList()[0]
except:
    raise ValueError('Couldn\'t get the list of serial numbers! Is your stage plugged in? Or is Thorlabs Kinesis open?')

#Function requests status of Rotation Stage, whether moving/not moving/homing/etc
def GetStatus(SN):
    p.ISC_RequestStatus(SN)
    #bits = p.ISC_GetStatusBits(SN)
    #print bin(bits)

#Function collects data from OSA with the x-axis unit of wavenumbers and the y-axis unit of mW
def GrabSpectrum():
    global grabbing, x, y, xcfloat, ycfloat, err1, err2
    grabbing = True
    #t3 = time.time()
    
    def retreiveSpectrum():
         global grabbing, x, y, xcfloat, ycfloat, err1, err2
         xcfloat = (c.c_float*16774)()
         ycfloat = (c.c_float*16774)()
         p = (c.c_float*16774)()
         f = c.c_int()
         l = c.c_uint()
         
         err1 = 10
         err2 = 10
 
         t2 = time.time()
         o.FTS_GetLastSpectrumArrays(sn, c.pointer(ycfloat), c.pointer(xcfloat), c.pointer(p), c.pointer(f), c.pointer(l))  #The spectrum has the x-axis unit of wavenumbers and the y-axis unit of mW. 
         err1 = time.time()-t2
  
         y = np.array(list(ycfloat))
         x = np.array(list(xcfloat))        
  
         
    #     o.FTS_GetLastSpectrumArrays(sn, x, y, p, fmt, length)
         
         
         grabbing = False
         
             
    def callback(p1,p2,p3):
        if p2 == 1:
            #print('Good news from the OSA: Interferogram now available!')
            a=0
        elif p2 == 2:
            #print('Good news from the OSA: Spectrum now available!')
            a=0        
            retreiveSpectrum()
        elif p2 == 3:
            #print('Good news from the OSA: One of the average spectra has been collected, but the averaged spectrum is not yet available.')
            a=0
        else:
            #print('Very strange news from the spectrometer. p1,p2,p3:',p1,p2,p3)
            a=0
        return a
    

    callback_format = c.CFUNCTYPE(c.c_voidp ,c.c_ushort, c.c_uint, c.c_uint)
    CALL_BACK = callback_format(callback)
    #t3 = time.time()
    #print('Acquiring spectrum...')    
    o.FTS_AcquireSingleSpectrum(sn, CALL_BACK)
   
    while grabbing==True:
        time.sleep(0.005)
    #err2 = time.time() -t3    
    return x,y,xcfloat,ycfloat, err1, err2
    #o.FTS_StopContinuousAcquisition(sn)


#Makes new data file folder if needed and run folder for each run that day    
#---Base Directory
today = datetime.datetime.now().strftime("%Y-%m-%d")
cwd = os.getcwd()
base_dir = os.path.join(cwd, today)
if not(os.path.isdir(base_dir)):
    os.mkdir(base_dir)

run_counter = 1 
run_folder = 'run %04i'%(run_counter)
while os.path.isdir(os.path.join(base_dir, run_folder)):
    run_counter = run_counter+1
    run_folder = 'run %04i'%(run_counter)
new_base_dir = os.path.join(base_dir,run_folder)
os.mkdir(new_base_dir)  

print 'Will save Data Log to run folder %i' %(run_counter);


#indexing of rotation stage
SN = c.c_buffer(serialNumber) #index of Rotation Stage to be used
try:
    p.ISC_Close(SN)
    print 'Previous Stage connection closed.'
except:
    pass
openval = p.ISC_Open(SN)

#indexing of Thor OSA
sn = c.c_int(0) # index of the OSA to be used. 
try:
    o.FTS_CloseSpectrometer(sn)
    print 'Previous OSA connection closed.'
except:
    pass
print('Initializing spectrometer(s)...')

#SN for rotation stage
#sn for Thor OSA



numSpectrometers = o.FTS_InitializeSpectrometers()
print('There is/are (%i) Thorlabs OSAs connected'%numSpectrometers)


o.FTS_OpenSpectrometer(sn)
o.FTS_SetAcquisitionOption_AverageSpectrum(sn, c.c_int(1))   

hardwareinfoval = getHardwareInfo(SN)
startpolval = p.ISC_StartPolling(SN,c.c_int(200))
loadsettingsval = p.ISC_LoadSettings(SN)


#Writes array of half-wave plate angles and calculated input powers to file directory 
listpowerarray = []
listanglearray = []
for count, (powerarray) in enumerate(powerarray):
    angle = -(np.arcsin( ((powerarray-MinPower)/(MaxPower-MinPower))**0.5))*(180/np.pi)/2.  + ExtinguishAngle
    listanglearray.append(str('%08.4f'%(angle)))
    listpowerarray.append((str('%08.4f'%(powerarray))))
  
  
time_now = datetime.datetime.now().strftime("%Y-%m-%d %X")
col_power = ["Instrument:"] + ["Time Stamp:"] + ["Min Pow Angle"] + ["Max Power"] + ["Min Power"] + ["Total Steps"] + ["", "Power (mW)"] + listpowerarray
col_angle = ["Thor OSA"] + [time_now] + [str(ExtinguishAngle)] + [str(MaxPower)] + [str(MinPower)] + [str(NumOfPoints)] + ["", "Angle (degree)"] + listanglearray

col_comb_LOG = zip(col_power, col_angle)
data_list_LOG = []
for data_row_LOG in col_comb_LOG:
    data_list_LOG.append('\t'.join(data_row_LOG))
data_string_LOG = "\n".join(data_list_LOG)

#Saving Data LOG  
#---Base File Name
file_name_LOG = 'DATALOG_THOR_osa-'+today+'_0001.txt'
#---Check for Pre-existing Files
file_counter_LOG = 1
while os.path.isfile(os.path.join(new_base_dir,file_name_LOG)):
    file_counter_LOG = file_counter_LOG+1
    file_name_LOG = ('DATALOG_osa-'+today+'_%.4i.txt')%file_counter_LOG
#---Save Data File
data_file_LOG = open(os.path.join(new_base_dir,file_name_LOG), 'w')
data_file_LOG.write(data_string_LOG)

data_file_LOG.close()


print 'Homing the rotation stage...'; sys.stdout.flush()
p.ISC_Home(SN)
while (p.ISC_GetStatusBits(SN))!=(-2147482624):  #While Motor is moving.  Stopped is -2147482624 (specific number to stage)
    #print 'in the process of homing'
    time.sleep(.05)
    
print 'Rotation Stage Has Been Homed'



stepsPerRev, gearBoxRatio, pitch = getMotorParamsExt(SN)
#print stepsPerRev, gearBoxRatio, pitch
microstepsPerFullstep = 2048 # from https://www.thorlabs.com/newgrouppage9.cfm?objectgroup_id=8750
conversion = stepsPerRev * microstepsPerFullstep * gearBoxRatio / pitch # convert to degrees
#print conversion


xaveall = []
yaveall = []
degree = anglearray
totalsteps = len(anglearray)

#o.FTS_SetAcquisitionOptions_Default(sn)
#o.FTS_SetAcquisitionOption_AverageSpectrum(sn, c.c_int(1)) 

#For loop goes through each angle data will be collected at - one run through for loop per angle
for count, (degree) in enumerate(degree[::]):
    DegreePosition = degree #value in degrees

    deviceUnits = abs(int(DegreePosition*conversion))#-deviceUnitsZero) #this involved rounding!
    print 'Moving to %5.3f degrees (%i Device Units)...'%(DegreePosition, deviceUnits); sys.stdout.flush()
    MoveToPosition(SN, deviceUnits)
    new_position = p.ISC_GetPosition(SN)
    new_degrees = new_position/conversion
    print 'Reported %5.3f degrees (%i Device Units).'%(new_degrees, new_position); sys.stdout.flush()

    
    count1=0
    xall=[]
    yall=[]
    print 'taking {} spectrum(s) at {} degrees to be averaged together.'.format(avenum,new_degrees); sys.stdout.flush()
    t1=time.time()
    #While loop goes through taking spectra and averaging it for the number of times to average that was set in first few lines of program
    while count1<avenum:
        
        x,y,xcfloat,ycfloat,err1,err2 = GrabSpectrum()
        
       # x,y,xcfloat,ycfloat = GrabSpectrum()  #Gets spectrum from previously defined function
        count1=count1+1
        y = y[:-1] # get rid of the zero mW (absolute power) component
        x = x[:-1] # get rid of the zero cm-1 (wavenumber) component        
        grid_size =  x[1]-x[0]
        
        
        #Converts data from cm-1 & mW to nm & dBm/nm
        y_mW_per_wavenumber = y/grid_size
        nmx=(1E7)/x  #cm-1 to nm
    
        y_mW_per_nm = y_mW_per_wavenumber /( 1e7/(x**2))
        #ynew=y/((5E-8)*nmx*nmx)  #absolute power to relative power
        dby=10*np.log10(y_mW_per_nm) #log scale of relative power
        xall.append(nmx)
        yall.append(dby)
        xa=x
        ya=y
        
    #averages all data collected from while loop 
    print 'the time taken to collect and average {} spectra was {} seconds'.format(avenum, (time.time()-t1))

    yave = np.mean(yall,axis=0,dtype=np.float64)
    xave = np.mean(xall,axis=0,dtype=np.float64)
    xaveall.append(xave)
    yaveall.append(yave)
    im = plt.plot(xave, yave)  #plots each itteration of for loop, shown all on top of eachother at end of program
    #plt.show(im)  
   
    #Saves each itteration for loop data to new text file (one file per half-waveplate angle)    
    #---Time Stamp
    time_now = datetime.datetime.now().strftime("%Y-%m-%d %X")  
    OSA_Thor = o.FTS_GetInstrumentNum()  #Gets Thor OSA Instrumentation
    file_name_OSA = 'thor-osa-data_'+today+'_0001.txt'
    #---Check for Pre-existing Files
    file_counter_OSA = 1
    while os.path.isfile(os.path.join(new_base_dir,file_name_OSA)):
        file_counter_OSA = file_counter_OSA+1
        file_name_OSA = ('thor-osa-data_'+today+'_%.4i.txt')%file_counter_OSA
    #---Save Data File
    with open(os.path.join(new_base_dir,file_name_OSA), 'w') as data_file_OSA:
        data_file_OSA.write('Thor OSA instrument: \t {} \n'.format(OSA_Thor))
        data_file_OSA.write('Wavelength (nm)\t Intensity (dBm/nm)\n')
        for x,y in zip(xave, yave):
            data_file_OSA.write('%.4f\t%.4f\n'%(x,y))    
    print 'Saved Trace to {} out of {}'.format(file_counter_OSA, totalsteps); sys.stdout.flush()



#returns rotation stage back to zero
#print 'Moving Back to Position Zero'
#MoveToPosition(SN,0)
print('Moving Back to Max Power')
MoveToPosition(SN, abs(int((ExtinguishAngle+45)*conversion)))

plt.show()

#Close all things at end of the program
o.FTS_ClearLastSpectrum(sn)
closeval = p.ISC_Close(SN)
o.FTS_CloseSpectrometer(sn)

print 'the total time taken was {} minutes'.format((time.time()-t)/60)