JJsimulationinPython3/SimpleJJinpython.py at master · zcan11/JJsimulationinPython3 · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
# -*- coding: utf-8 -*-
"""
Created on Sat Mar 04 22:33:55 2017
This is a line to line translation from Erik's Matlab simulation to python code

%Explanation of base program

%The program splits a single junciton up into xmax discrete sections.  Each
%section has a supercurrent density, and a phase difference that is
%contributed to by the field (PhaseF), and an arbitrary phase set at x=1
%called Phase1.  For each value of field and Phase1 the supercurrent at
%each point is calculated and the net supercurrent (sum along the junction)
%is calcualted to find the supercurrent carried at that phase and that
%field.  To find the maximum supercurrent for a given field, the max of
%that vector is taken, and can be plotted against the magnetic field.  This
%is the exact measurment of the critical current vs field that we do in the
%lab.

@author: can
"""

###################################
# Import libraries
###################################

import numpy as np
import operator
import matplotlib
import matplotlib.pyplot as plt


###################################
# Define and initialize variables
###################################

xmax=51
xArray = np.linspace(1, xmax, xmax, endpoint=True)

SCurrentMag=1.0
SCurrentNoiseMag=0.
# instead of defining a vector with zero and fill them in the loop
# i jsut create a vector with incremental value here
SCurrentDensity=(SCurrentNoiseMag*np.ones(xmax)+SCurrentMag*np.ones(xmax))/(xmax)
pi_true=np.pi

pmax=101
Phase1Min=-0*pi_true
Phase1Max=2*pi_true
Phase1=np.linspace(0, 2*pi_true, pmax, endpoint=True)


fmax=1001
FluxMin=-3
FluxMax=3
FluxinJunc= np.linspace(FluxMin, FluxMax, fmax, endpoint=True)

Phase1MaxSC=np.zeros(fmax)

SCurrent=np.zeros((xmax,fmax,pmax))
SCurrentNet=np.zeros(pmax)
MaxSCurrentNet=np.zeros(fmax)
IndexMax=np.zeros(fmax)


# And so on...1

###############################################################################
# Loop doing the work, with additional detail explained in comments within the loop
###############################################################################

#   Most of your (well-commented) code goes here.

for i in range(0,fmax):

    PhaseF=2*pi_true*FluxinJunc[i]*np.transpose(xArray)/xmax;

    for j in range(0, pmax):
        SCurrent=np.transpose(SCurrentDensity)*np.sin(Phase1[j]+PhaseF)
        SCurrentNet[j]=sum(SCurrent)

    IndexMax[i],MaxSCurrentNet[i]=max(enumerate(SCurrentNet),key=operator.itemgetter(1))

    Phase1MaxSC[i]=Phase1[int(IndexMax[i])]


#########################################3
######################################
# End of loop. Print the results.
print("all done. Result is... ")
plt.figure()
ax = plt.gca()

plt.subplot(1,2,1)
plt.plot(FluxinJunc,MaxSCurrentNet,'o', markersize=2)
plt.xlabel("FluxinJunc")
plt.ylabel("Ic")
plt.title("Fraunhoffer Diffraction pattern of Ic vs Flux-Can Zhang")


plt.subplot(1,2,2)
plt.plot(FluxinJunc,Phase1MaxSC,'o', markersize=2)
plt.xlabel("FluxinJunc")
# this arbitrary phase is the phase that maximized the critical current
plt.ylabel("Arbitrary Phase")
plt.title("Arbitrary Phase vs Flux-Can Zhang")

plt.savefig("Bext_vs_Ic_and_B_vs_Phase.png")

print("all done")


###############################################################################


################################################################################
# Make informative comments on the results you just found, like how quickly
# the series converges, and its comparison to the true value of pi.
################################################################################