Debugging

[SydFard]

Could you please help me debug this FORTRAN 77 code.
Thank you for your time.

C Grating.for - Grating with lenses for collimation and Focusing
C Derived based on the original Matlab file: mylensGR_200_asym_ausome_VER5.m
C
c Include library for random #s and time
INCLUDE 'FSUBLIB.FI'
INTEGER*2 HRS, MINS, SECS, HSECS
INTEGER NLP2
REAL
1 Ein(1045), yin(1045), yout(1045),
2 xFP(-47:47,-5:5), yFP(-47:47,-5:5),
3 K1, K2, K3, FW, FPW, dyin,
4 powG(5), powFL(5), powOUT(5), ambda(5),
5 y0, dyLin, xLB, xGB, xin, Phase,
6 Ndex1, Ndex2, Ndex3, xL(1045), xFL(1045),
7 yLB(1045), yL(1045), yGB(1045),
8 yFL(1045)
COMPLEX EL(-2:2,-522:522),
1 ELB(-2:2,-522:522), EGB(-2:2,-522:522),
2 EG(-2:2,-522:522), EFL(-2:2,-522:522),
3 Eout(-2:2,-522:522), J, JK1, JK2,
4 JK12x, JK3, JK32x, JK22x,
5 dxr, dxh
C
C Open files ...
OPEN (UNIT=4,file='Grating2.dat')
OPEN (UNIT=7,file='Grating.out')
OPEN (UNIT=8,file='Gscrn.out')
C Get clock time, calculate random seed
CALL GETTIM( hrs, mins, secs, hsecs )
ISEED = hrs*mins*secs*hsecs
Write(6,*) ('.............................................')
Write(8,*) ('.............................................')
Write(6,19) hrs,mins,secs
Write(8,19) hrs,mins,secs
19 Format('.............Run started at ',I2,':',I2,':',I2)
C
C Calculate useful constants
J = (0.,1.0)
Pi = 3.1415927
TwoPi = 2.*Pi
Piby2 = Pi/2.
C
C Read data
READ(4,93) f, d, R, R2, a, Ndex1, Ndex2, Ndex3, ambda0
READ(4,93) ambda1, dlambda, S, gwmin, ms, xGB
READ(4,91) n, NF, NFpts, NOP, NINP, inrange, ioutrange, NWL, NLP
C
C Calculate wavelengths
NWL = 2*((ambda0-ambda1)/dlambda)+1
Write(6,*) '# of WLs = ', NWL
Write(8,*) '# of WLs = ', NWL
Do 10 iwl=1,NWL
ambda(iwl) = ambda1+(iwl-1)*dlambda
10 Continue
Write(6,*) 'ambda = ', (ambda(iwl),iwl=1,NWL)
Write(8,*) 'ambda = ', (ambda(iwl),iwl=1,NWL)
Write(7,91) NWL
Write(6,*) ' '
Write(8,*) ' '

Write(6,*) ('.............................................')
Write(8,*) ('.............................................')
Write(6,11) ambda0
Write(8,11) ambda0
11 Format('.............ambda0 is ',F10.0)
CALL paused

C
C Generate normalized Gaussian mode field
C make sure NINP is odd
if (MOD(NINP,2).EQ.0) then
NINP = NINP+1
endif
NINP2 = (NINP-1)/2
dyin = ms*inrange/(NINP-1)
y0=0.
Do 20 inp=-NINP2,NINP2
yin(inp) = inp*dyin
Ein(inp) = exp(-(yin(inp)-y0)**2/(ms/2)**2)
20 Continue
xin = xI
C output plane
NOP2 = (NOP-1)/2
dout = ioutrange*S/(NOP-1)
Do 30 ioutp=-NOP2,NOP2
yout(ioutp) = ioutp*dout
xout = xD
30 Continue
C
C Calculate grating order, facet width,
n = ((ambda0/dlambda)/(NWL+1))
FW = (yup-ylow)/NF
Write(6,*) 'Grating Order n = ', n
Write(6,*) ' Facet Width = ', FW
Write(8,*) 'Grating Order n = ', n
Write(8,*) ' Facet Width = ', FW
C
C Display total Grating width, number of facets, Total # of points on lens and Total # of facet points
TotGW = yup-ylow
If((NF/2)*2.eq.NF) then NF=NF+1 !Make sure NF is odd
TFP = NF*NFpts
NLP = NF*NFpts
Write(6,51) TotGW
Write(8,51) TotGW
51 Format('Total grating width = ',F6.1,' microns')
Write(6,52) NF
Write(8,52) NF
52 Format('Total # of facets = ',I6)
Write(6,53) TFP
Write(8,53) TFP
53 Format('Total # of facet points = ',I6)
Write(6,54) NLP
Write(8,54) NLP
54 Format('Total # of lens points = ',I6)

C
If(TFP.gt.2**16) then
Write(6,*) 'Too many total facet points: Stop'
Write(8,*) 'Too many total facet points: Stop'
Stop
endif
C
C displacement from facet to facet
dx = n*ambda0/(Ndex1-Ndex2)
C
C Input output centres
xI = -d-3*R
yI = 0.0
xD = 8.3*R
yD = yI
C
C search lens boundary
yup = Sqrt(R**2/2)
ylow = -Sqrt(R**2/2)
xup = -d-(R-Sqrt(R**2/2))
xlow = xup
f2= f
C
C parabolic collimating lens, lens boundary, parabolic focusing lens
dyLin = (yup-ylow)/(NLP-1)
NLP2 = (NLP-1)/2
Do 60 iy=-NLP2,NLP2
yL(iy) = iy*dyLin
xL(iy) = -d-yL(iy)**2/(4*f)
yLB(iy) = yL(iy)
xFL(iy) = -d-yL(iy)**2/(4*f2)+d+2*R
yFL(iy) = yL(iy)
yGB(iy) = yLB(iy)
60 Continue
xLB = -10*R

C
C Set up Grating Grid based on central wavelength
xGB = -4E-3
FW = (yup-ylow)/NF
FPW = FW/NFpts
NF2 = (NF-1)/2
NFp2 = (NFpts-1)/2
c loop thru facets
Do 70 mf=-NF2,NF2 !facet # loop
c loop thru pt #
Do 65 mp=-NFp2,NFp2 !facet point loop
xFP(mf,mp) = mf*dx
yFP(mf,mp) = mp*FPW+mf*FPW*(NFpts-1)
65 Continue
70 Continue
C
call paused
C==========================================================
C..Loop thru wavelengths
C
Do 1000 iwl=1,NWL
K1 = 2*pi/ambda(iwl)*Ndex1
K2 = 2*pi/ambda(iwl)*Ndex2
K3 = 2*pi/ambda(iwl)*Ndex3
C----------------------------------------------------------
C Calculate field on collimating lens
sum=0.0
JK1=J*K1
Do 100 n=1,TFP
DelX=xL-xin
DelY=yL-yin
JK12x = JK1/(2.*DelX)
Phase = K1*DelX
DelY2 = DelY**2
dxr = dyin*Csqrt(JK12x/Pi)*Cexp(-J*Phase)
dxh = dxr*cexp(-JK12x*DelY2)
theta_d_in = atan((DelY)/(DelX))
term = dxh*(1+cos(theta_d_in))/2*Ein
sum = sum + term
EL(iwl,n) = sum
100 Continue

C Calculate field on boundary of flat surface behind first lens , using Huygens principle
JK3=J*K3
Do 200 n=1,TFP
DelX=xLB-xL
DelY=yLB-yL
JK32x = JK3/(2.*DelX)
Phase = K3*DelX
DelY2 = DelY**2
dxr = dyLin*Csqrt(JK32x/Pi)*Cexp(-J*Phase)
dxh = dxr*cexp(-JK32x*DelY2)
theta_d_in = atan((DelY)/(DelX))
term = dxh*(1+cos(theta_d_in))/2*EL(iwl,n)
sum = sum + term
ELB(iwl,n) = sum

!theta_d_in = atan((yLB-yL)/(xLB-xL))
!ELB(nWL,n) = sum(sqrt(i*k3/2/pi/(xLB-xL))*Cexp(-i*k3*(xLB-xL))*Cexp(-i*k3*((yLB-yL)**2)/(2*(xLB-xL)))*(1+cos(theta_d_in))/2*EL(nWL,*dyLin)
200 Continue
C Calculate field on grating boundary
JK1=J*K1
Do 300 n=1,TFP
DelX=xGB-xLB
DelY=yGB-yLB
JK12x = JK1/(2.*DelX)
Phase = K1*DelX
DelY2 = DelY**2
dxr = dyLin*Csqrt(JK12x/Pi)*Cexp(-J*Phase)
dxh = dxr*cexp(-JK12x*DelY2)
theta_d_in = atan((DelY)/(DelX))
term = dxh*(1+cos(theta_d_in))/2*ELB(iwl,n)
sum = sum + term
EGB(iwl,n) = sum

!theta_d_in = atan((yGB-yLB)/(xGB-xLB))
!EGB(nWL,n) = sum(sqrt(i*k1/2/pi/(xGB-xLB))*Cexp(-i*k1*(xGB-xLB))*Cexp(-i*k1*((yGB-yLB)**2)/(2*(xGB-xLB)))*(1+cos(theta_d_in))/2*ELB(nWL,*dyLin)
300 Continue

C Calculate field on grating
powGG = 0.0
JK2 = J*K2
Do 500 n=1,TFP
Do 400 mf=-NF2,NF2
Do 400 mp=-NFp2,NFp2
DelX=xFP(mf,mp)-xGB
DelY=yFP(mf,mp)-yGB
JK22x = JK2/(2.*DelX)
Phase = K2*DelX
DelY2 = DelY**2
dxr = FPW*Csqrt(JK22x/Pi)*Cexp(-J*Phase)
dxh = dxr*cexp(-JK22x*DelY2)
theta_d_in = atan((DelY)/(DelX))
term = dxh*(1+cos(theta_d_in))/2*EGB(iwl,n)
sum = sum + term
EG(iwl,n) = sum

!theta_d_in = atan((yFP-yGB)/(xFP-xGB))
!EG(nWL,n) = sum(sqrt(i*k2/2/pi/(xFP-xGB))*Cexp(-i*k2*(xFP-xGB))*Cexp(-i*k2*((yFP-yGB)**2)/(2*(xFP-xGB)))*(1+cos(theta_d_in))/2*EGB(nWL,*FPw)
powGG = powGG+abs(EG(iwl,n))**2*FPW
400 Continue
500 Continue
Do 550 n=1,TFP
EG(iwl,n) = EG(iwl,n)/sqrt(powGG)
550 Continue
powG(iwl) = powGG

C Calculate field on focusing lens
Do 700 n=1,TFP
Do 600 mf=-NF2,NF2
Do 600 mp=-NFp2,NFp2
DelX=xFL-xFP(mf,mp)
DelY=yFL-yFP(mf,mp)
JK12x = JK1/(2.*DelX)
Phase = K1*DelX
DelY2 = DelY**2
dxr = FPW*Csqrt(JK12x/Pi)*Cexp(-J*Phase)
dxh = dxr*cexp(-JK12x*DelY2)
theta_d_in = atan((DelY)/(DelX))
term = dxh*(1+cos(theta_d_in))/2*EG(iwl,n)
sum = sum + term
EFL(iwl,n) = sum


!theta_d_F = atan((yL-yFP)/(xFL-xFP))
!EFL(nWL,n) = sum(sqrt(i*k1/2/pi/(xFL-xFP))*Cexp(-i*k1*(xFL-xFP))*Cexp(-i*k1*((yL-yFP)**2)/(2*(xFL-xFP)))*(1+cos(theta_d_F))/2*EG(nWL,*FPw)
600 Continue
700 Continue
Do 750 n=1,TFP
field_intens = abs(EFL(iwl,n))**2*(yup-ylow)/(NLP-1)
sum = sum + field_intens
powFL(iWL) = sum
750 Continue


C Calculate field on output plane
Do 800 n=1,NOP
DelX=xout-xFL
DelY=yout-yL
JK32x = JK3/(2.*DelX)
Phase = K3*DelX
DelY2 = DelY**2
dxr = dyLin*Csqrt(JK32x/Pi)*Cexp(-J*Phase)
dxh = dxr*cexp(-JK32x*DelY2)
theta_d_in = atan((DelY)/(DelX))
term = dxh*(1+cos(theta_d_in))/2*EFL(iwl,n)
sum = sum + term
Eout(iwl,n) = sum

!theta_d_out = atan((yout-yL)/(xout-xFL))
!Eout(nWL,n) = sum(sqrt(i*k3/2/pi/(xout-xFL))*Cexp(-i*k3*(xout-xFL))*Cexp(-i*k3*((yout-yL)**2)/(2*(xout-xFL)))*(1+cos(theta_d_out))/2*EFL(nWL,*dyLin)
800 Continue

Do 850 n=1,TFP
field_intens = abs(Eout(iwl,n))**2*dout
sum = sum + field_intens
powOUT(iWL) = sum
850 Continue

C
CALL GETTIM( hrs, mins, secs, hsecs )
Write(6,901) iwl, Ein, EG, EFL, Eout,
1 hrs,mins,secs
Write(8,901) iwl, Ein, EG, EFL, Eout,
1 hrs,mins,secs
901 Format(I4,4(3x,f8.5),4x,I2,':',I2,':',I2)
C
C Save input field for this wavelength
WRITE(7, *) yin
WRITE(7,92) (yin(inp),inp=-NINP2,NINP2)
WRITE(7,92) (abs(Ein(inp)**2),inp=-NINP2,NINP2)
C Save output field for this wavelength
WRITE(7, *) yout
WRITE(7,92) ( yout(ioutp), ioutp=-NOP2,NOP2)
WRITE(7,92) (Cabs(Eout(iwl,ioutp)**2),ioutp=-NOP2,NOP2)
C
1000 Continue !End of wavelength loop
C

Write(6,*) ' '
Write(8,*) ' '
CALL GETTIM( hrs, mins, secs, hsecs )
Write(6,99) hrs,mins,secs
Write(8,99) hrs,mins,secs
CLOSE(4)
CLOSE(7)
CLOSE(8)
STOP
91 FORMAT(I4)
92 FORMAT(4(2x,F15.7))
93 FORMAT(F10.0)
95 Format(1x,I3,5x,F6.4,5x,f7.1)
99 Format('Run finished at ',I2,':',I2,':',I2,'.............')
End
C-----------------------------------------------------------------
C Subroutine to pause execution
SUBROUTINE Paused
Write(6,*) '...Paused '
Write(6,*) '...Enter 0 to Abort, other integer to go'
read*, igo
if (igo.EQ.0) then
Write(6,*) ' !!!!!!!ABORTED!!!!!!!'
Write(8,*) ' !!!!!!!ABORTED!!!!!!!'
STOP
ENDIF
RETURN
END
C-----------------------------------------------------------------

 

[SydFard]

I forgot to include the data file "Grating2.dat".

data file:

500 f focal length of the collimating lens
7500 d Distance from first lens vertex to device center
500 R Radius of the collimating lens
12000 R2 Radius of the focusing lens
-2500 a Distance from device center to grating center
1.50 Ndex1 Core index
1.33 Ndex2 Grating index
1.00 Ndex3 Index behind lenses
0.6 ambda0 design wavelength
0.5 ambda1 second wavelength
0.05 dlambda simulation wavelength step
40 S scale unit for output plane
1000 gwmin minimum width of grating area
40 ms input Gaussian field mode size
-4000 xGB grating boundary
0.6 CWL central wavelength
0.05 DWL wavelength spacing
-2 n grating diffraction order
95 NF Number of facets
11 NFpts number of points per facet, odd
1045 NOP number of output field points, odd
1045 NINP number of input points, odd
4 inrange input field calculation range, in number of mode size
4 ioutrange output field calculation range, in number of S
5 NWL
1045 NLP

 

[ArkM]

1. Use code tags to present your sources (see Process TGML link on the form).
2. What's your problem with this program? Regrettably, I (and many others on this forum) am not a wizard and have no telepathic capabilities...

 

[SydFard]

Hi,
My appology, I don't know what I was thinking!
My problem with this code is that it compiles with no errors and warnings but after all it dies out when I run the program. I am trying to output "Eout" against "yout" in a file and since I am a novice in FORTRAN I am not able to do this on my own. I greatly appreciate your input.
Thanks

Seyed

 

[ArkM]

Sorry, can't run this program on my installation: too many external dependencies (INCLUDEs, external functions etc)...

Try to print here, there and everywhere, for example, insert:

      ...
      PRINT *, 'Its me here! NINP = ', NINP
      ...