На языке FORTRAN - IV

 

program SUPER

c с учетом неоднородного уширения в единицах W^-1

c использован исправленный метод Эйлера

c единица измерения омега

REAL*8 AIZ(400,6),SPECTR(400),SPECTR0(400),ZL(200,6),st(400),

*DOP(60),DZ(60,5),Z(800,60),Z1(60),Z2(60),tt,am(100),INT(100),

*Z_(800,60)

COMPLEX*8 R(800,60),A(800),aa(800),r22(800,60),E1(800),

*R_(800,60),A_(800),A0(800),A1(800),

*R1(60),SPEC(400),SPEC1(400),RS(60),R2(60),E(800)

COMPLEX*8 A2,RK,RK1,W1,AT2,DJK

REAL*8 DD,GD,DGD,B,TP,SGM,ARG,AMP,nmax

open(2,file='SPECTR0.dat',status='unknown')

open(3,file='dopler.dat',status='unknown')

open(4,file='st.dat',status='unknown')

open(7,file='koger.dat',status='unknown')

open(8,file='Z.dat',status='unknown')

open(9,file='IN.dat',status='unknown')

OPEN(10,FILE='A0.DAT',STATUS='UNKNOWN')

OPEN(11,FILE='A1.DAT',STATUS='UNKNOWN')

open(12,file='spectr.dat',status='unknown')

open(13,file='neod1.dat',status='unknown')

open(14,file='neod2.dat',status='unknown')

OPEN(1,FILE='somega',STATUS='OLD')

 

READ(1,*) T,AL,DT,DS,GD,g2,TP,B,m,ND,t12,t23, Z0,X0,Y0

 

c входные параметры

c T- ВРЕМЯ ПРОЦЕССА, (20-500)

c AL- ДЛИНА ОБРАЗЦА, (1-20)

c DT - ШАГ ИНТЕГРИРОВАНИЯ ПО ВРЕМЕНИ И ДЛИНЕ, (0.005-0.01)

c DS - ШАГ ИНТЕГРИРОВАНИЯ ДЛЯ ВЫЧИСЛЕНИЯ ФУРЬЕ

c СПЕКТРА

c GD- ДИСПЕРСИЯ НЕОДНОРОДНОГО КОНТУРА,

c G2 - ВЕЛИЧИНА ОДНОРОДНОГО УШИРЕНИЯ,

c TP - ДЛИТЕЛЬНОСТЬ ВХОДНОГО

c ИМПУЛЬСА НА ПОЛОВИНЕ ВЫСОТЫ,

c 2*B - КОЭФФИЦИЕНТ ЗАДАЮЩАЯ ПЛОЩАДЬ ИМПУЛЬСА

c M - ПАРАМЕТР РАВНЫЙ ЕДИНИЦЕ,

C ND - ЧИСЛО УЗЛОВ ПО НЕОДНОРОДНОМУ КОНТУРУ, при

c отсутствии неоднородного контура ND=1

C T12,T23 - ВРЕМЯ ЗАДЕРЖКИ МЕЖДУ ПЕРВЫМ И

с ВТОРЫМ ВОЗБУЖДАЮЩИМИ ИМПУЛЬСАМИ

C Z0 - НАЧАЛЬНОЕ ЗНАЧЕНИЕ ИНВЕРСИИ АТОМНОЙ СИСТЕМЫ

C X0-НАЧАЛЬНОЕ ЗНАЧЕНИЕ ПОЛЯРИЗАЦИИ АТОМНОЙ СИСТЕМЫ

C Y0 - НАЧАЛЬНОЕ ЗНАЧЕНИЕ ПОЛЯРИЗАЦИИ АТОМНОЙ

с СИСТЕМЫ (R0=X0+i*Y0)

 

PI=3.1415926

n1=20./dt

n2=32./dt

n3=44./dt

n4=80./dt

T0=TP*2.5

dd=2.5*gd/nd

W1=CMPLX(0.0,1.0)

C=0.0

amax=0.0

m1=m

 

WRITE(6,2)T,AL,ND,DT,dd,DS,G2,GD,T0,TP,B,M,t2

2 FORMAT (' T=',g12.5/,' AL=',g12.5/,' ND=',I5/,' DT=',g12.5/,

*' DD=',g12.5/,' DS=',g12.5/,' G2=',g12.5/,' GD=',E12.5/,

*' T0=',E12.5/,' TP=',E12.5/,' B=',E12.5/,' M=',I5/

*' t2=',g10.5)

NL=AL/DT

NT=T/DT

IB=NT/200

NB=NT/IB

NLL=NL-1

J0=T0/DT

DO 3 K=1,2*ND-1

DGD=-(((K-ND)*DD)/GD)**2/2

DOP(K)=EXP(DGD)/(GD*SQRT(2.*PI))

IF (ND.EQ.1) DOP(K)=1.0

3 CONTINUE

if(ND.EQ.1) DD=1.0

c DO 55 LL=1,1

C задание начальных значений поля и поляризации

DO 4 I=1,NL

DO 5 K=1,2*ND

Z(I,K)=Z0

5 R(I,K)=CMPLX(X0,Y0)

st(i)=0.0

4 A(I)=CMPLX(0.0,0.0)

DO 52 I=1,NB

DO 51 K=1,6

51 AIZ(I,K)=0.0

52 CONTINUE

JGR=1

L=1

DO 6 J=1,400

SPEC1(J)=CMPLX(0.,0.)

6 SPEC(J)=CMPLX(0.,0.)

DO 62 I=1,400

DO 62 K=1,6

61 ZL(I,K)=0.0

62 CONTINUE

D2=1.-DT*G2

SE=0.0

SI=0.0

s0=0.0

S=0.0

SEIN=0.0

Z0=0.0

DC=1.-DT*C

NLL=NL-1

C вычисление

SGM=TP/2.355

SGM1=SGM

AMP=B*PI/(SQRT(2*PI)*sgm)

AMP1=AMP

AMAX=0.0

DO 7 J=1,NT

C A(1)=CMPLX(0.,0.)

ZT=0.0

J1=J0*2

j2=t12/dt

j3=j2+3*j0

j4=j2+4*j0

ARG=0.0

с задание граничных значений для поля

c arg=b*2.64/(tp*cosh(2.64*dt*(j-j0)/tp))

c a(1)=cmplx(arg,0.0)

IF(J.LE.J1) ARG=DT*(J-J0)/SGM

IF(J.LE.J1) A(1)=cmplx(AMP*EXP(-(ARG*ARG)/2),0.0)

C IF(J.GE.J1) A(1)=cmplx(AMP*EXP(-0.5*dt*(j-j3)*dt*(j-j3)

C */(sgm*sgm)),0.0)

C if(j.GE.j4) a(1)=cmplx(AMP1*EXP(-0.5*DT*(J-(J4+j1))*DT*(J-(J4+j1))

C * *SGM1*SGM1),0.0)

DO 8 K=1,2*ND-1

8 R1(K)=R(1,K)

DO 9 I=1,NLL

e1(i+1)=a(i+1)

A0(I)=CMPLX(0.,0.)

DO 10 K=1,2*ND-1

DK=DD*(K-ND)

RS(K)=R1(K)

R1(K)=R(I+1,K)

Z_(I+1,K)=Z(I+1,K)-DT*(REAL(R1(K)*CONJG(A(I+1))

*+A(I+1)*CONJG(R1(K))))

R_(I+1,K)=R1(K)+DT*(2*Z(I+1,K)*A(I+1)+W1*DK*R1(K))

 

10 A0(I)=A0(I)+RS(K)*DOP(K)*DD

A_(I+1)=A(I)+DT*A0(I)

9 CONTINUE

C ОДНОШАГОВАЯ КОРРЕКЦИЯ

DO 99 I=1,NLL

A1(I)=CMPLX(0.,0.)

A2=A_(I+1)

aa(i)=a2

ZLL=0.0

DO 11 K=1,2*ND-1

Z2(K)=Z_(I+1,K)

R2(K)=R_(I+1,K)

DKT=DD*(K-ND)

RR=R(I+1,K)+DT*(Z(I+1,K)*A(I+1)+Z2(K)*A2)+0.5*(W1*DKT*

*(R(I+1,K)+R2(K))-G2*(R(I+1,K)+R2(K)))

Z(I+1,K)=Z(I+1,K)-0.5*DT*REAL((R(I+1,K)*CONJG(A(I+1))+

*CONJG(R(I+1,K))*A(I+1)+R2(K)*CONJG(A2)+A2*CONJG(R2(K))))

R(I+1,K)=RR

IF(J.EQ.NT) GO TO 101

GO TO 102

101 IF (K.EQ.1) ZL(I+1,1)=Z(I+1,1)

IF (K.EQ.3) ZL(I+1,2)=Z(I+1,3)

IF (K.EQ.5) ZL(I+1,3)=Z(I+1,5)

IF (K.EQ.7) ZL(I+1,4)=Z(I+1,7)

IF (K.EQ.9) ZL(I+1,5)=Z(I+1,9)

IF (K.EQ.10) ZL(I+1,6)=Z(I+1,10)

c write(3,*) i*dt,zl(i+1,5),ZL(I+1,10),ZL(I+1,15),ZL(I+1,20)

ZLL=ZLL+DD*DOP(K)*(Z(I+1,K)+0.5)

102 CONTINUE

11 A1(I+1)=A1(I+1)+R_(I,K)*DOP(K)*DD

ZAL=ZLL

if (j.eq.nt) ZT=ZT+ZAL/NLL

st(i)=st(i)+dt*real(a(i))

99 E(I+1)=0.5*(A(I)+aa(i)+dt*(A0(I)+A1(I)))

DO 12 I=1,NLL

12 A(I+1)=E(I+1)

s0=s0+dt*cabs(a(1))**2/al

SI=SI+DT*CABS(A(NL))**2/al

SE=SE+2*DT*REAL(A(NL))

SEIN=SEIN+2*DT*REAL(A(1))

C ВЫЧИСЛЕНИЕ ГРАНИЧНЫХ ЗНАЧЕНИЙ

DO 13 K=1,2*ND-1

DKT=DD*(K-ND)

RK=R(1,K)

ZK=Z(1,K)

RK1=RK*D2+DT*(A(1)*ZK+W1*DKT*RK)

ZK1=ZK-DT*(REAL(CONJG(A(1))*RK+CONJG(RK)*A(1)))

R(1,K)=0.5*(RK+RK1)+0.5*DT*(W1*DKT*(RK1+RK)+2*ZK1*A(1))

Z(1,K)=0.5*(ZK+ZK1)-DT*0.5*(REAL(CONJG(A(1))*

*RK1+A(1)*CONJG(RK1)))

с вычисление ошибок ограничения

C f1=0.2*real(r(nll,nd)-r2(nd))

C f2=0.2*real(z(nll,nd)-z2(nd))

C f3=0.2*real(a(nll)-aa(nll))

C if(f1.ne.0.0) f11=f1

C if(f2.ne.0.0) f22=f2

C if(f3.ne.0.0) f33=f3

с формирование выходных файлов с расширением dat

if(j.eq.n1) write(3,*) dkt,(z(nll,k)+0.5)*dop(k)

c if(j.eq.nt) write(13,*) dkt,(z(nll,k)+0.5)*dop(k)

c if(j.eq.n3) write(14,*) dkt,(z(nll,k)+0.5)*dop(k)

c if(j.eq.n4) write(4,*) dkt,(z(nll,k)+0.5)*dop(k)

13 CONTINUE

C spectr

write(6,40) j,z(nll,nd)+0.5

C ВЫЧИСЛЕНИЕ ФУРЬЕ СПЕКТРА

DO 14 K=1,400

DJK=CMPLX(0.0,DT*J*DS*(K-200))

SPEC1(K)=SPEC1(K)+DT*A(1)*CEXP(DJK)

14 SPEC(K)=SPEC(K)+DT*A(NL)*CEXP(DJK)

C ФОРМИРОВАНИЕ МАССИВОВ ДЛЯ ВЫДАЧИ

JJ=MOD(J+IB,IB)

IF(JJ.EQ.0) GO TO 15

GO TO 7

15 AIZ(L,1)=L

AIZ(L,2)=REAL(A(1))

AIZ(L,3)=REAL(A(NLL))

AIZ(L,5)=cabs(A(NLL))**2/al

AIZ(L,4)=REAL(2*Z(NLL,ND)+1)*0.5

AIZ(L,6)=cabs(A(1))**2/al

c if(aiz(l,5).ge.amax) amax=aiz(l,5)

c if(aiz(l,5).lt.amax) amax=amax

L=L+1

7 CONTINUE

c AM(LL)=AMAX

 

c вычисление спектра излучения

DO 16 K=1,400

ts=(k-200)*ds

SPECTR(K)=CABS(SPEC(K))**2/(2*PI*al)

SPECTR0(K)=CABS(SPEC1(K))**2/(2*PI*al)

write(12,*) ts,spectr(k)

WRITE(2,*) TS,SPECTR0(K)

16 S=S+DS*SPECTR(K)/al

c WRITE(6,20) SI,SEIN,SE,ZLL,ZT,S

zc0=s0+1.

zct=si+zt

WRITE(7,20) zc0,zct,SEIN,SE,ZLL,ZT,S,g2,gd,b,al,M,tp,t0,

* nd,dd,ds

С ФОРМИРОВАНИЕ ДАТОВСКИХ ФАЙЛОВ ДЛЯ Пакетов

c GRAFER и ORIGIN

do 30 L=1,NB

tt=L*DT*IB

WRITE(8,*) tt,AIZ(L,4)

write(9,*) tt,AIZ(L,5)

WRITE(10,*) tt,AIZ(L,2)

WRITE(11,*) tt,AIZ(L,3)

c WRITE(4,*) tt,AIZ(L,6)

30 continue

С ФОРМИРОВАНИЕ ДАТОВСКИХ ФАЙЛОВ ДЛЯ ПРОГРАММЫ

c MICROCAL ORIGIN

 

do 34 i=1,nl

write(14,*) (z(i,k),k=1,2*nd)

34 continue

20 FORMAT (' zc0=',G10.5,

*' zct=',G10.5,

*' SEIN=',G10.4,

*' SE=',G11.5//,2x,

*' Z0=',G10.4,

* ' ZT=',G10.4,

* ' S=',G10.4/,

*' g2=',g10.4,

* ' gd=',g10.4,

* ' teta=',g10.4,

*' l=',g10.5/,

*' M1=',I4,

* ' tp=',g10.4,

* ' t0=',g10.4/,

*' nd=',i4,

*' dd=',g10.4,

*' ds='g10.4,

*' t2=',g8.2)

32 FORMAT (E5.3,4X,E12.5)

31 FORMAT (E5.3,4X,G10.5)

40 FORMAT (4X,i4,3E12.5)

STOP

END

 

Пример входных параметров для задачи когерентного усиления УКИ cвета (задается в файле SOMEGA).

20 7 0.01 0.05 0.5 0.05 1.6 0.25 1 30 1.0 2.0 0.5 0.0 0.0