fortran转换 matlab代码,大家帮帮忙！fortran转换成matlab

该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

Fc

c****************************************************************

c OUTPUT POWER SIMULATION FOR A STRAIGHT STRIPE SLED

c****************************************************************

DOUBLE PRECISION Nto,Ao,c,h,Nold,delta,qel,PI,n2,n1

&,v,Ncar,df,dr,gain,Psp,LHS,Qsp,Po,Qo,PN,QN,RHS

&,Pow,Imax,ISTEP,Iinj,Br,Conf,Clad,Scat,Jinj,lambda

&,Ro,RL,Pin,W,L,d,PL,QZ(500,500),PZ(500,500),z,DZ

&,gmth,gth,Nth,Ith

INTEGER itest,count,Icount,Insteps,IZ,countz

c---------INPUT DATA ----------------------------------------------

OPEN(51,FILE='in4.idat')

read(51,*) W

read(51,*) L

read(51,*) d

read(51,*) Ro

read(51,*) RL

read(51,*) Pin

read(51,*) Imax

read(51,*) Insteps

read(51,*) lambda

read(51,*) Conf

read(51,*) Br

read(51,*) Clad

read(51,*) Scat

read(51,*) Nto

read(51,*) delta

read(51,*) IZ

close(51)

c------------------------------------------------------------------

c... All lengths in um and current in A:

Ao=1.5d-08

c=3.0d+14

n2=3.38d0

n1=3.55d0

h=6.626d-34

qel=1.602d-19

PI=4.0d0*datan(1.0d0)

df=dacos(n2/n1)*(W/L)/(2.0d0*PI)

dr=df

v=c/n1

ISTEP=Imax/dble(Insteps)

DZ=L/dble(IZ)

c------------------------------------------------------------------

c...Threshold modal gain (um-1):

gmth=-dlog(Ro*RL)/(2.0d0*L)

c...Material gain at threshold (um-1):

gth=(gmth-(1.0d0-Conf)*Clad-Scat)/Conf

c...Carrier denisty at threshold (um-3):

Nth=Nto+gth/Ao

c...Threshold current(printed in mA):

Ith=W*L*qel*d*Br*Nth*Nth

print*,Ith*1.0d+03

c------------------------------------------------------------------

OPEN(30,FILE='licurve.dat')

OPEN(32,FILE='pz.dat')

c------------------------------------------------------------------

Nold=1.0d04

DO 100 Icount=1,Insteps

Iinj=dble(Icount)*ISTEP

Jinj=Iinj/(W*L)

LHS=Jinj/(qel*d*v)

call stripe(Conf,Ao,Nto,Br,Pin,Clad,Scat

&,df,dr,Ro,RL,v,L,LHS,delta,Ncar,Psp,Qsp,gain,count

&,Po,Qo)

PL=dexp(gain*L)*Pin/(1.0d0-Ro*RL*dexp(2.0d0*gain*L))

&+(dexp(gain*L)-1.0d0)*(Psp+Ro*dexp(gain*L)*Qsp)

&/(1.0d0-Ro*RL*dexp(2.0d0*gain*L))

Pow=(1.0d0-RL)*PL*v*W*d*h*c*1.0d+03/lambda

c..note Power Output in mW

WRITE(30,9991) Iinj,Pow,count

print *,Iinj,count,Ncar,dsqrt(Jinj/(qel*d*Br))

c &(1.0d0-Ro

c &*RL*dexp(2.0d0*gain*L)),gain

do 50 countz=1,IZ

z=dble(countz)*DZ

PZ(Icount,countz)=Po*dexp(gain*z)+Psp*(dexp(gain*z)-1.0d0)

QZ(Icount,countz)=Qo*dexp(-gain*z)+Qsp*(dexp(-gain*z)-1.0d0)

write(32,9992) z,PZ(Icount,countz),QZ(Icount,countz)

50 continue

write(32,*)

100 CONTINUE

9991 FORMAT(E16.5,4x,E16.5,4x,I8)

9992 FORMAT(E16.5,4x,E16.5,4x,E16.5)

CLOSE(30)

CLOSE(32)

END

c------------------------------------------------------------------------------c

SUBROUTINE stripe(Conf,Ao,Nto,Br,Pin,Clad,Scat

&,df,dr,Ro,RL,v,L,LHS,delta,Ncar,Psp,Qsp,gain,count

&,Po,Qo)

DOUBLE PRECISION Nto,Ao,Nold,delta,v,Ncar,df,dr,gain,LHS

&,Qsp,Po,Qo,PN,QN,RHS,Br,Conf,Clad,Scat,Ro,RL,Pin,L,Psp

INTEGER itest,count

itest=0

count=0

Nold=1.0d+05

do 50 while(itest.eq.0)

Ncar=Nold

gain=Conf*Ao*(Ncar-Nto)+(1.0d0-Conf)*Clad+Scat

Psp=df*Ncar*Ncar*Br/(gain*v)

Qsp=dr*Ncar*Ncar*Br/(gain*v)

Po=(Pin+Ro*(dexp(gain*L)-1.0d0)*(RL*dexp(gain*L)*Psp+Qsp))

&/(1.0d0-Ro*RL*dexp(2.0d0*gain*L))

Qo=RL*dexp(2.0d0*gain*L)*Po

&+RL*dexp(gain*L)*(dexp(gain*L)-1.0d0)*Psp

&+(dexp(gain*L)-1.0d0)*Qsp

PN=(dexp(gain*L)-1.0d0)*(Po+Psp)/(gain*L)-Psp

QN=(1.0d0-dexp(-gain*L))*(Qo+Qsp)/(gain*L)-Qsp

RHS=gain*(PN+QN)+Br*Ncar*Ncar/v

if(dabs(LHS-RHS).lt.delta)then

itest=1

Nold=Ncar

else

itest=0

if(LHS.lt.RHS)then

Nold=Ncar-Ncar*0.001d0

elseif(LHS.gt.RHS)then

Nold=Ncar+Ncar*0.001d0

endif

endif

count=count+1

50 continue

RETURN

END

c------------------------------------------------------------------------------c