c...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+ c File: mopagrat.f (Fortran 77 source for the mopagrat routine) c Author: Fredrik Jonsson c Date: December 26, 2005 c Last change: January 3, 2006 c c Description: The mopagrat routine calculates magneto-optically c parametrically amplified field envelopes of forward c and backward travelling wave components inside a Bragg c grating of sinusoidal spatial distribution of the c refractive index. c c Required external routines: c c ccubsolv - For solving the characteristic complex-valued cubic c polynomial equation that appear for the eigenvalues c of the system determinant of the wave propagation. c c zgaussj - For solving complex-valued linear systems of equations c by Gauss-Jordan elimination. c c Copyright (c) 2006, Fredrik Jonsson c...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+ subroutine mopagrat(alpha,beta,gamma,delta,eta,kappa,ll,z, + a1iplus,a2iplus,a1iminus,a2iminus, + a1famplplus,a1bamplplus,a2famplplus, + a1famplminus,a1bamplminus,a2famplminus, + a1fplus,a1bplus,a2fplus,a1fminus,a1bminus,a2fminus) real*8 alpha,beta,gamma,delta,eta,kappa,ll,z complex*8 a1iplus,a2iplus,a1iminus,a2iminus real*8 a1famplplus,a1bamplplus,a2famplplus real*8 a1famplminus,a1bamplminus,a2famplminus complex*8 a1fplus,a1bplus,a2fplus,a1fminus,a1bminus,a2fminus,c(3) complex*8 bbplus(3,3),xiplus(3,3),lambdaplus(3),ccplus(3) complex*8 bbminus(3,3),ximinus(3,3),lambdaminus(3),ccminus(3) real*8 tmp integer*4 k external ccubsolv,zgaussj c ********************************************************************* c Calculate the eigenvalues of the system matrix A of the magneto- c optical parametric interaction. This is done by solving the cubic c polynomial equation det(A-I*lambda)=0, where I=diag(1,...,1) denotes c the identity matrix, or equivalently of the form c c lambda^3+c[2]*lambda^2+c[1]*lambda+c[0]=0, c c by using the external ccubsolv() routine. Here the coefficients are c as follows: c c c(1)=c[0], c(2)=c[1], c(3)=c[2]. c c After the call, the three complex-valued eigenvalues are returned c in the vector lambda(1..3). c ********************************************************************* c(1)=cmplx(0.0,((abs(kappa)**2)-((alpha+delta)**2)) + *(alpha+delta+beta+gamma)-((eta)**2)*(alpha+delta)) c(2)=cmplx(((alpha+delta)**2)-((kappa)**2)-((eta)**2),0.0) c(3)=cmplx(0.0,-(alpha+delta+beta+gamma)) call ccubsolv(c,lambdaplus) c(1)=cmplx(0.0,((abs(kappa)**2)-((alpha-delta)**2)) + *(alpha-delta+beta-gamma)-((eta)**2)*(alpha-delta)) c(2)=cmplx(((alpha-delta)**2)-((kappa)**2)-((eta)**2),0.0) c(3)=cmplx(0.0,-(alpha-delta+beta-gamma)) call ccubsolv(c,lambdaminus) c ********************************************************************* c Having calculated the eigenvalues of the system matrix A, now c proceed with calculating the corresponding eigenvectors. The eigen- c vectors are stored in the matrix xi, in which the (j,k):th element c is the j:th element of the k:th eigenvector. c ********************************************************************* do 10,k=1,3 xiplus(1,k)=cmplx(1.0,0.0) xiplus(2,k)=cmplx(kappa,0.0) xiplus(2,k)=xiplus(2,k)/(cmplx(0.0,alpha+delta) + +lambdaplus(k)) xiplus(3,k)=(cmplx(0.0,alpha+delta)-lambdaplus(k)) xiplus(3,k)=xiplus(3,k)+cmplx((kappa**2),0.0) + /(cmplx(0.0,alpha+delta)+lambdaplus(k)) xiplus(3,k)=xiplus(3,k)*cmplx(0.0,1.0/eta) ximinus(1,k)=cmplx(1.0,0.0) ximinus(2,k)=cmplx(kappa,0.0) ximinus(2,k)=ximinus(2,k)/(cmplx(0.0,alpha-delta) + +lambdaminus(k)) ximinus(3,k)=(cmplx(0.0,alpha-delta)-lambdaminus(k)) ximinus(3,k)=ximinus(3,k)+cmplx((kappa**2),0.0) + /(cmplx(0.0,alpha-delta)+lambdaminus(k)) ximinus(3,k)=ximinus(3,k)*cmplx(0.0,1.0/eta) 10 continue c ********************************************************************* c Calculate the inverse of the matrix bb by using the Gauss--Jordan c elimination routine zgaussj(), as customized from the Numerical c Recipes' real-valued version of the same algorithm. c ********************************************************************* tmp=max(realpart(lambdaplus(1)), + realpart(lambdaplus(2)),realpart(lambdaplus(3))) do 20,k=1,3 bbplus(1,k)=xiplus(1,k) bbplus(2,k)=xiplus(2,k)* + cexp((lambdaplus(k)-cmplx(tmp,0.0))*cmplx(ll,0.0)) bbplus(3,k)=xiplus(3,k) 20 continue ccplus(1)=a1iplus ccplus(2)=cmplx(0.0,0.0) ccplus(3)=a2iplus call zgaussj(bbplus,3,3,ccplus,1,1) tmp=max(realpart(lambdaminus(1)), + realpart(lambdaminus(2)),realpart(lambdaminus(3))) do 30,k=1,3 bbminus(1,k)=ximinus(1,k) bbminus(2,k)=ximinus(2,k)* + cexp((lambdaminus(k)-cmplx(tmp,0.0))*cmplx(ll,0.0)) bbminus(3,k)=ximinus(3,k) 30 continue ccminus(1)=a1iminus ccminus(2)=cmplx(0.0,0.0) ccminus(3)=a2iminus call zgaussj(bbminus,3,3,ccminus,1,1) c ********************************************************************* c Calculate the logarithmic amplification of the transmitted and c reflected signal and idler, in measures of dB. c ********************************************************************* a1fplus=cmplx(0.0,0.0) a1bplus=cmplx(0.0,0.0) a2fplus=cmplx(0.0,0.0) a1fminus=cmplx(0.0,0.0) a1bminus=cmplx(0.0,0.0) a2fminus=cmplx(0.0,0.0) do 40,k=1,3 tmp=exp(realpart(lambdaplus(k))*ll) a1fplus=a1fplus+ccplus(k)*xiplus(1,k) + *cmplx(tmp*cos(imagpart(lambdaplus(k))*ll), + tmp*sin(imagpart(lambdaplus(k))*ll)) a1bplus=a1bplus+ccplus(k)*xiplus(2,k) a2fplus=a2fplus+ccplus(k)*xiplus(3,k) + *cmplx(tmp*cos(imagpart(lambdaplus(k))*ll), + tmp*sin(imagpart(lambdaplus(k))*ll)) tmp=exp(realpart(lambdaminus(k))*ll) a1fminus=a1fminus+ccminus(k)*ximinus(1,k) + *cmplx(tmp*cos(imagpart(lambdaminus(k))*ll), + tmp*sin(imagpart(lambdaminus(k))*ll)) a1bminus=a1bminus+ccminus(k)*ximinus(2,k) a2fminus=a2fminus+ccminus(k)*ximinus(3,k) + *cmplx(tmp*cos(imagpart(lambdaminus(k))*ll), + tmp*sin(imagpart(lambdaminus(k))*ll)) 40 continue a1famplplus=cabs(a1fplus/a1iplus)**2 a1bamplplus=cabs(a1bplus/a1iplus)**2 a2famplplus=cabs(a2fplus/a1iplus)**2 a1famplminus=cabs(a1fminus/a1iminus)**2 a1bamplminus=cabs(a1bminus/a1iminus)**2 a2famplminus=cabs(a2fminus/a1iminus)**2 c ********************************************************************* c Finally calculate the envelopes of forward and backward traveling c components of the signal (a1f,a1b) and idler (a2f) fields at spatial c location given by z (with 0 < z < L). c ********************************************************************* if (z.ne.ll) then a1fplus=cmplx(0.0d0,0.0d0) a1bplus=cmplx(0.0d0,0.0d0) a2fplus=cmplx(0.0d0,0.0d0) a1fminus=cmplx(0.0d0,0.0d0) a1bminus=cmplx(0.0d0,0.0d0) a2fminus=cmplx(0.0d0,0.0d0) do 50,k=1,3 a1fplus=a1fplus+ccplus(k)*xiplus(1,k) + *cexp(lambdaplus(k)*cmplx(z,0.0)) a1bplus=a1bplus+ccplus(k)*xiplus(2,k) + *cexp(lambdaplus(k)*cmplx(z,0.0)) a2fplus=a2fplus+ccplus(k)*xiplus(3,k) + *cexp(lambdaplus(k)*cmplx(z,0.0)) a1fminus=a1fminus+ccminus(k)*ximinus(1,k) + *cexp(lambdaminus(k)*cmplx(z,0.0)) a1bminus=a1bminus+ccminus(k)*ximinus(2,k) + *cexp(lambdaminus(k)*cmplx(z,0.0)) a2fminus=a2fminus+ccminus(k)*ximinus(3,k) + *cexp(lambdaminus(k)*cmplx(z,0.0)) 50 continue endif return end