HHBL

M|

mm

1

T~-, /TV.:. ••* WV*-;_ ; <•<* .'. '.»*>

li
NBS TECHNICAL NOTE 627

Computation of

Spectral Data for a

Josephson Junction Circuit
NATIONAL BUREAU OF STANDARDS

The National Bureau of Standards1 was established by an act of Congress March 3, 1901. The Bureau's overall goal is to strengthen and advance the Nation's science and

technology and facilitate their effective application for public benefit. To this end, the Bureau conducts research and provides: (1 ) a basis for the Nation's physical measure- ment system, (2) scientific and technological services for industry and government, (3) a technical basis for equity in trade, and (4) technical services to promote public safety. The Bureau consists of the Institute for Basic Standards, the Institute for Materials Research, the Institute for Applied Technology, the Center for Computer Sciences and

Technology, and the Office for Information Programs.

THE INSTITUTE FOR BASIC STANDARDS provides the central basis within the United States of a complete and consistent system of physical measurement; coordinates

that system with measurement systems of other nations; and furnishes essential services leading to accurate and uniform physical measurements throughout the Nation's scientific community, industry, and commerce. The Institute consists of a Center for Radiation Research, an Office of Measurement Services and the following divisions: Applied Mathematics—Electricity—Heat—Mechanics—Optical Physics—Linac

Radiation2—Nuclear Radiation 2—Applied Radiation 2—Quantum Electronics3—

Electromagnetics''—Time and Frequency3—Laboratory Astrophysics3—Cryogenics3

. THE INSTITUTE FOR MATERIALS RESEARCH conducts materials research lead- ing to improved methods of measurement, standards, and data (_a the properties of

well-characterized materials needed by industry, commerce, educational institutions, and

Government; provides advisory and research services to other Government agencies;

and develops, produces, and distributes standard reference materials. The Institute consists of the Office of Standard Reference Materials and the ff Mowing divisions: Analytical Chemistry—Polymers—Metallurgy—Inorganic Materials—Reactor

Radiation—Physical Chemistry.

THE INSTITUTE FOR APPLIED TECHNOLOGY provides technical services to pro- mote the use of available technology and to facilitate technological innovation in industry and Government; cooperates with public and private organizations leading to the development of technological standards (including mandatory safety standards), codes and methods of test; and provides technical advice and services to Government agencies

upon request. The Institute also monitors NBS engineering standards activities and

provides liaison between NBS and national and international engineering standards

bodies. The Institute consists of a Center for Building Technology and the following

divisions and offices: Engineering Standards Services—Weights and Measures—Invention and Innovation—Product Evaluation Technology—Electronic Technology—Technical

Analysis—Measurement Engineering—Fire Technology—Housing Technology4 —Federal Building Technology4 —Building Standards and Codes Services4—

Building Environment4—Structures, Materials and Life Safety4—Technical

Evaluation and Application4. THE CENTER FOR COMPUTER SCIENCES AND TECHNOLOGY conducts re- search and provides technical services designed to aid Government agencies in improving cost effectiveness in the conduct of their programs through the selection, acquisition,

and effective utilization of automatic data processing equipment; and serves as the principal focus within the executive branch for the development of Federal standards for automatic data processing equipment, techniques, and computer languages. The Center

consists of the following offices and divisions: Information Processing Standards—Computer Information—Computer Services —Systems Development—Information Processing Technology.

THE OFFICE FOR INFORMATION PROGRAMS promotes optimum dissemination and accessibility of scientific information generated within NBS and other agencies of

the Federal Government; promotes the development of the National Standard Reference Data System and a system of information analysis centers dealing with the broader

aspects of the National Measurement System; provides appropriate services to ensure

that the NBS staff has optimum accessibility to the scientific information of the world,

and directs the public information activities of the Bureau. The Office consists of the following organizational units: Office of Standard Reference Data—Office of Technical Information and

Publications—Library—Office of International Relations.

' TIi ' rters and Laboratories at Gaithersburg, Maryland, unless otherwise noted; mailing : Wasl i :on, D.C. 20234.

2 Part of the Center for Radiation Research.

' Located at Boulder. Colorado 80302.

1 Part of the Center for Building Technology.
91S

fc'H Computation of Spectral Data

62- for a Josephson Junction Circuit

i • ©

E.G. Johnson, Jr. Quantum Electronics Division

D.G. McDonald

Cryogenics Division

Institute for Basic Standards

, National Bureau of Standards

Boulder, Colorado 80302

NBS Technical Notes are designed to supplement the

Bureau's regular publications program. They provide

a means for making available scientific data that are

of transient or limited interest. Technical Notes may

be listed or referred to in the open literature.

U.S. DEPARTMENT OF COMMERCE, Peter G. Peterson, Secretary

NATIONAL BUREAU OF STANDARDS, Lawrence M. Kushner, Acting Director

Issued November 1972
National Bureau of Standards Technical Note 627

Nat Bur. Stand. (U.S.), Tech. Note 627, 63 pages (November 1972)

CODEN: NBTNAE

For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington. D.C 20402

(Order by SD Catalog No. C13.46627).
Table of Contents

Page

1. Introduction 1

2. The Basic Equation System 3

2. 1 Definitions 3

2. 2 Equation System 4

3. The Approximations Needed for Accurate Numerical Results . 7

3. 1 Constraints 7

3. 2 Part (a) 8

3. 3 Calculation Sequence 8

3. 4 Problems 9

3.5 Part (3) 10

4. The Fortran Program and Flow Chart 11

5. The Sample Collection of Output 11

5. 1 The Printed Output 11

5. 2 The Microfilm Output 23

6. The Summary 46

7. Appendix 49

7. 1 The Transient Equation 49

7. 2 The Use of the Transient Program 50

8. Acknowledgments 55

9. References 55

in

Computation of Spectral Data for a Josephson Junction Circuit

E. G. Johnson, Jr. and D. G. McDonald

A computer program has been developed to study power flow between different

frequency channels in a Josephson junction circuit. This paper discusses the mathematical assumptions used to get such results. They are the trapezoidal approximation from

spline theory and the use of a finite range of frequencies to characterize the frequency

spectrum. This paper describes the program and provides the FORTRAN listing, flow

charts, and discusses how to use the program. A discussion of possible sources of errors is also included.

Key Words: Differential equation; fast Fourier transform; Josephson junction; nonlinear

integral- differential equation; spline theory.

1. INTRODUCTION

The central problem of interest is studying the power flow between various frequency channels

in a Josephson junction circuit that is appropriate for harmonic mixing. To be realistic, the circuit

must include shunt capacitance and shunt resistance as well as source andload impedances all of whichmay

be frequency dependent. The general circuit we have chosen to investigate is shown in figure 1.

If the impedance Z., of figure 1, are not frequency dependent, it is relatively easy to solve the

differential equation for the circuit and this solution is included as an appendix. The result is what we

term the transient program, i. e. , it gives the behavior of the circuit as a function of time starting

from specified initial conditions. After a few oscillations, the transient behavior reaches a good

approximation to a steady state oscillation and thus provides an estimate of a consistent set

of currents, of currents, voltages, and phases which we use in the main program, the steady state program.

The steady state program imposes strict periodicity on the solution as required for Fourier

harmonic analysis. Steady state solutions are generated in two steps, first with the Z. frequency

independent and second, using the full integro-differential equation for general Z.. In the notation of

Josephson junction when we are given the applied signal V . The V is specialized to be

figure 1, we want the steady state power for each frequency present in the detector, Z_, and in the

.re given the applied signal V . The V is speci

V b + b T cos (NT s + 9 T ) + bTT cos (Ns + 9 TT ) O I O L Li L U U U J V = Vjb_ + b T cos (NT s + 9 T ) + bTT cos (NTTs + 9 TT ) | . (1)

Here V is the voltage scale factor, b is the dimensionless d. c. voltage, b is the lower applied r. f. O O L

amplitude, 8 L

is its phase, and b^ 9 y

are the same quantities for the upper applied r. f. . The integers,

"frequencies" NJj and N generate multiples of the phase s, o< s - W t ^ 2tt where — is the base period u o uu c

of the system. Since our strategy is to perform a discrete Fourier analysis, obviously we must adjust

the discrete spectrum to be consistent with the applied frequencies NT and NTT , the rotation frequency

2eV Nj of the self oscillation W = —— = N U) and all possible linear combinations of these. The b is n i o o adjusted to allow only periodic solutions with each frequency related to this base frequency, U) , by an

integer. Note that b

Q

is a function of b^ h^, 9^, e a nd N , as well as the impedances Z..
O UJ LU t—

I— OO

UJ >- Q CO 3ouUao

COA

a>

COOu
This report is subdivided into six additional sections. (2) The scale adjustment is defined and

the basic equation system is determined. (3) The approximations necessary to solve the basic equation

system are given. (4) The FORTRAN Program and flow chart used to solve the equation system are

described. (5) We present a sample collection of output from the computer showing the following for each case:

(a) the time plot of i , (b) the magnitude of the frequency spectra for i , (c)the magnitude of the frequency

spectraforV , (d) the relative phase of (i - V ) at each frequency, (e) the amplitude spectrum of -—, (f) the

out i. out as steady state power into the detector Z , (g) the steady state power across the junction device, and

finally (h) a summary chart showing such things as the scales of amperes, volts, resistance, frequency,

the b , b , b , 9 , 8 and the power and other items at the following selected frequencies, namely

d. c. , N , N , the direct beat frequency N , and the rotation frequency N . (6) The summary gives

general conclusions and recommendations. (7) Finally, we have the appendix describing the transient

program.

2 . THE BASIC EQUATION SYSTEM

2 . 1 Definitions

The following definitions are noted:

Z. are the frequency dependent impedances for the circuit. R is the shunt resistance of the Josephson Junction.

Z represents any residual impedance of the Junction device.

C is the capacitance of the Josephson Junction.

cp is the Josephson phase function.

J is the Josephson critical current. V is the voltage scale,

o

&

W<#,(uu )' is the radian frequency scale. R is the resistance scale,

o A is the current scale,

o P is the power scale. K is the impedance operator.

K is the d. c. part of the impedance operator.

D is the integral equation operator in the circuit system.

U is the net d. c. current needed to drive the system at the steady state rotation frequency N .

a is the characteristic radian frequency squared for critical current in dimensionless form.

ANU (v) is the characteristic frequency for loss due to RC type effects. ANUP (v ) is the characteristic frequency for loss due to the rest of the P

circuit. ANUT (v _) is the total loss characteristic frequency.

V, VA (v, v ) are the various voltages in dimensionless form.

*(Note that the first symbol is the FORTRAN one and the second is the algebraic form. This is only

shown for the less obvious cases. )
PHA (cp ) is the phase with the rotation and directly induced r. f. signal

removed.

PHB (cp ) is the directly induced r. f. signal phase plus the rotation and

the initial phase at s = 0. PHC (cp ) is the initial phase at s = 0. c THL, THU (9 , 9 ) is the phase of the lower, upper r. f. frequency.

2. 2. Equation System

We define

S = iu oo2e t, V = -r-^ , Aoooiio = V /R , Z. = z.R , J = A a, o

C = c,/R looU) , R = R r,, i, = A I, (<t = 1 to 6), and UO = 2tt/ T

o 1 t, o <C o

(2)

Note that a, R , and (JJ are fixed by requiring c to be unity in the dimensionless system, by requiring ar = — , and by using a near 1. This allows a good initial guess for starting slope and d. c. current

1 1 in the program. Near the hysteresis region a 1 would have to be carefully adjusted to allow study of

this region.

Further definitions are:

V = =— -rr- - y -r- = V V, i = A a sm cp , J 2e dt o ds o 4 o T

dV, dV, ... VT A , 2 „ „ 2 . dV , . L o dcp i_ = C -r—- = V C -=— Ul = A -— , and l = -— = -?- 5 dt o ds o o ds 6 R r ds

T . V dVl

Therefore l = 1 . + i_ + l, = A a sm cp + — + -— 145 6 ol r ds

J Now V s y V., i = i. + i_, i, = A I,, i. = A I. , s o 1 1 2 2 o 2 1 o 1

i = A I , V , = V V . o o out o 2 Therefore we have the dimensionless operator equations:

V*3 '[¥-]•

dV

I, = a sin cp + VI r. + -r— , 1 1 as

I =(l+z

3

- 1 Z2 )I1+ z

3

- 1 V,

V

1 = Z

1 I + V

2 ' and

V

2

= Z

3

J2

= Vl + V •

(3)

(4)

(5)

(6)
(7)

These are reduced to a single operator form.

Vj = Zj (1 + z

3

" z

2

) I

1 + ZjZ

3

" V + z 2

I + V, or

v

i = [

z

i + zz + Vs'S] [

asincp+ drJ

+ {[

Z

1 + Z

2 +Zl Z 3" lz

2]

/r

i +Z

l Z3

" 1+ ! J V -

We finally define

K= Zj+ z

2 + Sj^" 1 ^. V* 1/rj, DHK

'1(z

1 z

3

" 1 + 1) - ^- (*f(z°)~ * + l)> and Vp

^(z^ )" 1 + l) o o

< 8) as the d. c. part of operator, where z. is the D. C. value of z., and

v T

= v+ v p

to get the following pair of operator equations

-1 dV K (V n ) = a sin cp + —-+ v^V + DV, and

1 as 1 v.*.

<9) ds

As the reader can see, if the z. are frequency independent then the D operator is zero and we have a pure differential equation. If, as would usually be the case, there is frequency dependence, then D

4 <

s., namely z ~ 10 , z ~ 1, and z_ 4 < „8is not zero. However, if we put in estimates for the values of z., namely z ~ 10 , z ~ 1, and z_ ~ 10 , we have D ~ 10 . If a ~ 1, then our problem is dominated by the differential equation and the integral

equation part can be used in an iteration procedure. The acid test is to run the program.

As we are interested only in strictly periodic solutions, we assume the following forms for the Cp dcp dcPb and V. First we split cp = cpT + cp,; V = V + V, and note that cp =Ns+cp +cp, and V, = N, + -rp-d s . a b a b bled b 1 ds ds

We require cp to be defined by

i d ^H dCPH

dCp ^ K (v iA.c»^ + \-dT +D ^' (10)

d s and cp to be defined by

a dV

u = —;-2-+ a sin (cp + cp, ) + v V + DV . (11) ds a b T a a

Here u = K b - V^N, . o o 11

The purpose of this split is to remove a ramp term and to get the first order induced phase function, so

that the remaining part, cp V are periodic and have a chance for rapid convergence to the correct solution.

Since V , cp , cp are all periodic function of s, the exact solution would be periodic. We use the

a a d

Fourier space form to get an idea of the structure of D. Thus with -a> s t ^ m,
cp (s) = ,S e' lU

cp (*.), V (s) = ,2 e"

1,t/S a-o aa-o V {I), (note that V (0) =0.)

a

* a'

-iN s iN s " iN

Ti s iN

TT S

cp d (s)-e

L

5

d

(NL)+e 9* (N^) + e V/^+e U

^(Njj) -, (12)

., 2lT , and D(s) = E e*

S D(-t). ("We have made DV = f

-j- D(s - s') V(s').)

The equation (10) defines cp as

1 b L

_i9 L K (NL

)_T

e = ( " iNL

)+ V T

+ ° (NL

) ( " 1NL) {pd

tNL

) " (13)

The cp (N ) is gotten by U replacing L in the above expression, cp is the complex conjugate. We note

d U

that N and N cannot be zero. U i-i The equation (11) can be written as

K ? ^ d if u6 = (-it) + v(-il)+D(l) (-U) cp M)+af -^ e

1* 8

sin(cp + cp ) . (14) o T a I •'o

2tt """a " ^b'

where 6 is the Kronecker delta. We further note that D(-t) is given as

5(1) = (z

1 (l)z

3

~ 1 (l) + l)/KM) - (Zl (o)z

3

_1

(o)+ l)/K(o),

where

KM) = z (t)+ z (l)+ z (t) z At) z(l)~

l and

K(o) = K . o

It appears to be best to do the calculation in the time domain and then do the spectral transform.

Thus we have

r~ - iNL

s

cpd(s) = 2 Real ^(N^) e "iN s]

cp, (s) = cp (s)+ N s + cp , b d 1 c

d

/" 2TT d ' u = — V, + a sin (cp + cp (s)) + V_V +/ — D(s-s') V (s')

,

ds a ah T a •' o 2rr a

(15)

and v

a

= a^- ^a

as our basic equation system. To progress further, we must make numerical approximations.
3. THE APPROXIMATIONS NEEDED FOR ACCURATE NUMERICAL RESULTS

Wehave to choose a discrete set of s-points, N , because we must take a finite amount of computation time and must keep the entire data set in the machine to allow fast fourier transformation.

N

The set of points has to be equidistant and equal to 2 °, where N is some integer. Practical considerations restrict the range of N from 8 to 11.

3. 1 Constraints

A discrete set of points implies frequency folding, hence we require the following conditions

to minimize this effect. ..N

2

(a) The highest frequency of interest must be ^ — . All frequencies above this value are to

8 be treated as unreliable, due to the method of analysis. Computational noise must also be small.

(b) Our demand accuracy for analysis is near 5%. We make checks on the results by increasing N as the occasion demands. If the spectrum does not change to the demand accuracy, our analysis is treated as adequate and accurate to that order. Incidently, if the value of U is correct to 5%

or less, then the lower frequency spectra is also correct to that accuracy.

(c) This 5% computational accuracy level suggests using N_ = 256 or 512 with an occasional

extension to 2048 points for special studies. ,. 2

(d) We allow z.(-t-) to have arbitrary values for s l^——and require the remaining t values to be

fixed by z.*M) = z.(N - I).

(e) In D(-t) and cp ,(<£,), we use inplace of (-i-C) the periodic expression(e - l)/h, where h= 2rr/N is d <2 the interval between points in the time space. This allows folding and preservation of reality conditions. The sensitivity to this change is checked by the doubling of points.

(f) Because D(s - s

1 ) is a small correction to the differential equation system, we replace,

with no qualms, the integral by a sum using the trapezoidal rule. Thus

N

2 -l rds'D(s-s')V(s')= E rr D(j - j') V (j

1 )- (16)

/ ~ST

a

j'=0 N

2 a

•'o dV

(g) To complete the approximations, we must decide how to approximate —— type terms.

Finite differencing does not allow high enough precision in the fit. Since the spectral information for

dV 2 power uses —— in its current, we need to have smoothness in the results to at least order h . , Additionally, we want to avoid absolute instabilities and permit the application of quasilinear theory

to get the b value needed to allow strict periodicity. These requirements make it necessary to use

°

2 the spline method using at least the trapezoidal rule. We assume the trapezoidal rule everywhere in

the calculation. The results of this analysis should be smooth functions (accuracy order, h ), hence

a frequency separation of true results from noise due to the discrete analysis of the system and an

2 inherent overall accuracy of h . Periodicity should tend to force uniform error throughout the time

interval. Sharp spikes in the time dependence of cp wm negate the uniformity of the error and make

higher frequencies spectral information less accurate.

(h) Now the problem is split into a part (a) which neglects D and a part (B) that takes D into account.
3. 2 Part (a)

Here we have the simplified system of equations using the trapezoidal rule.

cp (j+ D = cp (j) + T V (j+ 1)+ V (j)

v

a (j+ D = v

a

(j)+ - |q,(j+ D+ q

|[va

(j+ i) + v

a (j)]

q^) = -a sinlcpa

(j)+ cpb

(j) - V T

V

a

(j)+ u.

V (o) = V (N,), and cp (o) = cp (N).

a a 2 a a 2

The u is adjusted to allow this result. The eq (19) can be rewritten with definitions

hj = (2 - hv

T

)/ (2 + hv

T

) , h

2

= - ah/ (2 + hv

T

) , h

3 = 2h/ (2 + hv T

),

The boundary conditions are

V

a

+ 1) = h

1 V

a

(j)+ h

3 u+ h I sincp(j + 1)+ sincp(j)

(17)

(18)

(19)

We generate two solutions starting with cp (o) fixed and unchanged for a given calculation run. To avoid

numerical noise problems, experience shows that it is absolutely necessary to fix the cp (o), when u is a

adjustable.

3. 3 Calculation Sequence

<P (j + 1 ) - cp (j) a » a

(1) We iterate for a given j, j+ 1 pair until of the n+ 1 iteration changes by

less than EPI =10".

(2) We compute two quasilinear solutions using the following starting conditions

6cp (1)

(o) = , 6V (1)

(o) = 1 , and 6cp (2) (o) = 6V (2) (o) = 0. a a a a

In case one 6u = and in case two 6u = 1. The two equation sets are, with q (j) = cos cp(j),

(i) ( 6cp

l

' (j + 1) = Sep

v a a

and

(i) (i), 6V 0+ D= h 6V x '(j)+ h. 6 a la 5 c.

i,

(J) + |[5Va

(1,

(j + l) + 6V

a

(»U)],

q

2 (j+ D6cp

(i) 0+ D+ q

2

(j) 6cp

(%]

(20)

These can be rewritten as

6cp (j + 1) = 6cp

(l) " a a

i)

(j) + |[6v

a

(i) a + i) + 6v

a

(i) a)].
6V

a

(i) (j+ l)=T

t

(i) /T

2

(i)

with ? l

(i) * [

h

l + h2l q 2 0+ U

J 6V

a

(l) ^ )+ h

3

6 21+ h

2

5cP

(l) <J> [q 2 0+ D+ q

2 0)l. (2D

(i) h

and T^'- l-h2l q

2 0+l) .

(3) We use the values computed at cp (N ), 6cp (N ), etc. to satisfy the following requirements and hence fix P and 6u:

cp (N ) - cp (o)+ 6cp (1) (N,)+ 6 U 6cp (2) (N, ) = 0, a 2 a a 2 a 2

and V (N ) - V (o)+ 0(6V (1) (N ) - 1)+ 6u 6V (2) (N ) = 0. (22) a 2 a a i. a 2

The new cp (j) , V (j) » an<* u are now defined by

a a

cp 0)= cp 0) (old) + 6cp (1) (j)+ 6u6cp

(2)

(j) , a a a a

V (j) = V (j) (old) + P 6V (U

0) + &u 6V (2)

(j) . (23) a a a a

and u = u (old) + 6u .

In the computer program the corrections to the initial guesses of u and V (o) are bounded because

initially the first guess may not be in the circle of convergence except when the computer program in the appendix is used. Thus, it is necessary to have built into the calculational procedure a means to search for the region where the quasilinear technique will work. The computer program has allowed

jumps in 6u and in 6 V (o) that shrink automatically as each computation is made with a new u and

a V (o) . If after ten iterations the convergence has not taken place, then the V (o) is set equal to

3. 3. V (N ) using the idea that it is more likely to be compatible with a periodic solution because the transient should be less. Then we repeat the above quasilinear process and the window shrinking

process for up to eight times. If the convergence is particularly difficult, the shift in V (o) can

take place up to ten times. Finally, a relaxed test for acceptance of the results as converged is made

so that a spectral analysis can be made.

3. 4 Problems

There are two main sources for convergence problems: (1) The number of computed points

is too small when there are sharp spikes in the V and -:— . This causes extensive computational

noise. (2) The initial V (o) and u are badly guessed. Both problems are particularly bad when

there are strong B , B terms and/ or when a, v are very large numbers.

V J-i T
3.5 Part (P)

3 Here we use the full quasilinear theory. First, we define

N2 -l

q (j)= -a sincp(j) - £ D(j - j') s V (j ') + u - v T V (j)

J = 2

q

2

(j) = coscp (j)

q

3 0) = h3H (Va^ + 1)_ V

a

0) >- °' 5 (^ i

<J + D+ q^))' (24)

q4 0)= hji (cpa

(j + 1) -q>a (j» " 0-5 (V

a

(j+ D+ V

a

(j))|

We note that q_, q are the errors due to not fitting the equation system to a quadratic spline. We

use these terms to define the error signal in the quasilinear theory. In addition to the terms needed

to correct the error signal, we must add effects due to needed changes in V (o), and u to allow pea

riodic solutions. Since the quasilinear theory integrates over the entire domain, the u is sensitive

to the global solution, and periodicity tends to force global accuracy; we should retain a global ac2 curacy of order h , and a smooth spectrum to that same order. This should be accurate enough to separate the true solution from the effects of the finite analysis. Again, spiking will make the higher

frequencies less accurate. The quasilinear differential equations are replaced by the trapezoidal rule. The basic equations become

. ( 6cpa

and

°G + D= 6cp a

(l) 0)+|["6V

a

(i) (j)+ 6V

a

(i) (j+ 1)1 - q

4

(j+ D *\

6V (i) (j+ D = h. 6 V

(l) (j)+h, [q,(j) 6 ^

(l) (j)+ q,(j+ D 6cP ^O + 1)1 ( 25 >

+ h

3

6

x2

- q

3 (j+ 1)6^ , for i= 1,2,3.

The initial conditions are 6cp ^'(o) = 0, for all i; the 6v l1 ' = for i = 1, 2, and the 6v (o) = 1. a a a The periodic conditions at j = N requires

6cp (1) (N ) + 6u 6cp (2) (N ) + 3 6cp (3) (N

2

) = 0,

6 V

(1) (N ) + 6u 6V (2) (N ) + P (6V

(3) (N ) - 1 ) = 0. (26)

c. 2 «

We use these equations to fix 6u and P. We compute for the different a, a new cp =

old cp + a(6 cp + 6u6 cp + P6 cp ) , and a similarily constructed new V which are then inserted

a a aa aN2-!

2 2 to compute a new q and q set for each a. The sum of these; namely, P(a) = N^ (q + q ), 3 4 jVo

gives us an error signal. Our choice for the correct cp , V is fixed by the requirement the correct

a a a has; P(a) <P(o). The smallest P(cp) improves the calculation the most and guarantees convergence.

10
The a chosen, is then used to construct the new cp and V and the entire part (3) quasilinear analysis

a a

2

is repeated. This process is continued until 10 iterations have been performed or until P(o)< 2h . The results of this process constitutes the best accuracy that can be obtained for a global solution

subject to the inherent errors of the trapezoidal method.

Next we compute by fast Fourier transform the desired spectra. This completes the discussion

of the assumptions made to get numerical results. The next section has the FORTRAN program listing

and a flow chart of the basic program flow. In particular, we show the input data requests.

4. THE FORTRAN PROGRAM AND FLOW CHART

Here we have listed the Fortran programs that give the main details of the computation. The

microfilm and fast Fourier transform programs are not listed. The former depends on the computer

installation and the latter should be available and optimized at each computer installation. After the Fortran list, figure 2, a flow chart, figure 3, shows the bulk data and computational flows between

programs.

5. THE SAMPLE COLLECTION OF OUTPUT

This section discussed the printed and the microfilm output for a single computed case. The

purpose is to give general information to allow use of the Fortran program by an interested party. No

attempt is made to report all the possible data that exists as a consequence of the computations.

Selected information generated as a consequence of this program will be described in separate reports.

To use the program correctly in detail, one must read the Fortran listing.

5. 1 The Printed Output

Figure 4 shows the printed output from a single computed case. The first row of data gives

the number of data points used, the number (11) corresponds to 2048 points in time and 1024 different

frequencies. The J (1. -4) is the supercurrent in amperes. The Al (1. 25-0) gives the scale adjustment

defined on page 4. C (8. 2277-15) is the capacitance in farads. R (1. ± 1) is the number of ohms.

EB1 (1. -8) is the demand accuracy for each step interation and scales the overall demand accuracies

in thevarious iterations of part (a) and part (3). The EP6 (1. -0) can be set to increase demand accuracy

of part (3). EP6 equal 1. uses the values stated on page 10 for the converged exit. EP6 less than one increases the accuracy demand and EP6 greater than one decreases accuracy demand.

The first Zl, etc. , row gives, with 19 equal 1, the complex values of Zl, Z2, and Z3 for all impedances as if they were frequency independent. Subsequent rows of Zl type data show specific

frequency exceptions to the frequency independence assumption. These are overlaid on the previous

values. Please note that 19 is one unit larger than the corresponding frequency number, thus 1 is the

d. c. case. In our sample output we have changed the first harmonic of the Z3 impedance from (1. + 8, 0. ) to (1. + 5, 0. ). The others have remained unchanged.

We continue the row by row analysis of the printed output. After the impedance spectrum is

defined, we next specify the lower frequency, NL(4), the upper frequency, NU(5), at which the two r. f

11
PROGRAM JOSPHC

C THIS PROGRAM COMPUTES FREQUENCY SPECTRA FOR A CIRCUIT MODELING

C THE POINT CONTACT JOSEPHSON JUNCTION. THE ANALYSIS USES A TIME

C PERIOD WITH N2 POINTS. NO IS THE ORDER 2**N0=N2 NO.LE. 11 . DATA IP 1=3.141 592653 59) .( IM=(0.,1. )

) , ( PHS=2.06778E-15

>

DIMENSION ZK1025) »Z2(1025) »Z3(1025)

COMMON /5/ BI»S1»S2»S3

1»A1

COMMON /16/ V16.P16.H8

TYPE COMPLEX Zl ,Z2 »Z3 .K ,DT ,ZF »ZH»ZG »KL »AN2» ATMP.KU.ATN.AN6.AN7.ZI

1»CBU»CBL »ZD TYPE REAL KO

COMMON /DAT/ PI.IM.PHS.U

DIMENSION K(1025> ,DT(1025> .KL(3),KU(3>

1 »I2T(1025> »VAT<1025>

EQUIVALENCE (K( 1) ,DPH(1 ) ) » ( DT ( 1 ) ,DVA(1>)

COMMON /l/ PHC»VA(2050).DPH(4»2050>»

1DVA( 4.20 50) »H,H1»H2.H3.H4,EP1.EP2.EP3.LBX(5> ,LBTR(5) .EP4(6)

EQUIVALENCE (DPH(6151) .ZKlD * ( DVA ( 6151 ) *Z2( 1 ) ) » ( Q4< 1 ) »Z3( 1 ) ) COMMON /2/D(2050) ,PHA(2050) .PHB(2050> »QK2050) *Q2(2050) .03(2050) » 1Q4(20 50> ,KEY(6) .NO ,N1 ,N2 »NH ,NP ,P0 ,V0» AO ,V» THL .THUtBO.BL »BU.NL ,NU

2.NB.KL.KU. AN5.H5, ANU , ANuP > ANUT »A .Rl .RO ,C. WO.T , R, PI2 »AN2 ,AN4»

3AN3»MM.N5»N6.N7»ISKPD»SC» OUT < 6. 11 ) »NR. ICONV.AJ. ICASE.

4POW1.POW2.POWJ.POW I N.K0.EP5.H7

5.CBU.CBL

EQU I VALENCE < X(l) ,DPH(1 ) )»(Y(1) »DPH(2051) ) . ( VAT( 1 ) »DVA( 1 ) )

»

1( I1T(1 ) .DVA (2051) ) .(VOUT(D .DVA (4 101) ) »(X1(D ,DPH( 1025) ) TYPE COMPLEX VAT. I IT »VOUT. I2T .IM

DIMENSION XQ025) . Y ( 1025 ) » I IT ( 1025 ) »VOUT(1025) ,X1 ( 1025 ) »Y1 ( 1025 ) EQUI VALENCE (Y1(1),DPH( 4101) ) , ( I2T ( 1 ) .Q3( 1 )

)

1 READ 2.N0.AJ.A1.C.R.EP1.EP6

2 FORMAT* I5.6E12. 4)

PRINT 200.N0.AJ.A1.C.R.EP1.EP6

200 FORMAT(* NO » J .Al »C »R »EPl .EP6 */I5»6E12.4>

N2= 2**N0 tNH=N2/2 $AN3=N2 SPI2= 2.*PI $AN4=PI2/AN3 $N5=N2+1

N6=N2+2 $N7=N2-1 $AN2=-IM*AN4 $MM=N0~1 $SC=.5/AN3 $NP=NH+1

XT=PI2*C S R7=PHS/(XT*AJ) $ A=R7/ ( R*A1*. 1 )**2

R0=SQRTF(A*R7) $T=R0*XT *W0=PI2/T $NR=NH/2

Rl= R/RO J ANU=1./R1 $VO=PHS/T $A0»V0/R0 $PO=VO*AO $AN8=Nh

EP2=EP1*A $EP5=.1*A

S1«R1*SQRTF(A) $S2=A1/10. $ S3=EXPF ( .2-Sl

)

C SET I8 = I9«=1 FIRST TIME TO GENERATE ALL THE SAME Z1.Z2.Z3.

C TO ALLOW SELECTED FREQUENCY DIFFERENCE USE 18=0 19 THE PARTICULAR

C HARMONIC. TO EXIT USE I8=-l AND 19 THE LAST PARTICULAR FREQUENCY

C HERE I9=NP AND Z1.Z2.Z3 ARE REAL AS MUST BE Zlll) ETC.

C 19=1 IS D.C. CASE ADD 1 TO HARMONIC ORDER TO GET CORRECT

C LOCATION IN DATA STORE.

5 READ 3.I8.I9.ZK 19) .Z2( 19) »Z3( 19) SPRINT199. 19 »Z1 ( 19) .Z2 ( 19) »Z3(I9)

3 FORMATf I3»I5»3C(E12.3»E12.3)

)

IF( 18)4,5.6

6 DO 7 I=2»NP $ZKI)=ZKlJ $Z2( I>=»Z2(1> $Z3(D=Z3(1)

7 CONTINUE $G0 TO 5 4 DO 8 1 = 1, NP $Z1( I)=Z1( IJ/RO $Z2 ( I ) =Z2 ( I ) /RO $Z3 ( I

) =Z3 ( I ) /RO

ATMP = 1.+Z1( D/Z3( I) SK( I )=Z1(

I

>+Z2< I )*ATMP $DT ( I ) = ATMP/K( I

)

8 CONTINUE SANUP=DT(1) $ANUT=ANUP+ANU $D0 9 I=2.NP

9 DTI D=DT( I)-DT(l) $KO = K(D $ DT(1>=0.

READ 10. NL.NU.ISKPD. ICASE

C NL IS LOWER APPLIED FREQUENCY NU IS UPPER

10 F0RMATUI5)

C ISKPD=+1 DO THE D STRUCTURES. ISKPD=0 SKIP D STRUCTURES HERE D=0.

PRINT 201.NL. NU.ISKPD. ICASE

201 FORMATJ* NL »NU, ISKPD. ICASE *»4I5>

C ICASE IS INCREMENTED TO LABEL EACH COMPUTED CASE LOAD LAST VALUE

TF(NL-NU) 11.12.12

12 PRINT 202»NL.NU SCALL EXIT

202 FORMAT ( * NL TOO LARGE * .215)

Figure 2. The steady state Fortran program listing.

12
11 IF(NU-NR>13»14.14

14 PRINT 203.NU.NR

203 FORMAT ( * NU TOO LARGE * 215)

13 NB=NU-(NU/NL)*NL

KEY(1)=0 SKEY(2»=NB $KEY(3»=NL $KEY(4)=NU $KEY<5>=NH

IF( IS<PD>15»19

15 DT(NP)=.5*DT(NP) $ZH=CEXP ( AN2 ) $ZG=CONJG(ZH>$DO 16 J-1.N2 $D(J)*0«

ZG=ZG*ZH $ZF=1. $DO 17 I=2»NP $ZF=ZF*ZG $ZD=DT<I)*ZF

17 D( J)=D(J)+REAL(ZO)

16 D(J)=D(J»/AN8 $ DT(NP)»2.*DT(NP)

READ 18.IDRED $60 TO ( 19.20 > IDRED

C DECIDE PRINT OF D IDRED=2 YES -1 NO

20 PRINT 204»<D(J) ,J»1.N2)

204 FORMAT ( * D SET * /(10E12.3>>

19 CONTINUE $H=AN4 $EP3=2.*H*H*EP6

H4=H*.5 $ HT=H*ANUT S HU=1./ ( 2.+HT > $H3=2.*H*HU $H2«=-A*HU*H

H1=(2.-HT)*HU »H9=H4*H2 $ H7=HU

AN7= AN2*flU $ATMP=(CEXP<AN7)-1.)/H $ATN=( ATMP+ANUT+DT(NU+1 )

KU= 0.5/(ATMP*K(NU+D*ATN) $KU( 2 >=ATMP*KU $KU( 3 ) =ATMP*KU< 2 ) AN6=AN2*NL $ATMP= ( CEXP ( AN6 > -1. ) /H $ATN= ( ATMP+ANUT+OT (NL+1 )

KL = .5/(ATMP*K(NL+1)*ATN) $KL ( 2 ) =KL*ATMP SKL (3 )=ATMP*KL< 2 > C Nl IS RAMP FREQUENCY RELATIVE TO BASE USE PHC=-PI IS Nl NOT ZERO

C PH0=PHC AT S=0. BL.THL IS LOWER APPLIED RELATIVE VOLT AND PHASE

C BU »THU ARE THE UPPER

CALL CHARTA

C NOTE BL=K0*BL USED IN TRANSIENT PROGRAM. SAME FOR BU. HERE Ko IS

C Zl/RO. HERE R0=R*A1*A/10.

2022 READ 2*N1 .PHC.BL, THL»BU»THU

BL« KO*BL $ BU=KO*BU

C ADD EACH TIME ICOuNT »1 IGUS=-1 FOR SMALL POWER INCREMENTS.

C ADD ICOUNT =1 IGUS=1 AND USE RESULTS OF TRANSIENT PROGRAM.

C ADD ICOUNT =1 IGUS=0 TO START NON RF GUESS SEQUENCE.

C ADO ICOUNT =0.-1 TO GET NEW POWER VALUES

2020 READ 2041 » ICOUNT» IGUS.BI »V16»P16

2041 FORMAT(2I5.3E20.9)

IF(EOF»60>22»23

22 CALL EXIT

23 IF( ICOUNT12022. 2022,2023

2023 AN5=H*N1 $ KEY(6)=N1

S4=N1*S2

PRINT 2111»N1»PHC,BL»THL»BU»THU

2111 FORMAT(* Nl.PHC. BL.THL. 8U.THU * I5»5E12.4)

21 FORMAT! I5.E20. 9 /2E20.9/2E20.9 ) CBU=BU*CEXP(-IM*THU) $ CBL=BL*CEXP <-IM*THL

ZF=CEXP(AN6) SZG=CEXP(AN7) $ZH=CONJG< ZF)*KL*CBL

ZI=CONJG(ZG)*CBU*KU *DO 24 I=1»N5 $AI=I~1

ZH=ZH*ZF $ZI=ZI*ZG $ E*2»*(ZH+ZI)

24 PHBU)=E+PHC+AI*AN5 $ICONV = 30 $ICASE= ICASE+1

IF( IGUS>302»300»303

302 BI=U $V16=VA(1) $P16=PHA(D $ GO TO 303

300 BI=A*S3*(SQRTF(1.+S4*S4>-S4>

301 V16=P16=0.

303 CONTINUE

PRINT 2113 .ICOUNT. IGUS.BI, V16.P16

2113 FORMAT(* ICOUNT . IGUS »BI .V16 .P16 */2 15 .3E20.9

CALL PARTA $IF (ICONV) 25.26.26

25 PRINT 205* Nl.ICASE

205 FORMAT ( * NO CONVERGENCE FOR Nl= * 15. * AND CASE NUMBER- * 1 5 ) GO TO 2020

26 ICONV=10 SCALL PART3 $IF( ICONV > 27 .31 .31

27 PRINT 205.N1.ICASE

31 CALL CHARTB

GO TO 2020

18 FORMAT (15)

199 FORMAT!* Zl .Z2.Z3 ( 19 )= * 15 »3C< E12.3 .E12.3 >

END

Figure 2. (Continued)

13
SUBROUTINE SCALE ( 19. 199

>

C PURPOSE IS TO CONVERT THE DC DOWN BY 10-7 AND FORM LOGlO(Yi)

C THE AMPLITUDE OR POWER SPECTRA 5 DECADES ONLY.

TYPE REAL MIN.MAX

COMMON /4/ ICLE(8.2»6> »IDC(8)

COMMON /l/ E(6151) » Yl < 1025

>

»F< 11299

>

GO TO (1.2.3.4.5.6)199

3 IDC(3>=ICLE( 19+1.1. 3>=ICLE( 19+1. 2. 3 > =0 $ RETURN

1 CONTINUE

2 CONTINUE

4 CONTINUE

5 CONTINUE

6 CONTINUE

C POWER CASE AND VOLT CASE

7 MIN=0. $ MAXaO.

IF( 19)8,9

9 IDC(I99»=-7 $ Y1(3)=1.E-7*Y1(3)

8 DO 11 I»1.1024

IF(Y1( I>)12»11»13

12 IF(MIN-Y1(I> >11.11.14

14 MIN = YK I) $ GO TO 11 13 IF(MAX-Y1(I> )15, 11.11

15 MAX = Y1U)

11 CONTINUE

IF(MIN)31.32.32

31 IMIN=AL0G10<-MIN) i MI N=-l.* < 10.

)

**IMIN SICLE ( 19+1 .1 » 199 ) = IMIN-5

32 IF(MAX»33»33.34

34 IMAX=AL0G10<MAX)$MAX=U0.)**IMAX SICLE ( 19+1 .2 » 199 )= IMAX -5

33 DO 21 1=1.1024 $ I F( Yl ( I ) ) 22 »21 .23

22 Yl( I)=-(ALOG10(Y1( I )/MIN)+5. ) $ IF ( Yl < I » ) 21 .21 .25

25 Yl( I)=0. $ GO TO 21

23 Yl( I)=ALOG10(YK I)/MAX)+5. $ IF( Yl ( I )) 24.21.21

24 Yl( I »=0. $ GO TO 21

21 CONTINUE $ RETURN $END

SUBROUTINE PART

DATA( ALPHA=0..1...2».1..01..005)

DIMENSION ALPHA (6>

DIMENSION 21(1025) .Z2 ( 1025

)

»Z3 ( 1025

>

COMMON /5/ BI, SI .52. S3

1.A1

COMMON /16/ V16.P16 .H8

TYPE COMPLEX Zl .Z2 .Z3 .K ,DT .ZF.ZH.ZG.KL .AN2.ATMP.KU.ATN.AN6.AN7 ,ZI

l.CBU.CBL

TYPE REAL KO

COMMON /DAT/ PI.IM.PHS.U

DIMENSION K<1025> »DT(1025) .KL(3) .KU(3>

1 .I2TU025) »VAT(1025)

EQUIVALENCE (K(l) ,DPH(1) ) . ( DT ( 1 ) ,DVA ( 1 )

)

COMMON /l/ PHC,VA(2050> .DPH( 4.2050 )

,

1DVA(4.2050) .H.H1»H2.H3»H4.EP1»EP2.EP3.LBX(5J .LBTR ( 5 ) .EP4 ( 6

)

EQUIVALENCE (DPHC615D »ZK1>) » ( DVA( 6151 ) »Z2( 1 ) ) .(Q4(l) .Z3«D ) COMMON /2/D(2050) »PHA(2050) .PHB<2050> .QK205 0) »Q2( 2050 ) »Q3 < 2050 ) . 1Q4(2050),KEY(6) .NO .Nl »N2 .NH.NP.PO.VO.AO.V » THL .THU.BO »BL »BU,NL »NU

2.NB.KL.KU. AN5.H5, ANU.ANUP .ANUT .A »R1 *R0 .C.WO, T . R.PI2 »AN2 ,AN4»

3AN3.MM.N5.N6.N7.ISKPD.se. OUT ( 6» 11

)

»NR . ICONV.AJ. ICASEt

4P0W1.P0W2.P0WJ.P0W I N.K0.EP5.H7

5.CBU.CBL

EQUIVALENCE! X(l) ,DPH( 1 ) ) . < Y( 1

)

»DPH(2051) ) » < VAT ( 1 ) .DVA( 1 ) ) K I1T(1) .DVA(2051) ) » ( VOUTI 1

>

.DVAI4101 )

)

»<X1(1) »DPH(1025))

TYPE COMPLEX VAT . I IT .VOUT. I2T »IM DIMENSION X ( 1025 ).Y( 1025) » I IT (1025) .VOUT (1025) .XI (1025) .YK1025)

EQU I VALENCE (YK1) ,DPH(4l0l) > . ( 12 T( 1 » »Q3( 1 > ) DIMENSION Tl<4) .T2<4)»T3«4>

Figure 2. (Continued)

14
ENTRY PARTA

C THIS PART USES PURE DIFFERENTIAL EQUATION AND TRAPAZOIDAL SPLINE

C APPROACH TO GET PHA.VA SHAPE IN TIME

DO 1 I-1»N6 $VA< I )=0.

1 PHA(I)=0.

EP77=EP5

EP78=.01*A $ZP=EP78*.3

EP74=EP78

IF(N1>2.3

3 U= A*SINF(PHB<1>) SGOTO 483

2 U«BI

VA<N5)=V16

483 ICONVl=12

DO 460 119=1,10

BET«6AM«XI»0.

PHA(D=P16 $ VA(1)=VA<N5) * VA(N5)=0.

ICONVl=ICONVl-2 $ IF( ICONV1 ) 61 »61 .620

61 ICONV=2 $ GO TO 4 620 ICONV=ICONVl

4 ICONV=ICONV-l $IF( ICONV ) 460 .461 »461

460 CONTINUE $ IF< ABSF ( PHA ( N5 > -PHA ( 1 > )-100.*EP2 >411 ,413 »413

411 IF(ABSF(VA(N5)-VA(1) )-100.*EP2 > 412 .413 »413

412 ICONV=l $ GO TO 405

413 ICONV=-l$ GO TO 405

461 U=U+GAM

VA<1)=VA( 1)+XI

TA=H3*U $TB = PHA(D+PHB(1) $Q3 ( 1 ' =SINF( TB)

401 Q2<l>=COSF(TB) $DPH( 1 » 1 )=DPH( 2 »1 )=DVA( 2.1 ) = 0«

DVA(1.1>=1. $ T3<2)=H3 ST3<1>=0.

DO 6 1 = 1. N2 $13=1 + 1 $TD=H1»VA(I >+TA-trH2*Q3II) $ IV=5

8 TB=PHA( I3)+PHB( 13) SQ2 ( I 3 > =COSF I TB) $ TC=1 .-H5*Q2( I3> IV=IV-1

IF( mio.ii.ii

10 PRINT 13 SCALL EXIT

13 FORMAT <* CONVERGENCE PROBLEMS IN PART A* ) 11 Q3( I3»=SINF(TB)

VA( I3>=TD+H2*Q3( I3» TB=PHA(I)-PHA(I3)+M4*(VA(I3)+VA( I>

DEL=TB/TC

PHA( I3>=PHA( I3)+DEL $ DEL=ABSF(DEL

IF(DEL-EP1)7»8»8

7 TB=PHA( I3>+PHB( 13)

Q3(I3>=SINF(TB>$TC= Q2 ( I > +Q2 ( 13 ) $VA< 13 )=TD+H2*Q3 < 13)

Q2( I3>=COSF<TB>

TB=H5*Q2(I3» $DO 5 11=1.2

Tl( I1»=DVA( I1.I)*(H1+TB>+T3( I1»+H2*DPH( I1,I>*TC

T2(I1)=1.-TB

DVA( I1.I3)=T1(ID/T2<I1)

5 DPH(I1.I3» = DPH( I1,I)+H4*(DVA< 1 1 » 13 >+DVA( II » I )

6 C9NTINUE $GO TO 12

12 B1=PHA(1)-PHA(N5) $B2=VA< 1 > -VA< N5

1450 A11=DPH(1,N5) $A12=DPH( 2.N5

A21=DVA(1»N5)-1. $A22=DVA<2»N5>

DET=A11*A22-A12*A21 $B1«B1/DET $ B2=B2/DET

XI =A22*Bl-A12*B2 $GAM=A11*B2-A21*B1

IF (ABSF( GAM J -EP5» 17,17.19

17 IF(ABSF(XI)-EP74)14,14.19

19 GAM=SIGNF(EP77.PHA( 1)-PHA(N5» ) IF(A8SF(PHA(1)-PHA(N5> ) -.1 ) 199. 199.198

198 XI=SIGNF(ZP»VA(4>*(VA(1)+VA(4)-VA(3 ) -VA(2) »> GO TO 200

199 XI=SIGNF(EP78.(VA(1)+VA(4)-VA(3)-VA(2) )*VA(4)

EP78=.9*EP78

200 EP77=.90*EP77

14 PRINT 1000»U.VA(1).PHA(1).GAM.BET.XI

1 »VA(N5>.PHA(N5»

1000 FORMAT!* U.VA.PHA * /(6E13.4>>

I F ( ABSF ( GAM > -EP2 ) 400 .4.4

400 IF(ABSF«XF)-EP1)40.4»4

Figure 2. (Continued)

15
40 RETURN

405 PRINT 999»<PHA< 17) ,VA( 17) ,I7=1,N5>

999 FORMAT(8E12«4)

ENTRY PARTB

C THIS PART GETS THE INTEGRAL EQUATION PART USING QUASILINEAR

C ANALYSIS AND THE TRAPAZOIDAL METHOD FOR THE ERROR ESTIMATE

DPH(1»1)»DPH<2»1)=DPH<3,1>=DVA(2»1 >=DVA« 1*1 )=0. $DVA<3*1)=1.

60 J6*l SGO TO 44

30 GO TO (21*28,28, 28, 28, 24*24, 60»35>IWAY

21 A12=-H4$A11=A22 = DET1= 1. $ T2 ( 2 ) =T2 ( 3 ) =0. $T3(2)=H3 $T3(3)»0.

DVA(2»1>=DPH<2»1>=0.

DO 22 1=1. N2 $12=1+1 $A21=-H2*Q2( 12) $DET=DET1-A12*A21

T3(1>=-Q3(I2)

T2(l>=-Q4< 12) $DO 23 16 = 1*3

B1=(DPH( 16* I)+H4*DVA(I6*D+T2( 16 > )/DET

B2=(H1*DVA( 1 6. 1 )+T3( I6>+H2*Q2( I )*DPH( 16*1 ) J/DET

DPH(I6»I2>= B1*A22 -B2*A12

DVA(I6»I2>= B2*A11 -B1*A21

23 CONTINUE

22 CONTINUE

Bl= -DPH(1»N5> $B2= -DVA(1»N5)

All=DPH(3»N5> $A12=DPH(2»N5>

A21=DVA(3*N5)-1. $A22=DVA( 2 »N5

>

DET=A11*A22-A12*A21 $B1=B1/DET $ B2=B2/DET

XI =A22*B1-A12*B2 $GAM = AH*B2-A21*B1

199=2 $ GO TO 42

24 TEMP=EP4( 1)SI99=1$IC0NV=IC0NV-1 $J6=1$IWAY=8 SIF ( ICONV) 39,26,26

26 DO 250 L6=2»6 $ IF ( EP4( L6)-TEMP ) 25 »250,250

25 TEMP=EP4(L6> $J6=L6

250 CONTINUE $ IF (J6-D 62.62*28

28 DO 32 I=1»N2

DPH(2»I»= ALPHA (J6>* DPH ( 1 * I >+PHA ( I

>

32 DVA(2»D= ALPHA( J6>* DVAd.I )+VA(I) $TC=U+ALPHA( J6) *GAM

GO TO <46»49> 199

29 EP4(J6)=0. $D0 34 I=1»N2

SMALL=Q4( I)**2+Q3( I )**2 $ IF(SMALL-EPl )341*341.344

344 PRINT 345»SMALL»U*I,J6

345 FORMAT(* ERROR AT POINT * »2E12. 3*215)

341 EP4< J6)=EP4( J6)+SMALL

34 CONTINUE

TEST=EP4(J6> $ IWAY=J6 $ J6=J6+1

IF( TEST -EP3)51»30.30

51 IWAY=9 SGO TO 46

62 PRINT 63»EP4»EP3 SRETURN

39 PRINT 41,EP4,EP3

41 FORMAT( * 10 LOOPS IN QUASI NO CONVERGENCE */(6E13.4))

35 RETURN

42 DO 43 I=1*N2

DPH(1»I)=DPH(1,I)+XI *DPH(3»I )+GAM*DPH(2» I

>

43 DVA(1»D=DVA(1,I>+XI *DVA ( 3 » I )+GAM*DVA( 2 » I > IWAY=J6=2 $ GO TO 30

44 DO 45 I=1*N2

DPH(2»I , =PHA( I) $DVA<2»D=VA( I) 45 CONTINUE $ TC=U SGO To 49

46 DO 47 I=1»N2 $PHA< I )=DPH( 2 » I ) $VA( I ) =DVA ( 2 I

)

47 CONTINUE $U=TC $GO TO 30

49 DO 52 1=1. N2 $TB=DPH( 2 » I> +PHB ( I ) SQ2 ( I > =COSF( TB)

Ql( D=-A*SINF(TB)+TC -ANUT*DVA( 2 . I

)

IF( ISKPD)54»52

54 DO 53 I3=1»N2

L=I-I3

IF(L)55.56.56

55 L=N2+L

56 CONTINUE

53 Ql( D=Q1( I)-D(L + 1)*DVA(2»I3)

52 CONTINUE

DVA(2»N5)=DVA(2.1)

Figure 2. (Continued)

16
DVA(2»N6)=DVA<2»2)

DPH(2*N5)»DPH(2»1J

DPH(2.N6)=DPH(2»2)

63 FORMAT (* NO CONVERGENCE POSSIBLE * /(6E13.4))

Q2<N5)=Q2<D SQ1(N5)=QKD $Q1 (N6 ) =Q1 (2

DO 50 I=1»N2 $ 12=1+1 $13=12+1

04

I2)=DPH(2.I2>-DPH(2»I>-H4*(DVA<2»I2>+DVA(2»I> ) 50 03( I2)=H7*(DVA(2»I2)-OVA(2»n-H4*(QK I2>+Q1( I) ) ) Q3<1> =Q3<N5)$Q4(1> =Q4(N5» SGO TO 29 SEND

SUBROUTINE CHART

DIMENSION Zl ( 1025 ) »Z2 « 1025 )»Z3( 1025)

COMMON /5/ BI.S1.S2.S3

1.A1

COMMON /4/ ICLE(8»2»6) »IDC(8)

COMMON /DDC/LU»LUC*IFL

COMMON /DD/ IN.IOR»lT,IS»IC»ICC»IX.IY

COMMON /TAB/ VI ( 60 ) , V2 < 121

»V3 ( 19

TYPE COMPLEX Zl .Z2 »Z3 .K »DT .ZF .ZH.ZG.KL »AN2 .ATMP.KU.ATN.AN6. AN7.ZI

l.CBU.CBL

TYPE REAL KO

COMMON /DAT/ PI.IM.PHS.U

DIMENSION K < 1025 ).DT (1025) »KL(3)»KU(3)

1 »I2T(1025> .VAT(1025)

EQUIVALENCE <K( 1 ) »DPH( 1 > > » (DT < 1 ) ,DVA( 1 >

COMMON /l/ PHC,VA(2050» »DPH( 4*2050 )

1DVA(4»2050)»H,H1,H2,H3*H4,EP1,EP2.EP3*LBX(5),LBTR(5),EP4(6)

EQUIVALENCE (DPH(6151) *Z1(D ) » ( DVA< 6151 ) »Z2 ( 1 ) ) , < Q4( 1 ) ,Z3( 1 >

COMMON /2/D( 2050 >»PHA( 2050' »PHB(205Q) »Q1 < 2050 > ,Q2( 2050 ) ,Q3(2o50)*

104(2050) »KEY(6) ,N0 *N1 *N2 »NH»NP ,P0»V0»A0,V, THL,THU»BO,BL*BU*NL,NU

2»NB»KL»KU. AN5.H5, ANU,ANUP»ANUT »A»R1 »R0 ,C» W0,T , R.PI2 ,AN2 ,AN4,

3AN3»MM»N5»N6»N7»ISKPD,SC. OUT ( 6» 11

»NR» ICONV. AJ, ICASE.

4P0W 1 » POW2 » POWJ » POW I N » KO , EP5 . H7

5.CBU.CBL

EQU I VALENCE ( X ( 1 > »DPH< 1 ) )»(Y(D .DPH(2051> ) »(VAT(D ,DVA(1») . 1( I1T(1) »DVA(205D ) »(VOUTU> »DVA(4101) ) » (XI ( 1 ) ,DPH( 1025 ) ) TYPE COMPLEX VAT* I IT »VOUT* I2T »IM DIMENSION X(1025) ,Y<1025) » I IT (1025) .VOUT(1025) ,X 1(1025) »Y1(1025)

EQU I VALENCE (YK1»,DPH( 4101) ) » ( I2T ( 1 ) »Q3( 1 ) ) ENTRY CHARTA

REWIND 1 WRITE TAPE 1*Z1,Z2»Z3

REWIND 1 IN=IOR=IT=IC=0 $IS=1 $LU=LUC=0 $IFL=30

CALL DDINIT(3»20H EGJOHNSON.EXT 3706 ) S1»S1/(R1*ANUT>

PRINT 100»A1»S1 $ PRINT 110»A»ANUT

110 FORMAT (* A»ANUT *,2'E12.4)

ENCODE ( 40* 100* VI )A1»S1

100 F0RMAT(16H Al »SQR( A) /NUT • / »2E12«4>

IX=89 $ IY=500 $ CALL DDTAB $ CALL DDTEXT15.V1)

ENCODE(480.1»Vl)NQ,N2»NB»NL»NU.ANU»ANUT,A.V0.A0,PQ.WQ»RO»C»R,H,AJ»

IT

1 FORMAT (89H THE JOSEPHSON SPECTRUM CONTROL PARAMETERS THAT REMAIN

1C0NSTANT '/FOR THE ENTIRE RUN.*/ .40H NO, DATA POINTS.BEAT ,LOWER»U

2PPER FREQ.'/ »5I5» 48H'/ NU.NU TOTAL.SCALED A.VOLTS. AMPERES* WATTS

3'/ » 6E12.4* 62H '/RADIAN FREQ, OHMS* CAPACITANCE. RESISTANCE* ST

4EP INTERVAL'/ . 5E12.4.36H '/SUPER CURRENT AMP. .PERIOD SEC.'/

5 2E12.4.24H '/CASES TO FOLLOW.'. > ENC0DE(32»99»LBX(1) ) 99 FORMAT (32H JOSEPHSON CURRENT IN TIME CASE ) IX=89$IY=950 JCALL DDTAB $CALL DDTEXT ( 60.V1 ) SCALL DDFR SRETURN

ENTRY CHARTB

READ TAPE 1»Z1.Z2»Z3

REWIND 1 C Q2(D IS A*SINF(PH>

Figure 2. (Continued)

17
DO 2 I=1.N2

2 Q2t I»«A*SINF(PHA( I)+PHB( I)

)

X1(D=0. $DO 3 1=2.256

3 Xl( D»XK I-D+4. $IFL=60

ENCODE* 8 .98. LBX (5) ) I CASE

98 FORMAT* 15)

NS = N2/257+l 500 4 I8=1»NS $I9 = 256*( 1 8-1 > $00 5 1 = 1*256 $11=1 + 19

5 Y1(I)=Q2(U» $CALL FILMGRAF ( XI » Yl *256.LBX. 1*0 ) 4 CONTINUE

DO 6 1=1. N2

6 Q2< I>=Q2< D+QK I>+ANU*VA«I>

DO 7 I=1»NH $X(I)=Q2« I

)

7 CONTINUE $D0 8 I=NP»N2

8 Y( I-NH)=Q2( I) $CALL REORDER ( X.Y.MM ) $CALL CFFTRC(X. Y.MM.SC. 1

)

CALL REALTRAN(X.Y.MM.l.l) $D0 901 = 1. NP $ I IT ( I

) =CMPLX(X< I

)

»Y< I )

)

90 CONTINUE

C CALL REORDER(A.B.M)

C REVERSIBLE PERMUTATION OF REAL SEQUENCE

C FROM FIRST-LAST-NORMAL SEQUENCE

C TO ODD-EVEN-REVERSE BINARY SEQUENCE.

C SEQUENCE LENGTH IS N = 2**M

C CALL REALTRAN(A.B.M.NE.INV)

C IF(INV.GT.O) UNSCRAMBLE THE TRANSFORM OF A REAL SEQU'-NCE.

C IF(INV.LT.O) SCRAMBLE THE TRANSFORM OF A REAL SEQUENCE.

C INPUT AND OUTPUT ARE IN NORMAL SEQUENCE,

C SEQUENCE LENGTH IS N = 2**M

C NE MUST AGREE WITH SIGN OF EXPONENT IN TRANSFORM DEFINITION.

C INNER LOOP SINES AND COSINES COMPUTED

C RECURSIVELY BY SINGLETONS 2ND DIFFERENCE ALGORITHM,

C INITIALIZED FROM A DATA TABLE.

C DISCRETE COMPLEX FAST FOURIER TRANSFORM.

C CALL CFFTRC(A.B.M.SC.NX)

C INPUT A(J) + I*B(J) IN REVERSE BINARY SEQUENCE.

C OUTPUT A<K> + I*B(K) IN NORMAL SEQUENCE.

C SEQUENCE LENGTH IS N = 2**M

C SC IS REAL SCALING MULTIPLIER.

C NX IS THE SIGN OF THE EXPONENT IN THE TRANSFORM DEFINITION.

C INNER LOOP SINES AND COSINES COMPUTED

C RECURSIVELY BY SINGLETONS 2ND DIFFERENCE ALGORITHM.

C INITIALIZED FROM A DATA TABLE.

I1T(1>= UT(1)+N1*ANU $IlT(NU+l) = IlT(NU+l ,

( KU( 2 ) *ANU+K.U< 3 > ) *CBU

1 1 T ( NL+1 > = I IT ( NL+1 ) + < KL ( 2 > *ANU+KL ( 3 > > *CBL

DO 10 1=1. NH $X(I)=VA( I

)

10 CONTINUE $D0 11 I=NP.N2

11 Y(I-NH)=VA( I) $CALL REORDER ( X.Y.MM ) $CALL CrFTRC(X. Y.MM.SC . 1

)

CALL REALTRAN(X.Y.MM.l.l) $DO 12 1=1. NP $VAT ( I )=CMPLX(X ( I ) . Y( I )

)

12 CONTINUE

VAT(1>=VAT(1)+N1 $VAT(NU+1>=VAT(NU+1>+KU(2> *CBU

VAT(NL+1)=VAT(NL+1)+KL(2) *CBL

VJUN=V0*VAT(1) $ AIJUN=A0*I1T(D $ ENCODE* 64 .9111, VI >VJUN,AUUN

9111 FORMAT ( 30H VJUN ( VOLTS ) . IJUN< AMPERES M / » 2El7«5>

DO 13 1=1. NP

VOUT(I)= Z2< I )*I1T( I )+VAT( I) $I2T ( I ) =VOUT( I ) /Z3 < I

)

13 CONTINUE $D0 14 1=1.1024 $X1(I)=I-1

14 Y1(I)=0.

C DROP NP POINT IN THESE PLOTS. SEE PRINT OUT FOR IT. POW1=POW2=POWJ=0.

NS = NH/129 + 1 $D0 15 18=1, NS $ I 9=128* ( 18-1 ) $K7=I8-1

DO 21 1=1,128 $ I2=(I-l>*7+3 $11=1 +19

21 Y1(I2)=CABS(I2T(I1) ) $ENC0DE(40,97,LBTR( 1 ) > ICASE, 18 97 FORMAT(26H AMPLITUDE I2TRAN. ICASE= 215. 4HPART ) CALL SCALE(K7.1)

CALL FILMGRAF(X1.Y1.1024»LBTR.1.0»

C PLOT OF CURRENT .12. AMPLITUDE THROUGH DETECTOR Z3. AS A FUNCTION

C FREQUENCY

DO 22 1=1.128 $I2=( I-l>*7+3 $11=1+19

22 Yl( I2>=CABS«V0UT( ID ) $ENCODE( 40.96.LBTR* 1 )

> ICASE. 18 96 FORMAT (26H AMPLITUDE VOuTTR. ICASE= 215. 4HPART )

Figure 2. (Continued)

18
CALL SCALE! K7. 2)

CALL FILMGRAF(X1»Y1»1024»LBTR.1»0)

C PLOT OF OUTPUT VOLTA6E ACROSS DETECTOR

DO 23 1=1.128 $I2=(I-l)*7+3 $11=1+19

23 Y1(I2>= CANG( I2T( ID )-CANG(VOUT( II)

ENCODE(40»85.LBTR(l) ) I CASE. 18 85 FORMAT (26H PHASE I2-VOUT TRAN ICASE= 215. 4HPART ) CALL SCALE(K7»3>

CALL FILMGRAF(X1»Y1»1024.LBTR.1.0)

C PHASE OF I2-VOUT

DO 24 1 = 1.128 $I2=(I-D*7+3 $11=1 + 19

24 Y1(I2)= CABS(VAT( II) ) $ENCODE (40,95 ,LBTR( 1 >) ICASE. 18 95 FORMAT (26H AMPLITUDE VPHTR. ICASE= 2l5» 4HPART ) CALL SCALE(K7.4)

CALL FILMGRAF(X1.Y1»1024»LBTR.1»0'

C PLOT OF AMPLITUDE OF VOLTAGE ACROSS JUNCTION. PH( DOT)

DO 25 1=1.128 $I2=(I-l>*7+3 $11=1+19

TEP= 2.*REAL< I2T( 1 1 )*CONJG( VOUT ( II > ) ) $POW2=POW2+TEP

25 Y1(I2)=TEP $ IF(I9)26»27

27 TEP= Y1(3) = .5*YK3» $POW2=POW2~TEP

26 ENCODE(40»94.LBTR(1> )ICASE»I8

94 FORMAT (26H POWER I2V0UT SSTR. ICASE- 2I5» 4HPART ) CALL SCALEIK7.5)

CALL FILMGRAF(X1»Y1»1024.LBTR.1.0)

C DETECTOR POWER FOR EACH FREQUENCY. P0W2 IS TOTAL POWER.

DO 28 1=1.128 $I2=( 1-11*7+3 $1 1=1+19

TEP=2.*REAL( I IT ( II )*CONJG ( VOUT ( II ) )) $POW1=POW1+TEP

28 Y1(I2>=TEP $IF(I9)29.30

30 TEP = Y1(3»=.5*YK3) $P0W1=P0W1-TEP

29 ENC0DE(40»93.LBTR<1> )ICASE»I8

93 FORMAT (26H POWER I1VOUT SSTR. ICASE= 215. 4HPART ) CALL SCALEIK7.6)

CALL FILMGRAF(X1»Y1»1024»LBTR.1.0>

C POWER IN 22+JUNCTION POW1 IS TOTAL POWER IN POWJ IS TOTAL IN JUN

C KO IS THE EFFECTIVE IMPEDENCE. NOTE ALL PLOTS ARE LOG10 EXCEPT

C PHASE CASE. THE POSITIVE .UP CASE IS SCALED INDEPENDENT OF

C THE NEGATIVE CASE, DOWN CASE. THIS IS DONE IS SCALE SUBROUTINE.

C FILMGRAF THEN DOES A SCALE CHANGE TO MAKE THE DATA FULL SCALE

15 CONTINUE

OUT ( 1 1 ) = ( U+N1*ANUT ) *K0 $OUT ( 2 » 1 ) =OUT (5.1) =OUT ( 6. 1 ) =0.

0UT(3»D=BL $0UT(4.D=BU $OUT ( 3»2 ) =THL $OUT( 1 »2 ) =PHC+PHA( 1

0UT(4»2'=THU $OUT(2»2)=OUT(5»2)=OUT(6»2)=0.

0UT(1.3)=0UT(1,D*( UT(1) + I2T(D) $OUT( 2 .3 ) =OUT ( 5 .3 ) =OUT( 6.3 ) =0.

OUT (3.3) =BL*REAL ( ( 1 1 T ( NL+1 > + 1 2 T ( NL+1 >

) *CEXP ( I M*THL )

0UT(4.3>=BU*REAL( ( I 1T( NU+1 >+I2T (NU+1 ) ) *CEXP( IM*THU)

POW I N =OUT ( 1 » 3 ) +OUT ( 3 » 3 > +OUT (4.3)

P0W1=P0W1+REAL( I1T(NP)*C0NJG(V0UT(NP) ) ) P0W2=POW2+REAL( I2T ( NP) *CONJG (VOUT (NP) )

POWJ = REAL( I1T(1)*C0NJG(VAT(1) ) +1 IT (NP)*CONJG( VAT(NP) )

DO 32 1=2. NH

32 P0WJ=P0WJ+2»*REAL( I IT ( I

) *CONJG( VAT ( I ) >

DO 34 14=1.6 $ I5=KEY( I4>+1

OUT( 14,4) = 2.*REAL(I1T( 15 )*CONJG( VOUT( 15) )

OUT( 14, 5)= 2.*REAL( I2T( 15 )*CONJG( VOUT( 15 > > > 0UT(I4.6)= REAL(Z1(I5)> $OUT ( 14.7 )=AIMAG(Z1 ( 15 >

0UT(I4.8)= REAL(Z2(I5)) $OUT ( 14,9 >=AIMAG(Z2( 15 )

0UT(I4,10)=REAL(Z3( 15) ) $OUT ( 14, 11 ) =AIMAG(Z3( 15 )

34 CONTINUE

DO 35 14=1,6,4 $ DO 36 17=4,5

36 OUT( I4»I7)=.5*0UT(I4.I7)

35 CONTINUE

ENCODE (968 ,9,V2 ) POWJ.P0W1 ,P0W2 .POWIN.KO. ( KEY ( I ) » (OUT( I.J) »J=1»11)

1.1=1.6)

9 F0RMAT(28H POWJ .P0W1 .POW2.P0WIN .KO, • / 5E12.3,2H«/ . 138HKEY.B.PHC.P0WIN.P0W1.P0W2. , /Z1,Z2»Z3»/ . 2 6(I4.5E12.3.2H«/,6E12.3»2H«/ )) ENCODE (152. 91 »V3>EP1.EP2»EP3.EP4

Figure 2. (Continued)

19
91 FORMATUOH ERROR LIMITS EPl»EP2»EP3 •/ EP4(6) •/ »3E12.3»2H«/

16E12.3. 2H«. ) IX=89 $ IY=950 $ CALL DDTAB $ CALL DDTEXK 121»V2

)

IX=89 $ IY=600 $ CALL DOTAB $ CALL DDTEXT(8tVl»

IX=89 $IY=500 $CALL DOTAB SCALL DDTEXT ( 19»V3

)

IX=89 $IY=400 $ ENCODE(480»911.V2> (IDC(I) *I=1*6) »( H ICLE(I.JtL) »I=

1 1.8) »J=1»2> »L=1.6)

CALL DDTAB $ CALL DDTEXT (60. V2 ) $ CALL DDFR

911 FORMAT (58H SCALE EXP ONLY DC SET THEN MIN.MAX SET FOR EACH GRAPH

1'/ » 6I4*2H«/ »6( 16I4»2H»/)

)

RETURN

END

Figure 2. (Continued)

20
r

THE MAIN PROGRAM JOSPHC

i f READ ) NO.

5( READ ~*) Z1 ,

~i

(MAY BE SKIPPED)

(REPEAT AS

NECESSARY)

(ICOUNT CONTROLS

PATHS)

(IGUS CONTROLS

INITIAL DATA

GUESSES)

L.

( IF(I8) ) 4.5.6

( READ ~) UL,

(IF(ISKPD)) 15,19

15 ( READ ) IRED

19(CONTINUE )

( CALL ) CHARTA

2022 ( READ ) ,N1

2020 ( READ ) ICOUNT . IGUS

CALL PARTA

PARTB

CHARTB

( GO TO ^) 2Q2Q

._l

THE MICROFILM PLOT

ROUTINES

L

DDINIT

FILMGRAF

DDTAB

DDTEXT

DDFR

SCALE J

fTTHE FAST FOURIER

TRANSFORM ROUTINES

REORDER

CFFTRC

REALTRAN

L.

r

THE CHART SUBROUTINE CALLS

THE MICROFILM

AND THE FAST

FOURIER TRANSFORM

TO GENERATE THE

SPECTRAL DATA

ON MICROFILM.

ENTRY CHARTA

ENTRY CHARTB

I l_

THE PART SUBROUTINE GENERATES

THE INITIAL

GUESS FOR PH AND THE QUASILINEAR PROCEDURE

FOR THE SELFCONSISTANT PH.

ENTRY PARTA

ENTRY PARTB

Figure 3. The flow chart for Josephson steady state program.

21

l>- iiyflftn 1
o o o o

4> «3 o

cd o O o 13 O • • O a a 4- 1 1 OOO

ao tO

O

a o • a o o •h fO CT o 4 1 o CD o o a a o a o a o o o 1 1 1 o o o ao o ao o • • o CVJ CM CD a •H v4 4- Lrt aO t4 o o (VI eo to a a CD • • • • CD CD CM ^ to

w*

Ooo O

CD

O•oCD

CD 1 1 1

CO 4- 1 CD CD CD CD o o O CD O CD CD o o CD o CD CD CD CD 1 o O o O 4- 4- 4- o • • o a CD CD O a o a 4- CD a O O a 1 1 CD IN CD O O o CD V0 CD CD CD • O • • • • H aoo oao o• CD 1 CD O CD

*H 4 1 O a a CD fO CT a o Q in e CD a

4- o O o CD CD CD CD o • • CD CD 1 1 o o o 4» CD N- O •4- a t4 4- i>»» S0 t4 o CD CD iH 4- CD • a CD vO ^ h1-4 o o J" CD • • • o o • CD i-i CM CT» o o CM CD UN 4 CD

«H CD

a CD

O o CD ina in CD CD

O 1 o o o en CD O CD N. • • CD • CD CD CD h- o a 4 4- 4 4> CM 1 1 Id CD CD CD

C\J

a 1-1 CD

a

O • 1 4- CM CM CM

CO iX\ in •H 10 vO \D

CD o • • • • a £9 ro i-i fO ro to o 4 4- in 1 o t 1 1 o o O CD o o CD 4- o O CD CD a CD CD o o o • • CD CD o O O o o o iH »4 4 CD O CD O CD a o

U\ vO CD 4 4 4 4- CM i-4 a CD CM ^k. O 00 CD • o a. CD CD cr> j- o CM CD vO «-i t4 CM X « CD in 4- CM CM CM CM

a. I— vD CD a» u> r» \0 r» vO UJ UJ •» i-t vO • • • • • • * 4- (/) 3 > • 4" fO j- to 4" fO •* o < 00 •k CM 1 1 1 a. o o •» M 1 UJ 1 II II •~» -1 CX) « CD M M « X •k i-« CD vl o *-l o

QC CD CP o>

a •— 00 CD CD CD CD CD CD

•k CD t~t

H Q. •k z> CD CD CD CD CD CD

O o ^^

w

m -1 CO iH 4 4> 4- 4- 4- 4-

•> • to to CO CO *-* 4 CD r- < (M fw <a«o

** ~* rsiM

M •k •k

X CD fvj

X 00 <M

X CO CM

«s •k •• •> o H- a. CD yt

a ro fNJ 0L tO CM t eg t\) 3 X z •» vO %<> • 4 N. • 4 N.

-> »-•

M INI 2 a r> i-« «x • • «a • • < • •

•» T-t •k » •> •» o > t\J -a- ;» CM .*• > CM 4 o 1-1 *-l -i •h o •» 1 • 1 » 1 ou

enso

M-lO

inM•rH

IM

M t-t => 22
signals are to be applied, whether we want to compute the integral equation part, (ISKP is zero means that the frequency spectra of impedance is to be ignored. ISKP equal 1 means that it is considered and the integral equation is now used.); and finally what the initial case number label for this computation run is to be, ICASD(5). The remaining cases are incremented by one for each new case. The Nl, PHC, etc. , row is repeated for each new case. Here Nl (4) is the chosen rotation

frequency for the Josephson oscillator, PHC (-3. 1416-0) is the chosen initial phase assuming PHA

is zero initial value. BL is the corrected value for the lower frequency r. f. applied voltage as defined on page 1. This BL (2. 4+5) is changed from the BL (7. 68+2)put into the program and used in the transient program in the appendix. The lower r. f. phase is THL(0. ) in this example. The BU(0. ) is likewise scaled from BU(0. ) and finally we have the upper initial phase THU(0. ). The next data row has two control variables ICOUNT and IGUS as well as input data that is used

as initial guesses for U, V (1), and cp (1). The ICOUNT(l) has two options--if it is zero, it causes

3. 3. the program to request a new set of Nl, PHC, BL, THL, BU, THU; and if it is equal to one it causes

the program to proceed to compute the case just read in. The IGUS has three options--if it is -1, it causes the program to use the computed previous values of U, VA(1), PHA(l) as the initial guesses

for the new set of data; (this option is useful if you are simply making small changes in previous data

cases), if IGUS is zero, it uses some computed guesses based on the approximation that the r. f. power

is not importantly affecting the initial data, hence the d. c. case can be a good basis for initial data

guesses; and finally, if IGUS is one, it uses the currently read in values BI for U, Vl6 for VA(1), and

Pl6 for PHA(l). These values are gotten from the program in the appendix or from other calculations.

The final rows of information in the printed output depends on how the calculation proceeds. If the calculation is normal, then you get what is shown in the sample case. If the calculation does not work normally, then you will get error messages and some detailed printouts. Please look at the

Fortran listings to understand these messages and printouts. The U, VA, PHA sequence is repeated

until there is convergence in the PARTA subroutine. The data printed out is the current U(-2.6+l),

VA(1) (4.96), PHA(l) (-3.62, 6 U (1.6177), change in cp (1), the change in VA(1), the VA(N5) (4.6127),

and finally PHA(N5) (-3. 6492). Please note that convergence occurs when the VA(1) equals VA(N5)

and PHA(l) equals PHA(N5). Also note that the phase PHA(l) is unchanged for the entire calculation

sequence. Here (1, N5) is equivalent to (0,N ) in the analysis.

5. 2 The Microfilm Output

We look at selected plots of the microfilm output of a single computer run. This data corres- ponds to the same run as the printed output. Figure 5 has information that remains constant throughout the different runs. The only information that may not be self evident is the SQR(A)/NUT term (. 5) which is the characteristic frequency ratio of the circuit. It isWa/vrr, and it measures the capacitance

level; the smaller it is the less important is the capacitance term in the circuit.

Figure 6 has 8 plots of 256 points each of a single period of computation. Please note that the

FILMGRAF routine will independently scale the individual plots so that the variations are full scale

in the y axis. If the range of y variables are both negative and positive then the scale is plus and

minus one. If the y variables are only one sign then the range is zero to minus one or one to zero

23
T*J J0t«»»t0>. l?IC'»u» C0«'*Ol »«»»»t'l»» '«»' »t«"l» COM'«V

'a* tmc txUM »u*i NI.O'T* »0!Vtt.M«'.l.O"t«. J»»l» »•!•

II till 1 • t HU.HU TOT«i. »C«l.te «.»Ol.'».«"»CH». •«"»

I ttM'ICI J »««»<C0! I IIIMII I HII'tM I )IM )T 4 Mft-IH

*tC!tN rK«.ew"S. C'MCPtHCt. *CI1S'"iCI »'f» !S'|t,«k s.?stt>tn i it»».is» • *i"-iti t tan. m i iMi-nf

iu»t» eumv •"» ,k»ioo sec MIIO-IM I (MJ-li:

:>$[S *: 'o»»o»

ii.miiinii

I MII'IM • IM1-II1

Figure 5. Sample microfilm output - the control parameters.

24
depending on the sign. We emphasize that it is possible for the different plots to have very different

scales for their respective data. In this sample case the figure 6 set has exactly the same scale for

all 8 plots.

Figures 7 and 13 are the first two plots of a series of 8 for the amplitude at each frequency of the current through the detector. The first figure gives the first 128 harmonics, the latter figure gives

the next 128 harmonics. The remaining 768 harmonics are not shown. The d. c. line amplitude is -7 reduced by 10 in the first plot. Then, in this set of 128 harmonics as well as each succeeding set

of 128 harmonics, the largest positive value is found, the l°g in is taken and the integer part is extracted.

Then all positive nonzero y values are scaled relative to this integer. Only the top five decades of log amplitude values are computed and put into the y axis. The remaining y values that are less than

this allowed 5 dacades are set equal to zero and is thus plotted that way. For those plots which have

negative nonzero y values the largest negative value is extracted, the integer part is removed, these y

values are scaled down by that number, and the log is taken and those numbers that are 5 decades

less than this negative maximum are set to zero. The numbers are then made negative. Both the negative and positive values that have been processed as above are now subject to further scaling in

the subroutine FILMGRAF. That scaling is the same as already discussed for the time series plots.

The above scaling takes place for all the figures 7 through 18 as well as those plots that have

not been shown. Figure 19 has at its bottom a set of scaling information which in conjunction with the

information given at the bottom of each figure permits the relative determination of the amplitudes for

each spectral line that has been plotted. The row of -7's with a zero mean that only the phase differ-

-7 ence series of plots does not scale the d. c. frequency line; the rest are scaled by this 10 . This was

done to allow detailed plot information of the higher frequencies without the d. c. values causing scaling

problems. After the -7 row are six rows of numbers. The first eight numbers of each row defines

the scaling used on the negative numbers in the spectrum. Note that only the last row which corres- ponds to the power spectra across the junction has negative values. The last eight values correspond

to the scaling of positive numbers. The rows are related to the plots as follows.

Row Plot

1 Amplitude of the current across detector (AMPLITUDE I2TRAN)

2 Amplitude of the Voltage across detector (AMPLITUDE VOUTTR)

3 Phase difference between detector current and voltage

(PHASE I2-VOUT TRAN)

4 Voltage across R (AMPLITUDE VPHTR)

5 Power across detector (POWER I2VOUT SSTR)

6 Power across Josephson junction (POWER I1VOUR SSTR)

We illustrate the reading of the plots to get an absolute number for the power across the

junction at the 4, 8 harmonics and at the 132, 133 harmonics respectively. The first two are gotten

from figure 12 and the latter from figure 18. The conversion to watts is given in figure 5 as (4. 8830 -11).

The fourth harmonic has value full scale (5. 918). The eighth harmonic has value (4. 7 X (4. 4410)/ 6. 2

= 3. 34).

25
JOSEPHSON CURRENT IN TIME CASE

Sc«:#

X - Mi«

o coo.oco 0.000*0001.020*003

1 020.005

Mm

•2.560*002

2.5(0*002

M«. 2.560*002

2.5(0*002

*-» W2 M«, 03 U24) «| J?:50:52 U.l.2" J. 256 P«'»t». Ff«»#

Figure 6. The next 8 plots give one complete period, each plot 256

points in time.

26
JOUPWMN CUMflfT \H HI* CASC

Of0*190 000*000! OtO*0H

1 Otf«Nf

2IH*0ff

*tM*0Ot

»-• !tt» *•, 01 -UK: *>i I? 10 *• t!.il« 1.

?1M*0tt

2fJM*00t

*Vm*>. '»«

Figure 6. (Continued)

27
JOttPMSON CUMCNT IN Tl?€ CASC

8 000*000

J Ml. MO

t wooes -2 K0O02

•2 940*002

2 140O92

2 HO.002 »* !f?2 «•» OS ;ltli «t !?:M-M .l.iK ». 3 2K »••»!•. r^flM

Figure 6. (Continued)

28
JOSEPHSON CURRENT IN TIME CASE

N»( t iflt sc«:#

Hiit 0.000*000

0.000*000

Ha. 1.020*003

1.020*003

• Hi*

-2.560*002

-2.560*002

2.560*002

2.560*002

B«i> !5"*2 H« y 03 1124) at 17:50:53 Ll.L2« I. 256 ^a<»t». "'•••

Figure 6. (Continued)

29
JOSEPHSON CURRENT IN TIME CASE

0.000*000

0.000*000

1.020*005

1

.

020*005

•2.599*002

•2.940*002

M.,

2.9(0*002

2.940*002

»«• W2 «•» 0! (124) at 17:90:55 Ll.l2« I. 294 ••••i», f'—

Figure 6. (Continued)

30
JOSEPHSON CURRENT IN TIME CASE

Sc.U

X • Hi*

0.000*000

0.000*000

1. 020*003

J. 020*005

• Hit

•2.553*002

•2.5(0*002

H«. 2.5(0*002

2.5(0*002

Rum 1972 M«, 05 (1241 at 17:50:55 LI.L2* J. 25( *»'M». TfMt

Figure 6. (Continued)

31
JOSEPHSON CURRENT IN TIME CASE

Sec It

X • Mi*

0. 000*000

0.000*000

1.020*009

1.020*00!

I • Hi*

•2.500*002

•2. MO'* 002

M«. 2.900*002

2.500*002

«U 1972 M«» 05 (1241 «t 17:50:55 U.l2» 1. 250 P>0iM». rr«M 5

Figure 6. (Continued)

32
JOSEPHSON CURRENT IN TIME CASE 6

M • Mi* M«i i • Mi* M« a F„»cti»*- 0.000*000 1.020*00! -2 559*002 2.5(0*302

Scala 0.000*000 1.020*003 -2.567*002 2560*002

Pur. !9?2 May 03 (1241 at !7;50:53 Lt.l2» I. 256 ^fiii, '•••• 2

Figure 6. (Continued)

33
ii.

AMPLITUDE I2TRAN, ICASE= 1PART

X Mi« M.i » • MiK N«!

f u»ct't» 0. 000*000 1. 025*005 0.000*000 4 540*000

$€•!• 0.000*000 1.025*003 0.000*000 4.540*000

»u» 19?2 M«y 03 1124) «t 17:50:5? L I . l2» t. 1024 ••mi. f' —

Figure 7. The first frequency block for the LOGjq amplitude of the

current through the detector. (Notice only £ = 4, 8, etc., are

non zero. ) t is the "frequency" relative to UU . 34
AMPLITUDE VOUTTR, ICASE=

J_L

tPART

X - Mi» M«. » • Mi» H«.

f u»c 1

#k 0. 000*000 1 . 029*009 0.000*000 5.059*000

Sc.la 0.000*000 1.025*00 J 0.000*000 5.055*000

*«•. 13?2 M«y 0J 1124) at t?;50:5? Lt.l2« I. 1024 P#i»t». * f— 10

Figure 8. The first frequency block for the LOGjq amplitude of the

voltage across the detector.

35
nr

PHASE I2-V0UT TRAN 1CASE= J PART

X • Mi* M«i l • Hit "•• r«»ctif» o.ooo*ooo 1 . 023*009 •5.M1-01! 2.910-0!!

Seal* 0.000*000 1.025*005 •5.O21-01! 2.910-0!!

P«i> !9T2 M«y 05 1124) «t 17:50:51 l!.l2« I. ! 024 P«iM». F»m« !!

Figure 9. The first frequency block for the relative phase between the

current and voltage of the detector.

36
h

AMPLITUDE VPHTR, I CASE

1PART

Sc.lt

Mi*

0.000*000

0.000*000

t. 025*005

1.025*005

Mil. 0.000*000

0.000*000

5.055*000

5.055*000

Pgft 1972 H« y 05 (124) 0| 17:50:59 It, 12* I. 1024 ••»!». ft— \2

Figure 10. The first frequency block for the LOG amplitude of the

voltage across the junction without z .

37
POWER I2V0UT SSTR, ICASE= 6 tPART

X • Hit, M#i I • Mi» M«i

r w«ct i«n 0.000*000 t.023«003 0.000*000 4.t%*000

Seal* 0.000*000 1.023*003 0.000*000 4.t?t«000

By* !9?2 H«y 03 (124) «i 17:50:53 Ll.lt" I. 1024 P|<ut, *'— 13

Figure 11. The first frequency block for the LOG

Q

power into the

detector.

38
•

•

•

• *

• •••••

•

•

•

• t 1

1

•

•

• ••••••

POWER ItVOUT SSTR, I CASE = 6 1PART

X • Mi* M«, I • HI* N«.

r u»ctit« 0.090*000 t . 023*005 -4.410*000 5.911*000

Seal* 0.000*000 1.025*00! -4.410*000 5.911*000

Rw« 1972 M«y 05 (124) «t 17:51:00 U.l2« 1. 1024 ?••»!•. '•#•• 14

Figure 12. The first frequency block for the LOG^ power into the

Josephs on junction.

39
AMPLITUDE I2TRAN, ICASE= 2PART

Mi*

0.000*000

0.000*000

1.023*003

1.023*003

• Mi*

0.000*000

0.000*000

4.??7*000

4.???*000

By* !9?2 H« y C3 (1241 «t 17:51:01 L1.L2" I. !024P»i»U. *'— 15

Figure 13. The second frequency block (128 £ l< 256) for the LOGjq

current through the detector.

40
AMPLITUDE VOUTTR, ICASE= 2PART

Sc«!#

o.ooo*oco

0.000*000

1. 023*003

t. 023*003

0.000*000

0.000*000

4.2?N0O0

4. 2^1*000

Ru» !9?2 H«v C3 (124) «t 17:51 : 01 Ll.l2* !. t024 P#im». fn !*S

Figure 14. The second frequency block for the LOG^ voltage across the

detector.

41
I L

i r 1^1

phase 12-vout tran icase*

Sc«!«

Kin 0.000*000

0.000*000

M«. 1.025*005

1.025*005

i • Mi»

-2.910-011

-2. 910-0!

I

2PART

M«,

5. 921-01

1

5.121-011

Ran 19?2 M«> 05 (124) «t 17:51:02 U.L2» I. 1024 Sim. Ff«*t 1'

Figure 15. The second frequency block for the relative phase.

42
AMPLITUDE VPHTR, ICASE=

FwftCtiffH Sc«!t

Mi.

0.000*000

0.000*000

M«.

I. 025*005

1.025*005

Mi«

0.000*000

0.000*000

2PART

4. 271*000

4.271,000

Ry» !9T2 M«y 05 1124) «t 17:51:05 U.l2« 1. 1024 ^iia, *'•*• II

Figure 16. The second frequency block for the LOG voltage across

junction without z .

43
POWER I2V0UT SSTR, ICASE= 6 2PART

r u«cti«« 0.000*000 t. 029*003 0.000*000 4.549*000

Scale 0.000*000 !. 025*005 0.000*000 4.549*000

»»» 1972 H#y 05 1124) at 17:51:05 U.l2« I, 1024 ••i»i». f»«M 19

Figure 17. The second 128 lines for LOG.q power into detector.

44
POWER I1VOUT SSTR, ICASE=

St.lt

0.000*000

0.000*000

M«i

1. 023*003

1.023*00!

Hi*

-4.390*000

4.H2*000

2PART

M«. 4.962*000

4.K2*000

»w» 1972 M«y 03 (124) «t 17:51:04 U.l2« t. 1024 Nun, *'— 20

Figure 18. The second 128 lines for LOG power out of junction.

(Noise pulses are seen as energy sources here. ) 45
(Here we have compared the length of the 8th harmonic to the 16th harmonic which is given as

full scale of (4. 410). We take the absolute value in this effort. The minus sign is simply a code to

tell us that the power is out instead of in. The 132th harmonic has value (5. 5 X (4. 350)/ 6. 2 = 3. 85).

The 133rd harmonic has value (3. 5 X (4. 963)/ 7. 1 = 2. 45). The above values are extracted from the

plots and the given scales in the plots. To complete the determination of the power we need to look at the conversion table of figure 1 9 to find that the positive 4th harmonic is reduced by -2, the negative

8th is reduced by -6, the positive 133rd is reduced by -24, and the negative 132 is reduced by -16.

Thus the power is given as Harmonic Power

4 10

3.918

8 -lo2 - 66

132 - 10

- 12' 15

133 lo"

21 ' 55

To scale these numbers to watts just multiply them by the 4. 8830 X 10 . Notice that the noise

signal in the 133 harmonic is approximately 9 decades less than the corresponding non- noise signal

of the 132 harmonic.

It should be obvious that if one wants extensive detailed tables of the various spectra, it is better to put in appropriate print statements into the programs. Our purpose has been mainly qualitative visual information of the power spectra. Figure 19 has a table of power and other information at selected harmonics. These numbers are in the dimensionless form used in the calculation. To convert to watts, volts, amperes, etc. , it is necessary to use the constants given in figure 5. We finish this section by describing the remaining details shown in figure 19. The first row

of data gives the total steady state power across the junction without the Z part (POWJ), the total steady state power across the full junction (POW1), the power across the detector (POW2), the power

put in by the applied voltages (POWIN), and finally the effective impedance KO. The next block of six

different row pairs give the harmonic number to be applied to the row pair information (KEY), the

applied voltages (B), the net initial phase at s = (PHC), the power in (POWIN), the power across

the junction (POW1), and the power across the detector (POW2). The second row of a given row pair

has the dimensionless values for the complex number impedances of that frequency. The frequencies

chosen to be printed out are d. c. , the direct beat frequency, the lower applied r. f. frequency, the

upper applied r. f. frequency, the highest frequency used as a consequence of the number of data

points computed, and finally the rotation frequency.

The remaining undiscussed information has the net d. c. voltage across, and the current through the

junction without the Z term, the data summarizing the error controls used in the program (EP1,

EP2, EP3) and finally the last computed P(a) set EP4(I6) . Please look at the program to see what

these errors limits mean in detail.

6. THE SUMMARY

In this report we are documenting the computer program, the mathematical techniques and

approximations used, and the procedures for use and interpretation of the computer program. The

46
»9»j.»e«l.»0M »MI*.i

1 'IMIt • *•!•• •i I "11104 » MUM* 1 llt'HI

»t>.».»«c.»oim »0»l.»0 t. /t.n.n

1 1 »».00. •4 '«1*11 I m>i t 4 14t.|9I * Ill-tit

i i»».ui i 900. 00 1 lll.lll 1 •*»•••» I ••••••I

• IH. 101 II 1*89 i s no-no « 09» 0*0

S 1M-IU 1 000.00 1 JOO.000 1 llt-OII 1 ••••199

4 I 411.111 •8 01 • •II ' I tit. 90S • fll-194

i tn.iii i 000.00 1. OO0.I0B 1 tit. I9t 9 ••9-919

* -0 199.091 •0 01 9.00 1 •! 119-911 t tlil<1

i tn-iii i 000-00 1 000.90! t IttOO* 9 990.990 :»M SIB. 01! 1. •• • •91 1 -4 ?4t-IM 1, MI-US

I tM.OIt C 080.00 • 390.010 I moot 9 111.999

4 011.10! 1. II • •99 C 919.0 • • IM.IIf * Wl-104

5 trt.oo* : 010.00 999*999 111.999 1 1«*00t 9 o»o. so:

.J'JH.lOtTJi . IJVMIIMItHSl

» 00099-114 4 Il'd-IH

(MO* ii«l'J I»!.t»| (M

CMfSi

! Ill-ill * 141-914 1 Mt-OOt

• *«»-9tl I 009.000 9 991.119 9 999*999 999.999 9 999*190

SCKI CI* OHi* DC Mi '€». « ln.mi ti* »:• i»c« 4»l»»

. T .? • -* .? -•»

: 9 • 1 9 9 -1 -It -10 -II -11 -II it to

: 9 • 1 9 9 9 1 -4 -1 -It -It •It It it -It

• • • • I • • • • t • • • I i I

• • • • • • • • -4 •• -It -It It It it It

• • • 1 • • • • -I II 1! II II -It it II

* 14 14 It - n -it it 14 -1 14 It It it n it It

Figure 19. Scaling information for a particular steady state run.

47
results using this program are to be shown in other reports.

Here the discussion is directed toward the procedures used to check the correctness of the

4 computer program. The strict d. c. cases have been selectively compared with the Stewart voltagecurrent model. The r. f. applied voltage case has been selectively compared with the Fack, et al, model and their reported results. Both cases agree within the accuracy of the techniques and the

graphs of the data reported in those papers. The final check has been a comparison with the voltage

driven model. Again the results are in agreement. The integral equation part has not been checked

against any other model or results since thereare no such published results. In terms of the computer

program, the program has been checked to be sure that the D is correctly generated. The only

remaining part for possible programming errors is in the integral equation section in subroutine PART.

The behavior and visual check shows no obvious errors. A detailed quantitative test is not known.

The most important numerical conclusion was to establish that it is very necessary to fix the

cp (o) and to allow the V (o) to be adjusted when one wants to pin the rotation frequency. If one wants

to use only time domain information, then one should use the much simpler programs such as used by

Fack, Stewart, and is used in the appendix and fix b instead of the rotation frequency. If one is interested in high order harmonic mixing, and detailed harmonic spectra, then this program is best.

If noise effects on the system of equations are of interest, then it is necessary to use the time sequence

procedures and to allow selected Fourier transforms on the various time series. That analysis will

be discussed in another report.

Because the analysis has been done for discrete uniformly spaced time points, we have discrete

induced frequencies that are rationally related to each other. The consequences of this situation are

additional steps induced in the d. c. current-voltage plots. The excursion of those steps is reduced

by increasing the number of time points used. To get an idea on the range of any d. c. step for a given

rotation frequency, it is necessary to make at least four calculations. First we chose a particular

PHA(l) = cp and compute the B for that case. If the program converges, then we have one point for

that case. To get a convergence case, we use the transient program with the chosen B and deduce an

adequate PHA(l). The steady state program just improves the accuracy of the results. Second, we

compute a B with PHA(l) now cp ± — . At least one of these will converge. If both converge, then

(B + B )/ 2 is the center and |(B - B )| is the range of the step. For example, if the B (B ) does not T - T - T - converge because the phase location denies steady state solution, then it is necessary to compute

cp + —, (cp - —) case to get the B _i , (B j_) value. We assume that the range of a step is defined by the

following simple formulae-- B = X + R sin 9, where X is the center of the step, R is the range and 9 is

the phase relation relative to the center of step. Normally 9 is not equal to cp and normally the formulae is not strictly accurate. It is only true if the dominate processes creating a step are due to a

single frequencing mixing. For estimating the range's order of magnitude, we can accept this approximate formulae to get the range and center of the step for the case considered, namely B (B ), B

respectively are given by

B = X + R sin 9

B_i = X + R sin (9 =F tt/ 4) B = X + R sin (9 =F tt/2)

48
X= (B, + B^ -"JT B_i)/(2 -yfT)

1 + +2

R = [(B

] - X)

2 + (B^ - X)2 ]'

B

l " X

— = sin 9 or 6 is gotten to give center. For more

accurate analysis, it is necessary to do empirical work.

7. APPENDIX

In this section we discuss briefly a simpler computer program called a transient program that can be used to get improved estimates for the effects of r. f. and d. c. applied voltages to the Josephson equations without the operator D. The equations used are closely related to those given in references 1 and 4. We use the trapezoidal spline method with variable step size so that the accuracy

per step is kept uniform throughout. Figure 20 gives the Fortran listing for this program. This

program simply generates a time series of selected points. This data can then be used to estimate the average rotation frequency, N , and to fix a PHA(l), VA(1), and a U.

7. 1 The Transient Equation

Here we show the equations used in the program and describe how they relate to equation 11. The transient equations are: *£= v

ds —dV — — —

= -v^V - a sincp+ B+ B T cos (NT s + 6 T )+ B TT cos (NTT s + 9 TT ) (Al) QS 1 Li Lt L U U U

Here cp = cp + cp + cp , a b c

where dcp dV

HT = V ST = " v T

V

b

+ *l cos (nl

s + 9 l

)+ *u cos ( V+ V '

and cp is the initial value of cp at s = 0; thus here cp = at s = 0. We note that B - b /K_, B t = b T /K.,

c a a O Li L, and B = b../K . The best method to fix a, v_,, and K is to use the main program first for a given

C, R, Z., a , J, etc. , with a d. c. case. We then use this program called the transient calculation for

the various B , B , 9 , 9 cases. Please note that the scale used for b , h necessarily gives

Li U Li U Li U

different phases and amplitudes in the main program if the K(N ), K(N ) are different from K . Li U \J When the programs are used, this difference needs to be considered. Before we look at the input and

output of this program, we briefly consider the criteria for the step size. It simply restricts the

absolute change in cp to be less than STEP/ (3V1-V2) and the absolute change in V to be less than

STEP'W [v T

(3Ql-Q2)]. The STEP is the control accuracy for this allowed change per step; the VI

dV and V2 are the last two v

1 s used; the Ql and Q2 are the last two q's used. The q is -— . as

49
7. 2 The Use of the Transient Program

This program is organized under the assumption that a given computer run will fix the r. f. and

selected control parameters and will sweep the applied d. c. current B over many values. Thus the

first read statement has a parameter which fixed the range of s. The initial V, and PH has no direct

correspondence to the expected output VA and PHA, hence they are used to just start the calculation.

The IPRINT defines how many steps are skipped before the usual print mode takes place. Figure 21 shows

the flow of input data. Figure 22 gives an example of printed output data. The first three rows of

printed output simply gives the summary of the input data use in this particular computer run. The

list of names S, VA, etc. , are the labels for the six columns of numbers that follow the Warning

statement. The S is the time variable. The VA is the cp& with the average steady state slope of the

ramp term not removed. PH is the total phase cp. PHB is the r. f. phase including the initial PHC.

VR is the integral of V; it is used to estimate the slope of the ramp. PHA is the cp a phase with the

ramp term still present. To illustrate how to get the necessary initial data for the main program

from this transient program, we run through the sample output to get the proper numbers. Please

note that the printed output far exceeds the necessary information for these variables. It may be

possible to suppress all printing except the few numbers of interest. The current calculation assumes

that the transients are unimportant after S exceeds 2tt, hence we let the computation only run about one

more 2tt interval. This last interval is used to estimate the ramp and the needed initial PHA(l), VA(1),

and U in the main program. Note that Pl6, Vl6, and BI are respectively used as the input variables

for this initial data. In our sample output data we have the following needed data:

s VA PHB VR PHA

6.28245 8.665 -.756 20.74 21.486

12. 57033 8.933 -.645 46.03 46.66

B = 102, and v = 32. We use subscript 1 and 2 to label each data now. We chose these points

because we want PHB to be nearly the same in both this program and in the main program. If you

change to other THL and THU then you can consider other data points. Keep in mind that it is not

necessary to have accurate values. The purpose is to get a guess that is at least accurate to 10 to 20%

so that the circle of convergence for the quasilinear procedure will operate effectively in the main

program and so that some information is available to estimate potential effects due to spiking. We

can then choose the necessary accuracy in the main calculation. The ramp slope is estimated as

N.J = (VR

2

- VRj)/(s

2

- Bj) + 4.

The initial data is thus estimated as

BI = B - V T

*N = -26

VI 6 = (VA + VA

2 )/2. - Nj = 4. 8

PI 6 = (PHA - N *2tt)+ (PHA

2

- Nj*4n) / 2.

= (PHA + PHA /2. - N *3tt = -3.626 . 1 Ct I

50
PROGRAM JOSPA

C TO GET INTIAL WAVE SHAPE FOR STARTING GUESSES IN JOSPHC

C PHDOT=V $ VDOT=B-V*ANUT-A*SINF(PH>+BL*COSF<NL*S+THL>+

C BU*COSF(NU*S+THU)

COMMON /l/ B.ANUT.BL.THL.BU.THU.NL.NU.A.STEP.STEP1

READ 1.RANGE»V.PH,BL.BU,THL.THU»STEP.ANUT,A.IPRINT,NL.NU

70 READ l.B

1 F0RMAT(5E12.3/5E12. 3/3110)

IF(EOF»60)71»72

71 CALL EXIT

72 PRINT 2.B.V.PH.BL.BU.THU.THL.STEP.ANUT.A.RANGE.IPRINT »NL»NU

2 F0RMAT(* B.V.PH.BL.BU/THU.THL.STEP.ANUT.A.RANGE/IPRINT .NL,Nu *

1 /5E12.3/6E12. 3/3110/* S.VA.PH.PHB. VR*PHA * //)

PRINT 99

99 FORMAT <* WARNING THE VA AND PHA HAVE NOT HAD THE RAMP TERMS REMOVE

ID *)

J=0 $STEP1=STEP*A/ANUT$IF(ABSF(V)-1.)26.26.25

26 H=STEP/100. $GO TO 270

25 H=STEP/( 100.*ABSF(V)

270 PHC=PH $ S=0.

Il=-2$RANGl=RANGE-7.

27 11=11+1 $ IF<S-RANGE)4,4»70

40 J=J-1 $IF(J)41»41»3

41 IF(S-RANG1>3. 417.417

417 IF(S)410»410,412

410 IF<NL*NU>70»70»415

415 Xl=-BL/(NL*SQRTF(NL*NL+ANUT*ANUT)

X2=-BU/(NU*SQRTF<NU*NU+ANUT*ANUT)

APH1=ATANF( ANUT/ND +THL $X3=-X1*NL

APH2=ATANF(ANUT/NU) +THU $ X4=-X2*NU

412 S8=NL*S+APH1 $ S9=NU*S +APH2

PHB=X1 *COSF ( S8 > +X2*C0SF ( S9 > VB=X3*SINF(S8)+X4*SINF(S9)

VA=V-VB

PHA=PH-PHB -PHC

413 PRINT 42.S.VA.PH.PHB.VR »PHA $J=IPRINT

3 CONTINUE S GO TO 27

42 F0RMAT(7E15«8)

4 IF( Il>14»15»16

14 V2=V $ PH2=PH $S=S2=0. $ CALL START < PH. S»QX ) $OT=QX-ANUT*V $60 TO

1 410

15 PHY=PH2+H*V2 SCALL COM ( PH.V.Q .S .PH2 .V2 .QX .0. .H.PHY)

VR=V*H

QT=Q-ANUT*V £ PH1=PH$ V1=V $Q2=QX $Q1=Q $S1=S $GO TO 40

16 CALL PR0J(PHY,H1.H»V2»V1»Q2»Q1»PH1)

CALL C0M(PH.V.Q,S,PH1»V1.Q1.S1.H1»PHY) $QT=Q-ANUT*V

DL1=PH-PH1 $DL2=PH1-PH2

Q2=Q1 $V2=V1 $PH2=PH1 $S2=S1

VR=VR+V*H1

Q1=Q $ V1=V $ PH1=PH $S1=S $H=H1 $IF (DL1*DL2 ) 412 .412 .40

END

SUBROUTINE COM(PH,V»Q.S.PHX . VX »QX .SX.HX.PHY)

COMMON /l/ B.ANUT.BL.THL.BU.THU.NL.NU.A.STEP.STEPl

DATA(ERO=l.E-5)

H=HX

H3=H*ANUT/2. $H4=1.-H3 $H5=1.+H3 $S=SX+HX

CALL QR(O.QY.QZ.PHY) $ H6= { H*H*.25 ) /H5 $H9=H4/H5 $H10=H/ ( 2.*H5

H7=1./<1.-H6*QZ) $H8=VX*H/H5+PHX+H6*QX $QZ1=QF(S) $QY=QY+QZ1

IV=5

4 IV=IV-1 $ IF( IV)5.6»6

6 DEL=(H8-PHY+H6*QY)*H7

PHY=PHY+DEL SCALL QR ( 1 »QY.O. .PHY ) SQY=QY+QZ1

IF(ABSF(DEL>-ERO)8»4»4

5 PRINT 7.S SCALL EXIT

7 FORMAT (* NO CONVERGENCE *.E12.3)

8 PH=PHY $Q=QY $V=(H9*VX+H10*(QY+QX) ) $RETURN $END

Figure 20. The transient Fortran program listing.

51
SUBROUTINE PROJ (PHY. HI ,H »V2 » VI ,Q2 *Ql »PH1

)

COMMON /l/ B»ANUT,BL»THL»BU»THU.NL.NU»A*STEP.STEPl

H1A=STEP/ABSF(3.*V1-V2)

H1B=STEP1/ABSF(3.*Q1-Q2>

H1=MIN1F(H1A,H1B)

PHY=PH1+H1*(VH-(V1-V2)*H1/(2.*H) ) SRETURN SEND

SUBROUTINE START ( PHX »SX »QX

)

COMMON /l/ B.ANUT,BL»THL»BU»THU»NL»NU»A»STEP»STEP1

CALL QR( 1»QX»0. .PHX) $QX=QX+QF ( SX > SRETURN SEND

FUNCTION QF(Z)

COMMON /l/ B.ANUT.BL.THL.BU.THU.NL.NU.A.STEP.STEP1

QF=B+BL*COSF(NL*Z+THL)+BU*COSF(NU*Z+THU) SRETURN SEND

SUBROUTINE QR( I .QA.QB.X)

COMMON /l/ B»ANUT.BL»THL»BU»THU»NL»NU»A.STEP»STEPl

IF( Dl»l»3

1 QB=-A*COSF(X)

3 QA=-A*SINF(X) SRETURNSEND

Figure 20. (Continued)

52
PROGRAM JOSPA

READ 1 , RANGE, ETC

70 READ 1 , B IF(EOF) EXIT

IF(S-RANGE)4,4,70

COMPLETE

PRINT

SUBROUTINES CALLED

COM

PROJ

START

QF

QR

Figure 21. The flow chart for Josephson transient program.

53 IB-in?«TM I
o

1-5

o CM 3 a

O o

+• a o

o U3 o U> 3 • • Z o M

•> 1 _Jz

~CM t4 a a

t— o n> Z 4- *

1-4 o o

C£ * (3 a. vo CM

(-4 • •

•» is. rO Ul

19z

<t o t-I « O O

N. O a

< 4- l •»(M o

t- J" o

=D -4 O Lfl Z . • <t rO r4

•> 1 0-

txJ h- O o

71 O o

» CD CD

_l 4- X o o J- «s

i— o CD X

«. CD 3 0. 3 • • •fc X -H O a:

r- 1 •>

s*

"3 •* ro cm CD CD

•» CD O T

_J o a o a. rn 4. 4- CO *. •. o CD X

X CM o n

0. o C3 •% * • • <i > -* o r>

t4O00Ot4»4t4t4»-|t4t4t-It4-4»-I«4t-It4t4-4t4t4t4t-ItHt-It4t4:t4t4 ooooooocdcd 3000000000000000000000 cdcdoocdoooooooooooooooooooooooooo I + + + +"+ '4>4>*4>4>+ + + 4>4> + 4> + 4>4> + ++ + 4>4> + 4>4> + roo3T403^r^uNN«oin'^fococMrv7cnLnincMin^:ocncors..»,-iT4T4incM ^-*orCN»^OiMX»rOeOeO»X>r\J-*Ors.vJDinCnrs.cO^-yeOT4<MiX>in\J30

coco-4CT>03ro.J> 03rs.invJ3X<v£»o3in^^rorT5CMCMT4^T4T4C\jCMCMiO»0.*

r^OOOOOOOOOOOOOOOOOOOOOOOOOOOOOO ooooooooooooooooooooooooocioonoo 4>4>4>4>4»4.*4.4.4>4.4- + 4«4» + 4.4>4.4..*.4»*4. + 4,4> + 4>4» +

oininoio'ocMroo^ooo303<7<OT4coininc7»c\j^rM.t, <»^T4rs,o3conjo

ororo*Hin<t3T4»-4co<<>rs.aors.^PO-*T4CMcOT4or^CMXt'T4XtCMT4"Oo(M oXt^njT40*CT>e003T4\XT^030rOgJeO<T>T4rgr\jroro.^.»^vOfs.O, '"» Ol^oS»03»»7 oin^eOT4lT3a3cOT4rOinN.CnrM.D*vXaOO<M-*vOcOOCMin

o (/I

o

1 CD

CM CM CM —

(

o CD o CD CD O o o o CD o a CD CD CD rs o o CD CD CD CD CD o o CD o CD T-t a a CD r> CD CD CD -3 CD CD CD ^3 CD CD o CD CD C3 -> o CD CD CD CD CD CD CD CD CD CD CD CD

i-j CD ID O C3 CD CD CD a CD CD CD CD o CD CD CD CD o CD CD CD CD CD CD CD CD CD O CD CD CD > 1 4. • + 4- 4- 4. 4- 4- 4- 4- 4- 4- 4- 4- +• 4- 4- 4- 4- + 4- 4- 4- 4- 4- 4- 4- 4- 1 o CO ro IT N» CD O 3^ CM CD ro CM f>- ro •J- r^ •O CM cn •S cr> 00 J" CD -r4 in ro ro ro fs. CO cn X ro X> CO J- J- J- CO CNJ in r>~ 03 ro -o "i CD CM in CD in CD S. O CO S. N- CM m in ro fs. j- uj W T-4 ^ ro J- o> J- J" T-« (M CD 03 ro CD in *o -•* cn CD r> t-t CO -J" N- j- t4 ro ao CT» ro o

ct t4 tH 33 r». ro tH J- N. J" ro cr» m <T> CP j- •r* CO J- a CO ^ ro -* T-l x> 03 cn 03 en is. rs. o CO in T> ro ro in ni CJ no N- CM ro T-^ CD ^ J- •o CD CM in CM T-l CO 03 C3 o CD j- CD CD

c/> J- J- 3D ~t (T> ro ro ro CNJ cr r^ ^O Cf in -J" 03 H in CM o CP <^ J" CT> in CD T-| CO CO t4 "t

5- 33 -o CD N. CO ir\ r^ in N- m CO (^ CM m vO >o rs. CO •H in CD CO 03 in J3 CO CD CM in cn -J"

UJ

"O -J- •o CO ro cr« ro cr> ro 03 j- CM '3 is. J" t4 T3 JN ro CD CO in ro -i CT> N- 03 J- CM o CM

is. if\ LTl LTi IT> LTi in in in in in m in j- J- J- o -o ro ro CM ,\j CM 00 t4 tH tH T-l t4 tH cn

• 1 1 l i I 1 t i i 1 1 I 1 1 I 1 • 1 1 1 t 1 1 1 1 I <D_ cd TD n CD CD •F-t CD X o ^5 CD CD CD CD O o o a CD CD CD CD CD CD CD CD CD o O a a CD CD O CD CD -D CD CD < O "3 CD CD '^D CD CD CD CD CD CD CD CD i-D CD CD CD CD CD CD a ~3 CD CD a CD CD CD CD CD CD a: 4- 4- 4. * 4- 4- * 4- 4- 4- 4- 4- 4- -». 4- 4- 4- 4- 4> 4- 4- 4- * 4- 4- 4- 4- 4- 4- 4- 4. o ro J- ITN C^ ro <M tH CO CM CD ro ro ro (V ro CD in in in N. vO CT> 03 T-l X) CT> 03 •T J- CD UJ CD r\j f^. —

<

H s. -o CM CO r\J r> N. J> CM -s. 03 T-| .0 T-l 00 tH in J- CD ro ro is. CM CM t4 in X o x> r»- (>J r^ r^ CD *-{ 'M »H CD ~3 CO J- CD CD CD CM 03 CM CM in CD J- "O o "O CD 'X CD <a i_ m >\J •o 1^ CO T« f^ CD •o (M CM T-4 -i- J- * O^ r\l T -o CO 03 CD T> m CT> >o a> • o m so J- in J- J- •o CO O fO 03 CM fM T-l rv) •H 03 in fs. J- J" CD fM T-l on CM •O CD in t4 CD AJ CO CD o T-t T-4 ~vj N. X> .3 CT« ao m rg N. CM »^. -H in CO T-l ro in 03 N. N. CO 33 CP C" CJ t4 M o o

<X J" ro ro (M ^O (\l 03 -D* (M r^ CT> CM ^t IS. CT> T-t J" 03 CO o CM J" 03 CO a oo in rs. o> t4 ro X -* CD N. vO rsj rO <T cr <H rM <M O -O ro -o J- J" -1- J- in >n m in in 03 03 03 03 >o N» N- o

z

rO

1-4

Jt CD

r«. »4

•

*-* J- T-t T-» t4 T-t T-l -H t4 tH •H T-| T-l -H •H T-l •rH H T-t -4 t4 t-4 t4 t4 t4 -4

CD CD -H -H -^ -4 CD CD tH T-l CD CD CD o CD C3 CD CD o -H T-l CD CD CD CD CD O O

CD C5 r-> rD CD f-D cr> CD CD rD --D m ~5 TO -5 C» m CD —

>

CD CD .-> CD CD CD -D -1 CD CD -1 -~> UJ O a CD o CD CD o CD CD o CD CD CD CD CD CD CD =3 CD CD CD CD CD CD o CD CD CD CD a CD > + + ¥ 4- 4- 1 4- I 4- 4- 4- 1 1 4- 4- 4- 4- 4- 4- 4- + 4- 1 1 4- 4- 4- 4- 4- 4- 4- <t CO ^H CD <T> J- rv. »\J N. in J- fM J- rr> CO rs. CT> ^ in ^ in O .D OO 0"> O CM IS. J" •X T-l CD X cn -D ro CD •o CO r\l CO j- "O —

1

n fv. Jt ro \C\ n> _r vO <M T-t m n o m J- in r\i IS. in

03 ro CO •H ^D CM r-( a r<^ 03 03 CD 03 ro CD cn tH in rv. rs. N- cn in ro CC 00 J- CO in CO ro el N- ro CD CO JD CD cn N. CM ro T-l J" T-» 'CD r-. J- CM CO \S\ ro 03 t4 ro 03 ro ro co T-l CD ro t4

T" l-S n <n tO r«. N. iO n 1^ ro f> CM r^ ro ro 'T. -H CM IJ-N T-l O (-J J- 00 in CD J- rs. ao f\i —

t

CL "O CO ro N* CO J- nj ~r (M ^ N. in ro J" CO -^ J" N. J" CM in r«- CO CO ro X ro •O ro m ro

x> LO (M 0D ro cn m 03 in (M O CD CD •o CD CD CM CD cn ro 03 ro ^0 CO ro ro CM 03 CM 03 03 Q

Z

CM CD (^ -H CD CP CD !X\ o 03 CO CO T-l T-l CM a^ CM CM CO ro in 03 CO J- K- cr t4 CM rO CM O

CM & J- J- •r-t •XI •^ 03 tH J- CM CP r^ CM ro ro J- J" ro "O CM t4 in >x» tH CM J" in 03 n. CO

1 1 1 1 1 1 1 I 1 t I I I • 1 I <I CD H a CD CD CD a CD CD CD CD CD CD O CD CD CD CD CD CD CD CD CD CD CD CD CD CD CD o CD > CD o o CD CD o o CD o CD CD CD CD CD CD CD CD CD CD CD CD CD CD CD CD O CD o CD CD CD

13 CD CD CD CD CD CD CD CD CD O CD CD CD CD CD CD CD CD CD CD CD CD O C3 o CD CD CD CD CD UJ 4- 1 * 4- 4- 4- + + * 4- 4- 4. + + 4- * + 4- 4- 4- 4- 4- 4- 4- 4- 4- 4- 4> 4- 4- + I CD CD tn ro LT» CO <NJ hs. •r-( m ro r^ -4 O J- J- CO CD 03 CM 03 ^ CM ro CD rs. J- CD rO X) CM r- CD cn ro CO •H cn m ^4 ro CO CP CD a J- CD 03 J- cr> 03 J" 03 ro a> tH rs. 03 03 rs. CO J- J- cd ao o^ T-« ro r» cr J" cr> in T-l J- fs* H in 00 CP N. ro T-l T-l CM 0~ -T in ,o T-l rs. X) 33 CM LD o >£ CO CP IT< s. m CD J- j- ro cn -r cn "O CO CO CD CP in —

1

T-l CO 03 m 03 CO t4 CM CP T-l Z CD iD -t- ao oo ^n _» ro •T-l jf ~4 en N. rr> ro J" O.I fs, N. 03 J" CO T-l J- -^ ro T-l s. CD t4 fj

>—

t

CD 05 J- •o o <M cr> ro CD CD T-l T-l <~0 CM CM CD N. CM 03 CP T-| 00 CM t4 CD CO 03 rO t4 CO in z CD ro vO in N. JD ro ro T-l CD CM J- 03 CO o CM ro m 03 N- cr o T-t CM ro ro J- m X 03 fs. £* CD <\j rg <r< ^ l^ (M •r-l CO CD O CD CD CD T-l T-l T-l T-l T-t t4 T-l CM CM CM CM 00 OJ CM CM CM CM

cd ITi •^-( —* 'M ro -3* in in 03 03 03 03 03 03 03 03 03 •X) 03 03 03 03 03 vO D 03 X) 03 D X>

cdu

CUD

ft 3:

0)

• H

co 1 a

(TJ u

4->

X3

U

<+H

o

cd

54
Note that only integer ft can be used in the main program and that the choice of B in the transient

program does not necessarily correspond to such integers.

8. ACKNOWLEDGMENTS

We wish to thank Stephen Jarvis and Clark Hamilton for useful discussions.

9. REFERENCES

1. W. C. Stewart, App. Phys. Letters 1_2, 8(1968), pp. 277-280.

2. Theory and Application of Spline Functions , Ed. T. Greville, Academic Press (1969). An

Introduction to the Application of Spline Functions to Initial Value Problems, Frank R. Loscalzo,

pp. 37-64.

3. A. Miele, "Modified Quasilinearization Method for Solving Nonlinear, Two Point Boundary Value

Problems," (Rice University, 1970), Clearinghouse Reports N71-1779.

4. H. Fack, V. Kose, andH. -J. Schrader, "Calculation of Characteristic Curves for Point- ContactJosephson- Elements Subjected to High Frequency Radiation," Messtechnik (February 1971) pp. 31-38.

55
FORM NBS-114A (1-71) U.S. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET

1. PUBLICATION OR REPORT NO.

NBS TN-627

2. Gov't Accession No.3. Recipient's Accession No.

4. TITLE AND SUBTITLE

Computation of Spectral Data for a Josephson

Junction Circuit

5. Publication Date November 1972

6. Performing Organization Coc

7. AUTHOR(S)

E. G. Johnson, Jr. and D. G. McDonald

8. Performing Organization

10. Project/Task/Work Unit N, 2710164

9. PERFORMING ORGANIZATION NAME AND ADDRESS

NATIONAL BUREAU OF STANDARDS, Boulder Labs

DEPARTMENT OF COMMERCE

Boulder, CO 80302

11. Contract/Grant No.

12. Sponsoring Organization Name and Address

Same as 9 above.

13. Type of Report & Period

Covered

Final (4/72-10/72)

14. Sponsoring Agency Code

15. SUPPLEMENTARY NOTES

16. ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a significant

bibliography or literature survey, mention it here.) A computer program has been developed to study power flow between differem

frequency channels in a Josephson junction circuit. This paper discusses the

mathematical assumptions used to get such results. They are the trapezoidal

approximation from spline theory and the use of a finite range of frequencies

characterize the frequency spectrum. This paper describes the program and

provides the FORTRAN listing, flow charts, and discusses how to use the

program. A discussion of possible sources of errors is also included.

17. KEY WORDS (Alphabetical order, separated by semicolons)

Differential equation; fast Fourier transform; Josephson junction; nonlinear

integral-differential equation; spline theory.

(

1

I

h

h

o

N

18. AVAILABILITY STATEMENT

PC] UNLIMITED.

|~~| FOR OFFICIAL DISTRIBUTION. DO NOT RELEASE

TO NTIS.

19. SECURITY CLASS

(THIS REPORT)

UNCLASSIFIED

20. SECURITY CLASS

(THIS PAGE)

UNCLASSIFIED

21. NO. OF PA

63

22. Price

60£

Ti U.S. GOVERNMENT PRINTING OFFICE: 1972—784-207/1103 REGION NO 8 USCOMM-OC 662-
NBS TECHNICAL PUBLICATIONS

PERIODICALS JOURNAL OF RESEARCH reports National

Bureau of Standards research and development in

physics, mathematics, and chemistry. Comprehensive

j scientific papers give complete details of the work,

including laboratory data, experimental procedures,

and theoretical and mathematical analyses. Illustrated

jj with photographs, drawings, and charts. Includes

listings of other NBS papers as issued.

Published in two sections, available separately:

• Physics and Chemistry

Papers of interest primarily to scientists working in these fields. This section covers a broad range of physical and chemical research, with major emphasis

1 on standards of physical measurement, fundamental

constants, and properties of matter. Issued six times

j a year. Annual subscription: Domestic, $9.50; $2.25

additional for foreign mailing.

• Mathematical Sciences

Studies and compilations designed mainly for the

mathematician and theoretical physicist. Topics in mathematical statistics, theory of experiment design,

numerical analysis, theoretical physics and chemistry, logical design and programming of computers

and computer systems. Short numerical tables. Issued

quarterly. Annual subscription: Domestic, $5.00;

$1.25 additional for foreign mailing.

TECHNICAL NEWS BULLETIN

The best single source of information concerning the Bureau's measurement, research, developmental, co- operative, and publication activities, this monthly

publication is designed for the industry-oriented

individual whose daily work involves intimate contact with science and technology

for engineers, chemists,

physicists, research managers, product-development

managers, and company executives. Includes listing of

all NBS papers as issued. Annual subscription: Domestic, $3.00; $1.00 additional for foreign mailing.

Bibliographic Subscription Services

The following current-awareness and literature- survey bibliographies are issued periodically by

the Bureau: Cryogenic Data Center Current

Awareness Service (weekly), Liquefied Natural

Gas (quarterly), Superconducting Devices and

Materials (quarterly), and Electromagnetic

Metrology Current Awareness Service (monthly) . Available only from NBS Boulder Laboratories. Ordering and cost information may be

obtained from the Program Information Office,

National Bureau of Standards, Boulder, Colorado 80302.

NONPERIODICALS Applied Mathematics Series. Mathematical tables, manuals, and studies. Building Science Series. Research results, test methods, and performance criteria of building materials, components, systems, and structures. Handbooks. Recommended codes of engineering

and industrial practice (including safety codes) de- veloped in cooperation with interested industries,

professional organizations, and regulatory bodies. Special Publications. Proceedings of NBS conferences, bibliographies, annual reports, wall charts,

pamphlets, etc. Monographs. Major contributions to the technical

literature on various subjects related to the Bureau's

scientific and technical activities. National Standard Reference Data Series.

NSRDS provides quantitative data on the physical

and chemical properties of materials, compiled from

the world's literature and critically evaluated. Product Standards. Provide requirements for sizes, types, quality, and methods for testing various industrial products. These standards are developed co- operatively with interested Government and industry groups and provide the basis for common understanding of product characteristics for both buyers and

sellers. Their use is voluntary.

Technical Notes. This series consists of communications and reports (covering both other-agency and

NBS-sponsored work) of limited or transitory interest. Federal Information Processing Standards

Publications. This series is the official publication

within the Federal Government for information on

standards adopted and promulgated under the Public

Law 89—306, and Bureau of the Budget Circular A—86

entitled, Standardization of Data Elements and Codes

in Data Systems.

Consumer Information Series. Practical information, based on NBS research and experience, cover- ing areas of interest to the consumer. Easily understandable language and illustrations provide useful background knowledge for shopping in today's tech- nological marketplace.

CATALOGS OF NBS PUBLICATIONS NBS Special Publication 305, Publications of

the NBS. 1966-1967. When ordering, include

Catalog No. C13.10:305. Price $2.00; 50 cents additional for foreign mailing.

NBS Special Publication 305, Supplement 1, Publications of the NBS, 1968-1969. When ordering, include Catalog No. C13.10:305/Suppl. 1. Price

$4.50; $1.25 additional for foreign mailing.

NBS Special Publication 305, Supplement 2,

Publications of the NBS, 1970. Whenord r- ing, include Catalog No. C13.10: 305/SuppI. 2. Price $3.25; 85 cents additional for foreign mailing.

Order NBS publications (except Bibliographic Subscription Services)

from: Superintendent of Documents, Government Printing Office, Washington, D.C. 20402.
U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Washington, D.C. 20234

POSTAGE AND FEES PAID

U.S. DEPARTMENT OF COMMERCE OFFICIAL BUSINESS U.S.M

Penalty for Private Use, $300