In [ ]:
## Benchmark testing the code for finding the optimal green value for evaporation

In [2]:
import optevap

In [3]:
%%timeit 

optevap.optimal( s0=7., wL=45., wC=28.)


1 loops, best of 3: 140 ms per loop

Code that optimizes the trap given s0, wL, wG

When optimizing the trap parameters only the ratio wL/wC is relevant so we do not need to run nested wL, wC loops, just fix wL at some value and vary wG.


In [46]:
import optevap
s=7.0
aS=650.



def ext( wIR, alpha ):
    if alpha <= 0.85:
        e = 1.10*wIR
    elif alpha <=0.90:
        e = 1.20*wIR
    elif alpha < 0.95:
        e = 1.30*wIR
    elif alpha <= 1.:
        e = 1.40*wIR
    elif alpha <= 1.05:
        e = 1.45*wIR
    elif alpha <= 1.10:
        e = 1.50*wIR
    elif alpha <= 1.13:
        e = 1.50*wIR
    elif alpha <= 1.16:
        e = 1.50*wIR
    elif alpha <= 1.18:
        e = 1.50*wIR
    elif alpha <= 1.20:
        e = 1.50*wIR 
    elif alpha <= 1.30:
        e = 1.45*wIR
    else:
        e = 1.30*wIR
        
    estr =  "wIR=%.2f, alpha=%.3f --> extents=%.2f"%(wIR,alpha,e)
    return e, estr

savefile = 'alpha_opt/001.dat'
open(savefile, 'w').close()

wIR = 47. 
wGRs = np.arange( 28., 72., 0.5) 
for wGR in wGRs:
    
    if np.abs( wIR/wGR - 1.17 ) < 0.0: 
        continue
    else:
        print "  wGR=%d"%wGR,
    #print "wGR=%d.." % (wGR),
    T=0.2; extents, extstr = ext( wIR, wIR/wGR) 
    try:
        res =  np.array( [ [s, wIR, wGR]+ \
                optevap.optimal( s0=s, wL=wIR, wC=wGR, a_s=aS, T_Er=T, \
                                 extents=extents) ] )

        f_handle = file(savefile, 'a')
        np.savetxt(f_handle, res, fmt='%.4g')
        f_handle.close()

    except Exception as e:
        print " wIR=%d,wGR=%d->ERR"%(wIR,wGR)
        print extstr
        print e 
        raise


  wGR=28   wGR=28   wGR=29   wGR=29   wGR=30   wGR=30   wGR=31   wGR=31   wGR=32   wGR=32   wGR=33   wGR=33   wGR=34   wGR=34   wGR=35   wGR=35   wGR=36   wGR=36   wGR=37   wGR=37   wGR=38   wGR=38   wGR=39   wGR=39   wGR=40   wGR=40   wGR=41   wGR=41   wGR=42   wGR=42   wGR=43   wGR=43   wGR=44   wGR=44   wGR=45   wGR=45   wGR=46   wGR=46   wGR=47   wGR=47   wGR=48   wGR=48   wGR=49   wGR=49   wGR=50   wGR=50   wGR=51   wGR=51   wGR=52   wGR=52   wGR=53   wGR=53   wGR=54   wGR=54   wGR=55   wGR=55   wGR=56   wGR=56   wGR=57   wGR=57   wGR=58   wGR=58   wGR=59   wGR=59   wGR=60   wGR=60   wGR=61   wGR=61   wGR=62   wGR=62   wGR=63   wGR=63   wGR=64   wGR=64   wGR=65   wGR=65   wGR=66   wGR=66   wGR=67   wGR=67   wGR=68   wGR=68   wGR=69   wGR=69   wGR=70   wGR=70   wGR=71   wGR=71

We create a pandas data frame with the optimized data


In [47]:
import pandas as pd 

dat = np.loadtxt('alpha_opt/001.dat')
df = pd.DataFrame( dat, columns=['v0', 'wIR', 'wGR', 'g0', 'EtaF', 'N', 'SN', 'HWHM', 'HWHMwL', 'DeltaF'])
df['alpha'] = df['wIR'] / df['wGR'] 

df


Out[47]:
v0 wIR wGR g0 EtaF N SN HWHM HWHMwL DeltaF alpha
0 7 47 28.0 1.964 7.641 216900 1.624 19.73 0.4197 4.8340 1.678571
1 7 47 28.5 2.034 7.453 226400 1.626 20.11 0.4279 4.6970 1.649123
2 7 47 29.0 2.106 7.261 236600 1.628 20.49 0.4361 4.5570 1.620690
3 7 47 29.5 2.179 7.065 247600 1.629 20.49 0.4361 4.4140 1.593220
4 7 47 30.0 2.254 6.865 259300 1.631 20.88 0.4442 4.2690 1.566667
5 7 47 30.5 2.330 6.661 271700 1.633 21.26 0.4524 4.1200 1.540984
6 7 47 31.0 2.407 6.454 285000 1.635 21.64 0.4605 3.9690 1.516129
7 7 47 31.5 2.485 6.243 299200 1.638 22.03 0.4687 3.8160 1.492063
8 7 47 32.0 2.564 6.028 314500 1.640 22.41 0.4768 3.6600 1.468750
9 7 47 32.5 2.658 5.770 340400 1.635 22.79 0.4850 3.4740 1.446154
10 7 47 33.0 2.740 5.549 358900 1.638 23.18 0.4931 3.3120 1.424242
11 7 47 33.5 2.824 5.324 379000 1.641 23.56 0.5013 3.1490 1.402985
12 7 47 34.0 2.908 5.095 400900 1.645 24.33 0.5176 2.9820 1.382353
13 7 47 34.5 2.994 4.862 425100 1.649 24.71 0.5257 2.8130 1.362319
14 7 47 35.0 3.081 4.627 451800 1.654 25.09 0.5339 2.6410 1.342857
15 7 47 35.5 3.169 4.387 481600 1.659 25.47 0.5420 2.4670 1.323944
16 7 47 36.0 3.259 4.144 514900 1.665 26.24 0.5583 2.2900 1.305556
17 7 47 36.5 3.350 3.898 552600 1.671 26.70 0.5682 2.1100 1.287671
18 7 47 37.0 3.442 3.647 595900 1.678 27.56 0.5864 1.9280 1.270270
19 7 47 37.5 3.535 3.392 646300 1.687 27.99 0.5955 1.7420 1.253333
20 7 47 38.0 3.630 3.134 705500 1.696 28.84 0.6136 1.5540 1.236842
21 7 47 38.5 3.726 2.872 776000 1.708 29.70 0.6318 1.3640 1.220779
22 7 47 39.0 3.824 2.607 862000 1.723 30.55 0.6500 1.1710 1.205128
23 7 47 39.5 3.907 2.382 925800 1.752 31.16 0.6630 1.0060 1.189873
24 7 47 40.0 4.007 2.111 1054000 1.775 32.49 0.6912 0.8085 1.175000
25 7 47 40.5 4.107 1.836 1224000 1.806 33.81 0.7194 0.6084 1.160494
26 7 47 41.0 4.210 1.557 1464000 1.845 36.02 0.7665 0.4058 1.146341
27 7 47 41.5 4.268 1.401 1510000 1.929 35.58 0.7571 0.2915 1.132530
28 7 47 42.0 4.270 1.397 1236000 2.061 32.93 0.7006 0.2889 1.119048
29 7 47 42.5 4.268 1.405 982500 2.187 29.84 0.6348 0.2951 1.105882
30 7 47 43.0 4.274 1.393 801600 2.296 27.18 0.5784 0.2859 1.093023
31 7 47 43.5 4.273 1.398 637700 2.380 24.97 0.5313 0.2898 1.080460
32 7 47 44.0 4.277 1.390 519800 2.441 23.65 0.5031 0.2840 1.068182
33 7 47 44.5 4.279 1.389 424500 2.473 21.88 0.4655 0.2829 1.056180
34 7 47 45.0 4.271 1.412 341600 2.454 20.72 0.4409 0.3002 1.044444
35 7 47 45.5 4.280 1.390 293300 2.456 19.44 0.4136 0.2836 1.032967
36 7 47 46.0 4.266 1.432 241500 2.406 18.59 0.3955 0.3143 1.021739
37 7 47 46.5 4.264 1.439 208500 2.372 17.73 0.3773 0.3197 1.010753
38 7 47 47.0 4.304 1.414 196500 2.383 17.53 0.3730 0.3010 1.000000
39 7 47 47.5 4.304 1.497 173400 2.344 16.71 0.3555 0.3616 0.989474
40 7 47 48.0 4.304 1.588 154800 2.307 16.30 0.3467 0.4279 0.979167
41 7 47 48.5 4.304 1.678 139700 2.274 15.88 0.3379 0.4935 0.969072
42 7 47 49.0 4.304 1.767 127300 2.244 15.47 0.3292 0.5585 0.959184
43 7 47 49.5 4.304 1.856 116600 2.207 15.13 0.3219 0.6227 0.949495
44 7 47 50.0 4.304 1.943 107900 2.187 14.75 0.3138 0.6861 0.940000
45 7 47 50.5 4.304 2.029 100300 2.171 14.37 0.3056 0.7489 0.930693
46 7 47 51.0 4.304 2.114 93830 2.156 13.98 0.2975 0.8109 0.921569
47 7 47 51.5 4.304 2.198 88160 2.144 13.98 0.2975 0.8721 0.912621
48 7 47 52.0 4.304 2.281 83180 2.134 13.60 0.2893 0.9326 0.903846
49 7 47 52.5 4.304 2.364 78770 2.123 13.26 0.2821 0.9924 0.895238
50 7 47 53.0 4.304 2.445 74850 2.116 13.26 0.2821 1.0510 0.886792
51 7 47 53.5 4.304 2.525 71340 2.109 12.91 0.2746 1.1100 0.878505
52 7 47 54.0 4.304 2.604 68190 2.104 12.91 0.2746 1.1670 0.870370
53 7 47 54.5 4.304 2.681 65340 2.099 12.55 0.2671 1.2240 0.862385
54 7 47 55.0 4.304 2.758 62750 2.094 12.55 0.2671 1.2800 0.854545
55 7 47 55.5 4.304 2.834 60380 2.090 12.16 0.2586 1.3350 0.846847
56 7 47 56.0 4.304 2.909 58220 2.086 12.16 0.2586 1.3890 0.839286
57 7 47 56.5 4.304 2.983 56240 2.082 12.16 0.2586 1.4430 0.831858
58 7 47 57.0 4.304 3.056 54420 2.079 11.83 0.2517 1.4960 0.824561
59 7 47 57.5 4.304 3.128 52730 2.076 11.83 0.2517 1.5480 0.817391
... ... ... ... ... ... ... ... ... ... ...

88 rows × 11 columns

Then we plot the optimization results


In [55]:
from matplotlib import rc
rc('font',**{'family':'serif'})
rc('text', usetex=True)


color = ['blue', 'green', 'red', 'black']


dat = np.loadtxt('alpha_opt/001.dat')
df = pd.DataFrame( dat, columns=['v0', 'wIR', 'wGR', 'g0', 'EtaF', 'N', 'SN', 'HWHM', 'HWHMwL', 'DeltaF'])
    
T = 0.2 
    
fig = plt.figure(figsize=(8.5,2.4))
gs  = matplotlib.gridspec.GridSpec(1,4, wspace=0.36, hspace=0.35,\
        left=0.05, right=0.99, bottom=0.18, top=0.93)

axL = [ fig.add_subplot( gs[i] ) for i in range(4) ] 

axL[0].plot( df['wIR']/df['wGR'], df['g0'], 'o-', color='black', ms=3)
axL[1].plot( df['wIR']/df['wGR'], df['EtaF'], 'o-', color='black', ms=3)
axL[2].plot( df['wIR']/df['wGR'], df['DeltaF'], 'o-', color='black', ms=3)
axL[3].plot( df['wIR']/df['wGR'], df['SN'], 'o-', color='black', ms=3)


YLAB = ['$g_{0}$', '$\eta_{F}$', '$\Delta_{F}$', '$S/N$']        
for i,ax in enumerate(axL):
    ax.set_xlabel(r'$\alpha_{\mathrm{w}}$')
    ax.set_xlim(0.85,1.22)
    ax.set_ylabel( YLAB[i] ) 
    ax.grid(alpha=0.3)
    ax.xaxis.set_major_locator( matplotlib.ticker.MultipleLocator(0.1))
    ax.xaxis.set_minor_locator( matplotlib.ticker.MultipleLocator(0.05))
    ax.tick_params(axis='both', which='major', length=2.5)
    ax.tick_params(axis='both', which='minor', length=1.5)
    
    ax.axvline(1.136, color='0.2', alpha=0.5, lw=2.)
    ax.text( 0.78, 1.01,  r'$\alpha_{\mathrm{evap}}$', \
            ha='center', va='bottom', transform=ax.transAxes)
    
    
    
axL[0].set_ylim(3.2,4.6)
axL[1].set_ylim(0.8, 3.2)
axL[2].set_ylim(0.,1.8)

#axL[4].set_ylim(0.,10.)
#axL[4].legend( bbox_to_anchor = (0.,1.1), loc='upper left', \
#        numpoints=1, labelspacing=0.2, prop={'size':8},\
#        handlelength=1.1, handletextpad=0.5 )




fig.savefig('alpha_opt/001.png', dpi=300)


Looking at atom number at optimized $g_0$

We now consider the length scale of the system by varying $w_{L}$. The optimization does not need to be run again, becaue the results are lengthscale independent and we have them available from above.


In [36]:
import optevap
s=7.0
aS=650.

import pandas as pd
dat = np.loadtxt('alpha_opt/001.dat')
df = pd.DataFrame( dat, columns=['v0', 'wIR', 'wGR', 'g0', 'EtaF', 'N', 'SN', 'HWHM', 'HWHMwL', 'DeltaF'])
df['alpha'] = df['wIR'] / df['wGR']

from scipy.interpolate import interp1d
g0f = interp1d( df['alpha'].values[::-1], df['g0'].values[::-1] )


def ext( wIR, alpha ):
    if alpha <= 0.85:
        e = 1.10*wIR
    elif alpha <=0.90:
        e = 1.20*wIR
    elif alpha < 0.95:
        e = 1.30*wIR
    elif alpha <= 1.:
        e = 1.40*wIR
    elif alpha <= 1.05:
        e = 1.45*wIR
    elif alpha <= 1.10:
        e = 1.50*wIR
    elif alpha <= 1.13:
        e = 1.50*wIR
    elif alpha <= 1.16:
        e = 1.50*wIR
    elif alpha <= 1.18:
        e = 1.50*wIR
    elif alpha <= 1.20:
        e = 1.50*wIR 
    elif alpha <= 1.30:
        e = 1.45*wIR
    else:
        e = 1.30*wIR
        
    estr =  "wIR=%.2f, alpha=%.3f --> extents=%.2f"%(wIR,alpha,e)
    return e, estr

savefile = 'alpha_opt/002.dat'
open(savefile, 'w').close()
    
wIRs = [45., 50., 55., 60., 65., 70.]
for wIR in wIRs:
    print "  wIR=%d"%wIR,
    wGRs = np.arange( 28., 72., 1.) 
    for wGR in wGRs:
        
        alpha = wIR/wGR 
        if alpha < 0.8 or alpha > 1.3:
            continue 
    
        T=0.2; extents, extstr = ext( wIR, wIR/wGR) 
        try:
            g0 = g0f(alpha) -0.02
            res =  np.array( [ [s, wIR, wGR]+ \
                    optevap.get_trap_results( s0=s, wL=wIR, wC=wGR, a_s=aS, T_Er=T, \
                                     extents=extents,  g0=g0) ] )

            f_handle = file(savefile, 'a')
            np.savetxt(f_handle, res, fmt='%.4g')
            f_handle.close()

        except Exception as e:
        
            print "\n\nwIR=%d,wGR=%d->ERR"%(wIR,wGR)
            print alpha
            print g0f(alpha) - 0.02
            print optevap.optimal( s0=s, wL=wIR, wC=wGR, a_s=aS, T_Er=T, \
                                     extents=extents)
            print extstr
            print e 
            raise


  wIR=45   wIR=50   wIR=55   wIR=60   wIR=65   wIR=70

In [65]:
from matplotlib import rc
rc('font',**{'family':'serif'})
rc('text', usetex=True)


color = ['blue', 'green', 'red', 'black']


dat = np.loadtxt('alpha_opt/002.dat')
df = pd.DataFrame( dat, columns=['v0', 'wIR', 'wGR', 'g0', 'EtaF', 'N', 'SN', 'HWHM', 'HWHMwL', 'DeltaF'])

    
T = 0.2 
    
fig = plt.figure(figsize=(5.3,2.3))
gs  = matplotlib.gridspec.GridSpec(1,2, wspace=0.36, hspace=0.35,\
        left=0.08, right=0.99, bottom=0.16, top=0.93)

axL = [ fig.add_subplot( gs[i] ) for i in range(2) ] 

for wIR in np.unique(df['wIR']):
    subdf = df[ df['wIR'] == wIR ]
    axL[0].plot( subdf['wIR']/subdf['wGR'], subdf['N']/1e5, 'o-', ms=3, label=r'${:d}$'.format(int(wIR)))
    axL[1].plot( subdf['wIR']/subdf['wGR'], subdf['HWHMwL'], 'o-', ms=3)



YLAB = ['$N/10^{5}$', '$\mathrm{HWHM}/w_{L}$']        
for i,ax in enumerate(axL):
    ax.set_xlabel(r'$\alpha_{\mathrm{w}}$')
    ax.set_xlim(0.85,1.22)
    ax.set_ylabel( YLAB[i] ) 
    ax.grid(alpha=0.3)
    ax.xaxis.set_major_locator( matplotlib.ticker.MultipleLocator(0.1))
    ax.xaxis.set_minor_locator( matplotlib.ticker.MultipleLocator(0.05))
    ax.tick_params(axis='both', which='major', length=2.5)
    ax.tick_params(axis='both', which='minor', length=1.5)
    
    ax.axvline(1.136, color='0.2', alpha=0.5, lw=2.)
    ax.text( 0.76, 1.01,  r'$\alpha_{\mathrm{evap}}$', \
            ha='center', va='bottom', transform=ax.transAxes)
    
    
    
#axL[0].set_ylim(3.2,4.6)
#axL[1].set_ylim(0.8, 3.2)


#axL[4].set_ylim(0.,10.)
axL[0].legend( bbox_to_anchor = (0.05,0.95), loc='upper left', title='$w_{L}\,(\mu\mathrm{m})$',\
        numpoints=1, labelspacing=0.2, prop={'size':8},\
        handlelength=1.1, handletextpad=0.5 )




fig.savefig('alpha_opt/002.png', dpi=300)


From here on, tests....


In [7]:
import scubic
import lda

s=7.
g=4.35
wIR=71.
wGR=80.
T = 0.2
extents=50.

pot = scubic.sc(allIR=s, allGR=g, allIRw=wIR, allGRw=wGR)
test = lda.lda(potential = pot, Temperature=T, a_s=650., extents=extents, globalMu='halfMott', verbose=True, ignoreExtents=True)
print test.Number
fig = lda.plotLine( test , line_direction='111', extents=extents) 

print test.getRadius() / wIR


OK: Bottom of the band has positive slope up to r111 = 10 um
OK: Chemical potential is below evaporation threshold.
OK: Chemical potential is below the bottom of the band along 100
OK: Radial density profile along 111 decreases monotonically.
266875.838949
-8.66116121158 6.98591809963
0.280365579054

In [10]:
import scubic
import lda

s=7.
g=4.3479
wIR=60.
wGR=60.

pot = scubic.sc(allIR=s, allGR=g, allIRw=wIR, allGRw=wGR)
test = lda.lda(potential = pot, Temperature=0.26, a_s=650., extents=50., globalMu='halfMott', verbose=True)
print test.Number
fig = lda.plotLine( test , line_direction='111', extents=50.)


OK: Bottom of the band has positive slope up to r111 = 10 um
OK: Chemical potential is below evaporation threshold.
OK: Chemical potential is below the bottom of the band along 100
OK: Radial density profile along 111 decreases monotonically.
545520.541319

In [32]:
import optevap
wIR = 47. 
wGR = 40.
s=7.
results = []
results.append( [s, wIR, wGR]+optevap.optimal( s0=s, wL=wIR, wC=wGR, T_Er=0.383, a_s=650., extents=100. ))
print results


[[7.0, 47.0, 40.0, 4.006665678552408, 1.9545011948819773, 1659347.2424227342, 2.4501077005510155]]

In [34]:
import scubic
import lda

s=7.
g=4.0066
wIR=47.
wGR=40.

pot = scubic.sc(allIR=s, allGR=g, allIRw=wIR, allGRw=wGR)
test = lda.lda(potential = pot, a_s=650., globalMu='halfMott', Temperature=0.383, extents=100., verbose=True)
fig = lda.plotLine( test , line_direction='111', extents=50.)


OK: Bottom of the band has positive slope up to r111 = 10 um
OK: Chemical potential is below evaporation threshold.
OK: Chemical potential is below the bottom of the band along 100
OK: Radial density profile along 111 decreases monotonically.

In [37]:
import optevap
wIR = 32. 
wGR = 27.
s=7.
results = []
results.append( [s, wIR, wGR]+optevap.optimal( s0=s, wL=wIR, wC=wGR, T_Er=0.383, a_s=650., extents=80. ))
print results


[[7.0, 32.0, 27.0, 3.938097014635825, 2.1596415418707413, 473357.98910100898, 2.4784877360872914]]

In [39]:
import scubic
import lda

s=7.
g=3.938
wIR=32.
wGR=27.

pot = scubic.sc(allIR=s, allGR=g, allIRw=wIR, allGRw=wGR)
test = lda.lda(potential = pot, a_s=650., globalMu='halfMott', Temperature=0.383, extents=80., verbose=True)
fig = lda.plotLine( test , line_direction='111', extents=50.)


OK: Bottom of the band has positive slope up to r111 = 10 um
OK: Chemical potential is below evaporation threshold.
OK: Chemical potential is below the bottom of the band along 100
OK: Radial density profile along 111 decreases monotonically.

In [ ]: