Contents of file 'mopagrat.f':
1 c...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
2 c File: mopagrat.f (Fortran 77 source for the mopagrat routine)
3 c Author: Fredrik Jonsson <fj@phys.soton.ac.uk>
4 c Date: December 26, 2005
5 c Last change: January 3, 2006
6 c
7 c Description: The mopagrat routine calculates magneto-optically
8 c parametrically amplified field envelopes of forward
9 c and backward travelling wave components inside a Bragg
10 c grating of sinusoidal spatial distribution of the
11 c refractive index.
12 c
13 c Required external routines:
14 c
15 c ccubsolv - For solving the characteristic complex-valued cubic
16 c polynomial equation that appear for the eigenvalues
17 c of the system determinant of the wave propagation.
18 c
19 c zgaussj - For solving complex-valued linear systems of equations
20 c by Gauss-Jordan elimination.
21 c
22 c Copyright (c) 2006, Fredrik Jonsson
23 c...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
24 subroutine mopagrat(alpha,beta,gamma,delta,eta,kappa,ll,z,
25 + a1iplus,a2iplus,a1iminus,a2iminus,
26 + a1famplplus,a1bamplplus,a2famplplus,
27 + a1famplminus,a1bamplminus,a2famplminus,
28 + a1fplus,a1bplus,a2fplus,a1fminus,a1bminus,a2fminus)
29
30 real*8 alpha,beta,gamma,delta,eta,kappa,ll,z
31 complex*8 a1iplus,a2iplus,a1iminus,a2iminus
32 real*8 a1famplplus,a1bamplplus,a2famplplus
33 real*8 a1famplminus,a1bamplminus,a2famplminus
34 complex*8 a1fplus,a1bplus,a2fplus,a1fminus,a1bminus,a2fminus,c(3)
35 complex*8 bbplus(3,3),xiplus(3,3),lambdaplus(3),ccplus(3)
36 complex*8 bbminus(3,3),ximinus(3,3),lambdaminus(3),ccminus(3)
37 real*8 tmp
38 integer*4 k
39 external ccubsolv,zgaussj
40
41 c *********************************************************************
42 c Calculate the eigenvalues of the system matrix A of the magneto-
43 c optical parametric interaction. This is done by solving the cubic
44 c polynomial equation det(A-I*lambda)=0, where I=diag(1,...,1) denotes
45 c the identity matrix, or equivalently of the form
46 c
47 c lambda^3+c[2]*lambda^2+c[1]*lambda+c[0]=0,
48 c
49 c by using the external ccubsolv() routine. Here the coefficients are
50 c as follows:
51 c
52 c c(1)=c[0], c(2)=c[1], c(3)=c[2].
53 c
54 c After the call, the three complex-valued eigenvalues are returned
55 c in the vector lambda(1..3).
56 c *********************************************************************
57 c(1)=cmplx(0.0,((abs(kappa)**2)-((alpha+delta)**2))
58 + *(alpha+delta+beta+gamma)-((eta)**2)*(alpha+delta))
59 c(2)=cmplx(((alpha+delta)**2)-((kappa)**2)-((eta)**2),0.0)
60 c(3)=cmplx(0.0,-(alpha+delta+beta+gamma))
61 call ccubsolv(c,lambdaplus)
62 c(1)=cmplx(0.0,((abs(kappa)**2)-((alpha-delta)**2))
63 + *(alpha-delta+beta-gamma)-((eta)**2)*(alpha-delta))
64 c(2)=cmplx(((alpha-delta)**2)-((kappa)**2)-((eta)**2),0.0)
65 c(3)=cmplx(0.0,-(alpha-delta+beta-gamma))
66 call ccubsolv(c,lambdaminus)
67
68 c *********************************************************************
69 c Having calculated the eigenvalues of the system matrix A, now
70 c proceed with calculating the corresponding eigenvectors. The eigen-
71 c vectors are stored in the matrix xi, in which the (j,k):th element
72 c is the j:th element of the k:th eigenvector.
73 c *********************************************************************
74 do 10,k=1,3
75 xiplus(1,k)=cmplx(1.0,0.0)
76 xiplus(2,k)=cmplx(kappa,0.0)
77 xiplus(2,k)=xiplus(2,k)/(cmplx(0.0,alpha+delta)
78 + +lambdaplus(k))
79 xiplus(3,k)=(cmplx(0.0,alpha+delta)-lambdaplus(k))
80 xiplus(3,k)=xiplus(3,k)+cmplx((kappa**2),0.0)
81 + /(cmplx(0.0,alpha+delta)+lambdaplus(k))
82 xiplus(3,k)=xiplus(3,k)*cmplx(0.0,1.0/eta)
83 ximinus(1,k)=cmplx(1.0,0.0)
84 ximinus(2,k)=cmplx(kappa,0.0)
85 ximinus(2,k)=ximinus(2,k)/(cmplx(0.0,alpha-delta)
86 + +lambdaminus(k))
87 ximinus(3,k)=(cmplx(0.0,alpha-delta)-lambdaminus(k))
88 ximinus(3,k)=ximinus(3,k)+cmplx((kappa**2),0.0)
89 + /(cmplx(0.0,alpha-delta)+lambdaminus(k))
90 ximinus(3,k)=ximinus(3,k)*cmplx(0.0,1.0/eta)
91 10 continue
92
93 c *********************************************************************
94 c Calculate the inverse of the matrix bb by using the Gauss--Jordan
95 c elimination routine zgaussj(), as customized from the Numerical
96 c Recipes' real-valued version of the same algorithm.
97 c *********************************************************************
98 tmp=max(realpart(lambdaplus(1)),
99 + realpart(lambdaplus(2)),realpart(lambdaplus(3)))
100 do 20,k=1,3
101 bbplus(1,k)=xiplus(1,k)
102 bbplus(2,k)=xiplus(2,k)*
103 + cexp((lambdaplus(k)-cmplx(tmp,0.0))*cmplx(ll,0.0))
104 bbplus(3,k)=xiplus(3,k)
105 20 continue
106 ccplus(1)=a1iplus
107 ccplus(2)=cmplx(0.0,0.0)
108 ccplus(3)=a2iplus
109 call zgaussj(bbplus,3,3,ccplus,1,1)
110 tmp=max(realpart(lambdaminus(1)),
111 + realpart(lambdaminus(2)),realpart(lambdaminus(3)))
112 do 30,k=1,3
113 bbminus(1,k)=ximinus(1,k)
114 bbminus(2,k)=ximinus(2,k)*
115 + cexp((lambdaminus(k)-cmplx(tmp,0.0))*cmplx(ll,0.0))
116 bbminus(3,k)=ximinus(3,k)
117 30 continue
118 ccminus(1)=a1iminus
119 ccminus(2)=cmplx(0.0,0.0)
120 ccminus(3)=a2iminus
121 call zgaussj(bbminus,3,3,ccminus,1,1)
122
123 c *********************************************************************
124 c Calculate the logarithmic amplification of the transmitted and
125 c reflected signal and idler, in measures of dB.
126 c *********************************************************************
127 a1fplus=cmplx(0.0,0.0)
128 a1bplus=cmplx(0.0,0.0)
129 a2fplus=cmplx(0.0,0.0)
130 a1fminus=cmplx(0.0,0.0)
131 a1bminus=cmplx(0.0,0.0)
132 a2fminus=cmplx(0.0,0.0)
133 do 40,k=1,3
134 tmp=exp(realpart(lambdaplus(k))*ll)
135 a1fplus=a1fplus+ccplus(k)*xiplus(1,k)
136 + *cmplx(tmp*cos(imagpart(lambdaplus(k))*ll),
137 + tmp*sin(imagpart(lambdaplus(k))*ll))
138 a1bplus=a1bplus+ccplus(k)*xiplus(2,k)
139 a2fplus=a2fplus+ccplus(k)*xiplus(3,k)
140 + *cmplx(tmp*cos(imagpart(lambdaplus(k))*ll),
141 + tmp*sin(imagpart(lambdaplus(k))*ll))
142 tmp=exp(realpart(lambdaminus(k))*ll)
143 a1fminus=a1fminus+ccminus(k)*ximinus(1,k)
144 + *cmplx(tmp*cos(imagpart(lambdaminus(k))*ll),
145 + tmp*sin(imagpart(lambdaminus(k))*ll))
146 a1bminus=a1bminus+ccminus(k)*ximinus(2,k)
147 a2fminus=a2fminus+ccminus(k)*ximinus(3,k)
148 + *cmplx(tmp*cos(imagpart(lambdaminus(k))*ll),
149 + tmp*sin(imagpart(lambdaminus(k))*ll))
150 40 continue
151 a1famplplus=cabs(a1fplus/a1iplus)**2
152 a1bamplplus=cabs(a1bplus/a1iplus)**2
153 a2famplplus=cabs(a2fplus/a1iplus)**2
154 a1famplminus=cabs(a1fminus/a1iminus)**2
155 a1bamplminus=cabs(a1bminus/a1iminus)**2
156 a2famplminus=cabs(a2fminus/a1iminus)**2
157
158 c *********************************************************************
159 c Finally calculate the envelopes of forward and backward traveling
160 c components of the signal (a1f,a1b) and idler (a2f) fields at spatial
161 c location given by z (with 0 < z < L).
162 c *********************************************************************
163 if (z.ne.ll) then
164 a1fplus=cmplx(0.0d0,0.0d0)
165 a1bplus=cmplx(0.0d0,0.0d0)
166 a2fplus=cmplx(0.0d0,0.0d0)
167 a1fminus=cmplx(0.0d0,0.0d0)
168 a1bminus=cmplx(0.0d0,0.0d0)
169 a2fminus=cmplx(0.0d0,0.0d0)
170 do 50,k=1,3
171 a1fplus=a1fplus+ccplus(k)*xiplus(1,k)
172 + *cexp(lambdaplus(k)*cmplx(z,0.0))
173 a1bplus=a1bplus+ccplus(k)*xiplus(2,k)
174 + *cexp(lambdaplus(k)*cmplx(z,0.0))
175 a2fplus=a2fplus+ccplus(k)*xiplus(3,k)
176 + *cexp(lambdaplus(k)*cmplx(z,0.0))
177 a1fminus=a1fminus+ccminus(k)*ximinus(1,k)
178 + *cexp(lambdaminus(k)*cmplx(z,0.0))
179 a1bminus=a1bminus+ccminus(k)*ximinus(2,k)
180 + *cexp(lambdaminus(k)*cmplx(z,0.0))
181 a2fminus=a2fminus+ccminus(k)*ximinus(3,k)
182 + *cexp(lambdaminus(k)*cmplx(z,0.0))
183 50 continue
184 endif
185
186 return
187 end
188
Generated by ::viewsrc::