C $STORAGE:4 C $NODEBUG C $NOFLOATCALLS COMMON/GM1/GMBAR(2) DIMENSION VOL(2),AREA(2), RVOL(2),GBAR(2) 1 , VOL0(2),AREA0(2),RVOL0(2),GBAR0(2),GMBAR0(2) 2 , RATV(2),RATA(2),RATR(2),RATG(2),RATGM(2), 3 omega(2) DATA PI /3.141593/ 1 WRITE(*,2) 2 FORMAT(' ENTER Q,F: (F=0. KILLS)',$) READ(*,*)Q,F IF (F.LE.0.) STOP WRITE(*,3) 3 FORMAT(' ENTER NX,NT:',$) READ(*,*)NX,NT X1 = 0.6 CALL XL1(Q,C1,X1,IER1) X2 = 2.0 - X1 CALL XL1(Q,C2,X2,IER2) C = C1 + (F-1.0)*(C2-C1) IF (F.LT.1.0.AND.F.GT.0.) C = C1/F T=0. X=0. Z11 = RHO (C,Q,X,T) T=PI*0.5 Y11 = RHO (C,Q,X,T) X=1.0 Y12 = RHO (C,Q,X,T) T = 0. Z12 = RHO (C,Q,X,T) WRITE(*,311)Y11,Z11,Y12,Z12 311 FORMAT(' Y1,Z1,Y2,Z2 = ',4F11.5) CALL CENVOL(Q,F,C,VOL,AREA,GBAR,X1,X2,NX,NT) QP1 = Q+1.0 QPD = 1.0/QP1 QQPD = Q*QPD RVOL0(1)=2.0*QPD/(C-2.0*QQPD-QQPD*QQPD+4.0*QPD*QPD/C**2) RVOL0(2)=2.0*QQPD/(C-2.0*QPD-QPD*QPD+4.0*QQPD*QQPD/C**2) GBAR0(1) = 2.0/(QP1*RVOL0(1)**2) GBAR0(2) = 2.0*Q/(QP1*RVOL0(2)**2) DO 10 IK=1,2 omega(ik) = (C-QQPD)*0.5*qp1 RVOL(IK) = (0.75*VOL(IK)/3.141593)**0.3333333 GMBAR0(IK) = 1.0 / GBAR0(IK) AREA0(IK) = 4.0*3.141593*(RVOL0(IK))**2 RATR(IK) = RVOL(IK)/RVOL0(IK) RATA(IK) = AREA(IK) / AREA0(IK) RATG(IK) = GBAR(IK) / GBAR0(IK) RATGM(IK) = GMBAR(IK) /GMBAR0(IK) 10 CONTINUE IUNIT = 6 4 WRITE(*,6)Q,F,C,IER1,IER2,C1,X1,C2,X2,VOL,AREA, 1 RVOL,NX,NT, 1 GBAR,GMBAR,omega,(RATR(IK),RATA(IK),RATG(IK),RATGM(IK),IK=1,2) 6 FORMAT(' Q=',F6.3,12X,'F=',F6.3,12X,'C=',F9.5,9X, 1 'ERR.CODES=',2I2 / ' C1=',F9.5,9X,'X1=',F9.5, 2 9X,'C2=',F9.5,9X,'X2=',F9.5 / 3 ' VOLUMES=',2F10.6,T40,'AREAS=',2F10.6 / # ' VOL.RADII= ',2F13.7 ,T40,'NX,NT=',2I5 / 4 ' AVE.GRAV.=',2F13.5,T40,'AVE.INV.GRAV=',2F13.6/ Z ' OMEGA VALUES = ',2F13.6/ A ' PR: ', 5 ' R/R0 =',F9.5,' S/S0=',F9.5,' G/G0=',F9.5,' GM/GM0=',F9.5/ B ' SEC:', 5 ' R/R0 =',F9.5,' S/S0=',F9.5,' G/G0=',F9.5,' GM/GM0=',F9.5) WRITE(*,7) 7 FORMAT('---------------------------------------------', 1 '------------------------------------'/) GO TO 1 C IF(IUNIT.EQ.3)GO TO 1 C IUNIT = 3 C GO TO 4 END SUBROUTINE CENVOL(Q,F,C,VOL,AREA,GBAR,X1,X2,NXIN,NT) C C Computation of Roche common (and detached) envelope volumes by C double Simpson's Rule quadrature. Surface areas and mean gravities C also computed. C C This version: Inst. of Astron. 12 July 1983. C Modified for MS-FORTRAN 77 v. 3.2 June 1984 C C Written by: Stefan W. Mochnacki C REAL*8 VSUM,DV,CARE,SARE,DDT,DGB,ARSUM,GSUM,DA,DG, 1 RSQ,DS,TE,GMSUM,DGBM,DGM,DXF DIMENSION VOL(2),XLE(2),DXS(2),RMA(2),XSS(2),NXS(2), 1 R(61),DS(61),TE(61),CARE(400),SARE(400),AREA(2),GBAR(2), 2 DGB(400), COST(61),SINT(61),GMBAR(2),DGBM(400),RSQ(61) 3 , DXF(400) COMMON/GM1/GMBAR DATA PI/3.141593/, PI2/1.570796/ C Compute End points of components along x-axis in cylindrical scheme. IF (F.GT.1.99999)GO TO 80 XA = XEND(C,Q,-0.25) XB = XEND(C,Q,1.10) GO TO 81 80 XB = X2 IF(Q.EQ.1.0)XA=1.0-XB IF(Q.NE.1.0)XA=XEND(C,Q,-0.25) XAI=X1 XBI=X1 GO TO 82 81 IF(F.LT.0.99999)GO TO 83 XAI=X1 XBI=X1 GO TO 82 83 XAST = MIN(0.49,-XA) XAI = XEND(C,Q,XAST) XAST = MAX (0.51,2.0-XB) XBI = XEND(C,Q,XAST) C Set up for integration. 82 XIN = 0. WRITE(*,900)XA,XAI,XBI,XB 900 FORMAT(' XA,XAI,XBI,XB=',4F12.5) FU = 2.0 / (1.0+Q) C WRITE(3,900)XA,XAI,XBI,XB XSS(1) = XA XSS(2) = XBI FUQ = Q*FU XLE(1) = XAI - XA XLE(2) = XB - XBI NTP = NT/4 C FORCE NTP ODD> NTP = ntp/2 ntp = ntp*2+1 NTT = NTP - 1 NTT = NTP-1 NX = NXIN DX = (XLE(1) +XLE(2)) / NX NTM = NTP-2 C FORCE NXS(I) ODD = NO. OF POINTS,INCL. ENDS IN X-COORD. DO 4 I=1,2 NXS(I) = XLE(I)/DX ITE = NXS(I) / 2 ITE = ITE*2 + 1 ITE = MAX0(ITE,31) C NXS(I) HAS MIN VAL. = 31. ITM = ITE - 1 DXS(I) = XLE(I) /ITM 4 NXS(I) = ITE + 14 NX = NXS(1) + NXS(2) DT = PI2 / NTT DDT = DT DO 40 J=1,NTP RJ = J - 1 T = DT*RJ SINT(J) = SIN(T) 40 COST(J) = COS(T) C Loop for each component IK = 1 , 2 DO 8 IK=1,2 IBEG = (IK-1) * NXS(1) + 1 NIX = NXS(IK) + IBEG - 1 C Calculate cross-section slices CARE, SARE, DGB. DO 12 I=IBEG,NIX CARE(I) = 0. SARE(I) = 0. DGB(I) = 0. DGBM(I)=0. IF(IK.EQ.1.AND.I.EQ.IBEG) GO TO 12 IF(IK.EQ.1.AND.I.EQ.NIX.AND.F.LE.1.0)GO TO 12 IF(IK.EQ.2.AND.I.EQ.IBEG.AND.F.LE.1.0)GO TO 12 IF(IK.EQ.2.AND.I.EQ.NIX)GO TO 12 IBG = I-IBEG INX = NIX-I IF (IBG.LE.7.OR.INX.LE.7) GO TO 41 RI = IBG - 7 X = RI*DXS(IK) + XSS(IK) DXF(I) = DXS(IK) GO TO 42 41 DXF(I) = DXS(IK)*0.125 IF(IBG.LE.7) X = XSS(IK) + IBG*DXF(I) IF(INX.LE.7) X = XSS(IK) + XLE(IK) - INX*DXF(I) 42 CONTINUE ARSUM = 0. GSUM = 0. VSUM = 0. GMSUM = 0. XSQ = X*X XM = X - 1.0 XMSQ = XM*XM C Angular-dependant values. DO 9 J=1,NTP RJ = J-1 T = DT*RJ RR = RHO(C,Q,X,T) R(J)=RR IF (RR.LE.1.0E-30.OR.RR.GT.1.0) GO TO 90 GO TO 92 90 R(J) = 0.0 TE(J) = 0.0 DS(J) = 0.0 WRITE(*,91) I,J,T,Q,X,F,RR 91 FORMAT(' ERROR IN RHO: I,J,T,Q,X,F,RHO=',2I5,5F10.5) 92 RRS = RR*RR RSQ(J) = DBLE(RRS) R1S = X*X + RRS XM = X-1.0 R2S = XM*XM + RRS R1 = SQRT(R1S) R2 = SQRT(R2S) E1 = -FU/(R1S*R1) E2 = -FUQ/(R2S*R2) FXS = E1*X + E2*XM + 2.0*X - FUQ FYS = RR*SINT(J) * ( E1 + E2 + 2.0 ) FZS = RR*COST(J) * ( E1 + E2 ) Z = SQRT( FXS*FXS + FYS*FYS + FZS*FZS ) TE(J) = Z DS(J) = 0. DIVS = FYS*SINT(J) +FZS*COST(J) IF(DIVS.EQ.0.)WRITE(6,93)J,T,Q,X,F,RR 93 FORMAT(' DS ERROR: J,T,Q,X,F,RR=',I5,5F10.5) IF(DIVS.EQ.0.)GO TO 9 DS(J) = -Z*RR / DIVS 9 CONTINUE C Integrate w.r.t. Theta by Simpson's Rule. DO 10 J=1,NTM,2 JP = J + 1 JPP = J + 2 DV = RSQ(J) + 4.0*RSQ(JP) + RSQ(JPP) DA = DS(J) + 4.0*DS(JP) + DS(JPP) DG = TE(J)*DS(J) + 4.0*TE(JP)*DS(JP) + TE(JPP)*DS(JPP) DGM = DS(J)/TE(J) + 4.0*DS(JP)/TE(JP) + DS(JPP)/TE(JPP) GMSUM = GMSUM + DGM ARSUM = ARSUM + DA GSUM = GSUM + DG 10 VSUM = VSUM + DV SARE(I) = 1.333333333333*ARSUM*DDT DGB(I) = 1.333333333333*GSUM*DDT DGBM(I) = 1.333333333333*GMSUM*DDT CARE(I) = 0.6666666666667*VSUM*DDT 12 CONTINUE C EXTRAPOLATE S,G AT ENDS. NIXM = NIX-2 NIXM1 = NIX - 1 IB1 = IBEG + 1 IB2 = IBEG + 2 IF(IK.EQ.2.AND.F.GT.1.0) GO TO 120 SARE(IBEG) = 2.0*SARE(IB1) - SARE(IB2) DGB(IBEG) = 2.0*DGB(IB1) - DGB(IB2) DGBM(IBEG) = 2.0*DGBM(IB1) - DGBM(IB2) 120 IF(IK.EQ.1.AND.F.GT.1.0) GO TO 121 SARE(NIX) = 2.0*SARE(NIXM1) - SARE(NIXM) DGB(NIX) = 2.0*DGB(NIXM1) - DGB(NIXM) DGBM(NIX) = 2.0*DGBM(NIXM1) - DGBM(NIXM) 121 CONTINUE C INTEGRATE W.R.T. X BY SIMPSON'S RULE. ARSUM = 0. GSUM = 0. GMSUM = 0. VSUM = 0. GBAR(IK) = 0. GMBAR(IK) = 0. DO 15 I=IBEG,NIXM,2 IP = I + 1 IPP = I + 2 DA = SARE(I) + 4.0*SARE(IP) + SARE(IPP) ARSUM = ARSUM + DA*DXF(IP) DG = DGB(I) + 4.0*DGB(IP) + DGB(IPP) DGM = DGBM(I) + 4.0*DGBM(IP) + DGBM(IPP) GSUM = GSUM + DG*DXF(IP) GMSUM = GMSUM + DGM*DXF(IP) 15 VSUM = VSUM + (CARE(I)+CARE(IP)*4.0+CARE(IPP))*DXF(IP) IF(DABS(ARSUM).LE.1.0D-30)GO TO 150 GBAR(IK) = GSUM / ARSUM GMBAR(IK) = GMSUM / ARSUM 150 AREA(IK) = 0.333333333333*ARSUM 8 VOL(IK) = 0.3333333333333*VSUM RETURN END SUBROUTINE XL1(Q,C1,X1,IER) C CALCULATION OF INNER & OUTER LAGRANGIAN POINTS, AND CRITICAL POTENTIALS. IMAX=15 EPS=2.E-6 A0=-2./(1.+Q) U=Q/(1.+Q) UU=U*U IT=0 IER=0 1 IT=IT+1 XX=X1*X1 XM=X1-U XN=ABS(X1-1.) E1=A0/(XX*X1) E2=-2.*U/(XN*XN*XN) FX=X1*(2.+E1+E2)-E2-2.*U D=2.*(1.-E1-E2) IF(ABS(D)-EPS)3,3,2 2 DX=-FX/D X1=X1+DX CG=ABS(DX) IF(IMAX-IT)3,5,5 5 IF(CG-EPS)4,4,1 3 IER=1 RETURN 4 XM=X1-U XN=ABS(X1-1.) C1=-A0/X1+2.*U/XN+XM*XM RETURN END FUNCTION XEND(CIN,QIN,XIN4) C EVALUATION OF NON-SINGULAR END-POINTS. IMPLICIT REAL*8 (A-H,O-Z) REAL*4 XEND,CIN,QIN,XIN4 C=DBLE(CIN) Q=DBLE(QIN) XIN=DBLE(XIN4) A0=-2./(1.+Q) U=Q/(1.+Q) IMAX=25 UU=U*U IT=0 EPT=1.E-30 EPS = 5.0E-6 / C C *** SPECIAL SECTION TO GET APPROXIMATE INITIAL VALUE... BRUTE FORCE. DX = -0.03 / C IF(XIN.GT.1.0.OR.(XIN.GT.0.0.AND.XIN.LT.0.5)) DX=-DX XTEST = 0. IF (XIN.GT.0.5) XTEST = 1.0 10 XTEST = XTEST + DX XM = XTEST - U XA = 1.0/ABS(XTEST) XB = 1.0/ABS(XTEST-1.0) F = -A0*XA + 2.0*U*XB +XM*XM -C IF (F.GT.0.)GO TO 10 C WRITE(*,936)XTEST,DX,F C936 FORMAT(' XTEST,DX,F(DIFF)=',3F12.5) C WRITE(3,936)XTEST,DX,F XEND8 = XTEST - 0.5*DX C NEWTON-RAPHSON ITERATION FOLLOWS. 1 IT=IT+1 XM=XEND8-U X1=XEND8-1. XA=1./ABS(XEND8) XB=1./ABS(X1) X2=XA*XA XT=XB*XB F=-A0*XA+2.*U*XB+XM*XM-C FX=A0*X2*XA*XEND8-2.*U*X1*XB*XT+2.*XM IF(ABS(FX)-EPT)3,3,2 2 DX=-F/FX XEND8=XEND8+DX IF(IMAX-IT)3,3,5 5 IF(ABS(DX)-EPS)4,4,1 3 XEND8=0. 4 XEND = SNGL(XEND8) RETURN END FUNCTION RHO(CIN,QIN,XIN,THIN) IMPLICIT REAL*8 (A-H,O-Z) REAL*4 CIN,QIN,XIN,THIN,RHO C ROCHE CYLINDRICAL RADIUS CALCULATOR. C=DBLE(CIN) Q=DBLE(QIN) X=DBLE(XIN) TH=DBLE(THIN) C IMPROVED VERSION: AUGUST 1983, IOA. IMAX=40 IC=0 RHODST=0.5/C STH=SIN(TH) STH=STH*STH RK=1./(1.+Q) RL=Q*RK RM=X-RL RM=RM*RM C DAMPENED NEWTON-RAPHSON ITERATION. RHOD = RHODST EPS = RHODST*5.0E-5 DELMAX = 0.2 C WRITE(*,400)RHODST,EPS C 400 FORMAT(' RHO(START)=',F10.6,' EPS=',1PE14.6) 1 RR=RHOD*RHOD RS=X*X+RR RT=RS-2.*X+1. R1=SQRT(RS) R2=SQRT(RT) IC=IC+1 F = 2.*RK/R1 + 2.*RL/R2 + RM + RR*STH - C RS = R1*RS RT = R2*RT FRHOD = -2.0*RHOD*(RK/RS+RL/RT-STH) IF(ABS(FRHOD)-1.0E-30)3,3,4 4 DEL = -F/FRHOD SDEL = DEL ADEL = MIN(0.8*RHOD,ABS(DEL)) DEL = SIGN(ADEL,DEL) DEL = MIN(DELMAX,DEL) RHOD = RHOD + DEL C WRITE(*,100)IC,RHOD,DEL,SDEL C 100 FORMAT(' ITER=',I3,' RHO=',F9.5,' DEL=',1PE12.5,' SDEL=',1PE12.5) IF(IMAX-IC)3,3,5 5 IF(ADEL-EPS)2,2,1 3 RHOD=0. 2 RHO = SNGL(RHOD) RETURN END