Fortran Error: Run-time error M6201-MATH -atan2: Domain Error

[logi]

Hi All,



I have some problem about my Fortran program. When I execute my program, it says



Run-time error M6201-MATH



-atan2 : Domain error







I don’t know what to do. Please help me out from this error.







Thanks



Logi

 

[zeitghost]

You haven't linked some maths library or other into your build.

It would be amazingly useful if you'd told us which compiler & platform your FORTRAN program is running (or not running) on.

 

[alphysicist]

Hi logi,

Can you post at least a small program that recreates the error? Sometimes it's much easier to read the program rather than think about error messages.

However, in this case what I think your program is probably doing is calling atan2 with both arguments of zero:

atan(0.,0.)

which is not allowed (at least with my compilers).

If that's not it, maybe you're program is calling atan2 with integers instead of real values, which will also give an error. It's just hard to tell without a program to look at.

 

[logi]

Hi all,
Thank you for the reply. I sending my program( not all because long).
Please check that for me.
Thanks again
Logi



C PROGRAM silocoldfront
REAL*8 GR,GRMX,K,MU,A,B,D1,L,D2,S2,S3,SEMI,ST,EN,KMAX,LMAX,MUMIN,
&CI
COMPLEX*16 CGUESS(100)
INTEGER N,ISZ
COMMON/CB1/ K,MU,A,B,D1,D2,L,S2,S3
C
C OPENING GKS AND OPENING AND ACTIVATING THE WORKSTATION
C A LIST OF DEFAULT PARAMETER VALUES
C THE DEFAULT ALONG-SHELF WAVE NUMBER
C THE MOST UNSTABLE MODE FOR MU=1.5
C K=1
C C
C THE MOST UNSTABLE MODE FOR MU=1.43

K=1.0D0
CC
C CURRENT WIDTH
L =2.5D0
C
C THE DEFAULT CURRENT HALF WIDTH
C (A=1.65 CORRESPONDS TO AN INTERNAL DEFORMATION RADIUS)
C
A=1.65D0

C A1=A-L/2

C A2=A+L/2
C
C THE DEFAULT SHORE DISTANCE
C
B=1.88D0*A
C
D1=0.624D0*A

D2=1.255D0*A

C THE SLOPE2
S2=4.13D0
C THE SLOPE3
S3=0.83D0

C THE DEFAULT INTERACTION PARAMETER
C (MU IS THE RATIO OF THE BAROCLINIC STRETCHING TO THE SLOPE)

C
MU=1.43D0
C
C THE DEFAULT SEMI-CIRCLE RADIUS
C
SEMI=DSQRT(4.0D0*MU/(L*K**2))
C
C THE MAXIMUM WAVENUMBER
C
KMAX=DSQRT(2.0D0*MU/L)+DSQRT(l+4.0D0*MU/L)
C
C THE MAXIMUM CURRENT HALF-WIDTH
C
IF(K.NE.1.0D0)LMAX=16.0D0*MU*(K/(K**2-1.0D0))**2
IF(K.LE.1.0D0)LMAX=614.603175D0
C
C THE MINIMUM INTERACTION PARAMETER
C
MUMIN=L*((K**2-1.0D0)/K)**2/16.0D0
IF(K.LE.1.0D0)MUMIN=0.0D0
C PARAMETER INTERVAL END POINTS
C
ST=0.0D0
EN=1.0D0
THE CURRENT BOUNDARY PERTURBATION
C
EP=CDEXP(CDSQRT(K**2-1.0D0/C)*(B+A-L/2))
EM=CDEXP(-CDSQRT(K**2-1.0D0/C)*(B+A-L/2))
SN=(EP-EM)/2.0D0
EAP=CDEXP(CDSQRT(K**2-1.0D0/C)*(A+L/2))
EAM=CDEXP(-CDSQRT(K**2-1.0D0/C)*(A+L/2))
EDP=CDEXP(CDSQRT(K**2+S2/C)*D1)
EDM=CDEXP(-CDSQRT(K**2+S2/C)*D1)
C
C AT Y = A-L/2
C
PHI1=MU*AL3*SN/(1.0D0-C)
PHI1R=DREAL(PHASE*PHI1)
PHT1=A-L/2+PHI1R
C
C AT Y = A+L/2
C
PHI2=MU*AL4*EAM/(1.0D0-C)+MU*AL5*EAP/(1.0D0-C)
PHI2R=DREAL(PHASE*PHI2)
PHT2=A+L/2+PHI2R

C CROSS-SHELF SECTION
C
DO 1000 I=1,1001
Y=YL+(YR-YL)*DFLOAT(I-1)/1000.D0
ETAY=DCMPLX(1.65D0,1.65D0)
C
C ISOPCYNAL HEIGHT AT Y=A1 AND Y=A2 RESPECTIVELY
C
IF(Y.LT.1.65D0)HYM=8.0D0*PHI1/(A-L/2)
IF(Y.GT.1.65D0)HYM=8.0D0*PHI2/(A+L/2)
C
C INITIALIZING MEAN ISOPCYNAL HEIGHT AND TOTAL HEIGHT
C
H0=1.0D0
HT=1.0D0
C
C THE SOLUTION FOR A-L/2 < Y < A+L/2
C
YABS=DABS(Y)
IF(YABS.GT.A)GOTO 701
100 FORMAT(' **** AIRY ROUTINE ERROR **** ',2I6,2F10.5)
200 FORMAT(' **** BIRY ROUTINE ERROR **** ',I6,2F10.5)
EXCP=(C-1.0D0)*C*L**2/(8.0D0*MU)
EXCPR=DREAL(EXCP)
EXCPI=DIMAG(EXCP)
MAG=(DSQRT(EXCPR**2+EXCPI**2))**(2.0D0/3.0D0)
ARG=2.0D0*(DATAN2(EXCPI,EXCPR))/3.0D0
EXCP=MAG*CDEXP(DCMPLX(0.0D0,ARG))
EXCP=EXCP*(K**2-1.0D0/C+8.0D0*MU*(Y-A)/(C*(C-1.0D0)*L**2))
ZTA=(2.0D0/3.0D0)*EXCP*CDSQRT(EXCP)
ZTAR=DABS(DREAL(ZTA))
EXCPR=DREAL(EXCP)
EXCPI=DIMAG(EXCP)
ARG=DABS(DATAN2(EXCPI,EXCPR))
IF(ARG.GE.PI)KODE=1
IF(ARG.LT.PI)KODE=2
CALL ZAIRY(EXCPR,EXCPI,0,KODE,AIR,AII,NZ,IER)
IF(IER.NE.0.OR.NZ.NE.0)WRITE(7,100)NZ,IER,C
CALL ZAIRY(EXCPR,EXCPI,1,KODE,DAIR,DAII,NZ,IER)
IF(IER.NE.0.OR.NZ.NE.0)WRITE(7,100)NZ,IER,C
IF(KODE.EQ.1)AIP=DCMPLX(AIR,AII)
IF(KODE.EQ.2)AIP=DCMPLX(AIR,AII)*CDEXP(-ZTA)
IF(KODE.EQ.1)DAIP=DCMPLX(DAIR,DAII)
IF(KODE.EQ.2)DAIP=DCMPLX(DAIR,DAII)*CDEXP(-ZTA)
CALL ZBIRY(EXCPR,EXCPI,0,2,BIR,BII,IER)
IF(IER.NE.0)WRITE(7,200)IER,C
CALL ZBIRY(EXCPR,EXCPI,1,2,DBIR,DBII,IER)
IF(IER.NE.0)WRITE(7,200)IER,C
BIP=DCMPLX(BIR,BII)*CDEXP(DCMPLX(ZTAR,0.0D0))
DBIP=DCMPLX(DBIR,DBII)*CDEXP(DCMPLX(ZTAR,0.0D0))
ETAY=AL1*AIP+AL2*BIP
HY=-8.0D0*(Y-A)*MU*ETAY/(L**2*(C-1.0D0))
H0=1.0D0-4.0D0*((Y-A)/L)**2
GO TO 703
701 CONTINUE
C
C THE SOLUTION FOR -B < Y < A-L/2
C
IF(Y.GT.A)GOTO 702
EP=CDEXP(CDSQRT(K**2-1.0D0/C)*(B+Y))
EM=CDEXP(-CDSQRT(K**2-1.0D0/C)*(B+Y))
SN=(EP-EM)/2.0D0
CN=(EP+EM)/2.0D0
ETAY=AL3*SN
HY=HYM+8.0D0*(Y-(A-L/2))/(L*PHASE)
GO TO 703
702 CONTINUE
C
C THE SOLUTION FOR A+L/2 < Y < 0
C
EAM=CDEXP(-CDSQRT(K**2-1.0D0/C)*Y)
EAP=CDEXP(CDSQRT(K**2-1.0D0/C)*Y)
ETAY=AL4*EAM+AL5*EAP
HY=HYM-8.0D0*(Y+A+L/2)/(L*PHASE)
703 CONTINUE
C

C
C THE SOLUTION FOR 0 < Y < D1
C
EAM2=CDEXP(-CDSQRT(K**2+S2/C)*Y)
EAP2=CDEXP(CDSQRT(K**2+S2/C)*Y)
ETAY=AL6*EAM2+AL7*EAP2
HY=HYM-8.0D0*(Y-D1)/(L*PHASE)
704 CONTINUE
C

C THE SOLUTION FOR D1 < Y < D2
EP3=CDEXP(CDSQRT(K**2+S3/C)*(Y-D2))
EM3=CDEXP(-CDSQRT(K**2+S3/C)*(Y-D2))
SN3=(EP3-EM3)/2.0D0
CN3=(EP3+EM3)/2.0D0
ETAY=AL8*SN3
HY=HYM-8.0D0*(Y-D2)/(L*PHASE)
705 CONTINUE

C COMPUTING THE REAL PART OF THE SOLUTIONS

C
ETA=DREAL(ETAY*PHASE)
H=DREAL(HY*PHASE)
HT=H0+H
IF(HT.LE.0.0D0)HT=0.0D0
IF(HT.LE.0.0D0)H=-H0
PRES0=Y+MU*H0
PRESP=MU*(H+ETA)
PREST=Y+MU*(HT+ETA)
RATIO=0.0D0
IF(ETA.NE.0.0D0)RATIO=H/(ETA*DSQRT(8.0D0*MU/L**2))
C

 

[alphysicist]

Hi logi,

I was asking if you could write up a short program that recreated the error you're seeing, by just isolating the part that causes the error.

There is one thing that caught my eye. The expression that you send to atan2 depend on different expressions, but many seem to depend upon a variable C. But I don't see where C is declared, and I don't see where it is given a value. Was the C in the equation supposed to be the CI that was declared? Even then, I don't see where CI is given a value. (But perhaps it's taken care of somehow in the part of the program that's not shown.)

If that doesn't fix the problem, I would suggest adding this line right before the two lines where the datan2 is used:

print*,excpi,excpr

When the program fails, you'll be able to tell what values had just been sent to datan2. If those values are both zero, then that will be the problem.

 

[logi]

Hi alphysicist,
Thank you so much for your reply. You are right. I didn't declared any value for cguess, that is why I got the problem. Now I am running my program. If I get any problem, I will reply to you.

Thanks again.
Regards
Logi

 

[logi]

Hi alphysicist ,
I have some problem again. I found the values for cguess.Now I have to find the exact C using this values. How can i do that? Pl help me out.
Thanks Logi

 

[alphysicist]

Hi logi,

Unfortunately I have no idea what cguess and C represent, and I don't know what your program is doing in its calculation. I'm not sure how I can help.


Is your program some kind of iteration process, where you plug back in the value of CR and CI to make a new cguess, and then run it again with the new cguess? If so, I think you would just keep repeating the process until you get to a cguess that doesn't change.