C********************************************************************* C* ***** RtStClb2.for ***** C Многослойные кулоновские структуры с дифференциальным вращением C Joseph J. Smulsky C Последнее изменение 18 ч. 00 м. 02.02.2015 г. C********************************************************************* C--------------------------------------------------- C Нужно задать значения параметров! PARAMETER (l2=1000,l4=6000) ! C --------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A1(l2,l2+1),fi(l2,l2),cm(l2), *R(l2),Ro(l2),qo(l2),Prd(l2),wo(l2), *omo(l2*l2+1),Xo(l2*l2+1,3),Uo(l2*l2+1,3),X1(l2*l2+1,3), *U1(l2*l2+1,3),Pns(3),Xc(3),AM1(3),Pn(3),X(l2*l2+1,3), *U(l2*l2+1,3),Sp(l2*l2+1,3),Et(l2*l2+1),Ra(l2*l2+1),dUp(l2*l2+1,3) INTEGER*4 N,Nstr,N2,N3,Idta,q CHARACTER*12 FIMA1,FIMA2,FIMA3,FIMA4,FIMA5,FIMA6,FIMA7 FIMA1 ='RtStClb2.dat' FIMA3 ='ErRtClSt.dat' FIMA4 ='RtQNcREr.dat' ! Номер кольца, радиус кольца, ошибка FIMA6 ='RtStCfi0.dat' FIMA7 ='RtStCfi1.dat' C--- Считывание файла RtStClb2.dat. OPEN(UNIT=1,FILE=FIMA1,STATUS='OLD') READ(1,39,ERR=40)N2 READ(1,39,ERR=40)N3 READ(1,39,ERR=40)nA READ(1,55,ERR=40)Tmp READ(1,54,ERR=40)cmp READ(1,54,ERR=40)cme READ(1,54,ERR=40)cmne READ(1,54,ERR=40)ee READ(1,54,ERR=40)ed READ(1,54,ERR=40)UnT READ(1,55,ERR=40)Tmp READ(1,54,ERR=40)okr READ(1,39,ERR=40)Kl1 READ(1,'(59X,A12)',ERR=40)FIMA2 READ(1,38,ERR=40)Ivr READ(1,38,ERR=40)Kl3 READ(1,54,ERR=40)Prdi READ(1,54,ERR=40)eps0 READ(1,54,ERR=40)dt READ(1,'(59X,A12)',ERR=40)FIMA5 READ(1,39,ERR=40)Idta CLOSE(1) 38 FORMAT(59X,I5) 39 FORMAT(59X,I12) 53 FORMAT(59X,F11.4) 54 FORMAT(59X,D11.5) 55 FORMAT(2X,D10.4) GO TO 1 40 WRITE(*,*)'*01Rt Error reading file RtStClb2.dat *** ', *N2,N3,nA,cmp,cme,cmne,ee,ed,UnT,okr,Kl1,Ivr,Kl3,Prdi, *eps0,FIMA5,Idta OPEN(UNIT=3,FILE=FIMA3) WRITE(3,*)'*01Rt Error reading file RtStClb2.dat *** ', *N2,N3,nA,cmp,cme,cmne,ee,ed,UnT,okr,Kl1,Ivr,Kl3,Prdi, *eps0,FIMA5,Idta CLOSE(3) GO TO 150 1 pi = 3.1415926535897932385D+00 p2 = 2.0D+00*pi DZI = 0.0D+00 PND = 1.5D+00 PND2 = 1.0D+00/3.0D+00 dfi0 = p2/N3 C--- Constants from Handbook by Yavorsky & Detlaf: pp.749, 910, 912, 913. C cmp = 1.67252D-27 C cme = 9.10910D-31 C cmne = 1.67482D-27 C ee = 4.80298D-10 C Prd0 = 3.294713055979042D-17 C Prd0 = 1.0D-15 Prdsc = Prdi*UnT nZ = N2*N3 qo0 = nZ qo1 = -1 N23 = N2*N3 Nne = nA - nZ cm0 = nZ*cmp + Nne*cmne Rn0 = 1.5D-15 Re = 1.5D-15 C An = 2*Zn Rn = Rn0*(nA**PND2) TMP = 1.0D-9*ee*ee*Prdsc*Prdsc/(4*pi*pi*ed*cme) R0 = TMP**PND2 C okt = 1/(100*Prdsc) okt = 1/UnT C WRITE(*,*)'*01*** ',R0,Prdi,Prdsc ssm = cm0 DO J = 1,N2 cm(J) = cme ssm = ssm + cm(J)*N3 !!-1str C WRITE(*,*)'*01 cm0,ssm*** ',cmp,cmne,Nne,nZ,cm0,ssm END DO ! J TMP = 1.0D-9*ee*ee/(ed*ssm*okt*okt) Am = TMP**PND2 okv = 1/(okt*Am) N = N2*N3 + 1 Nstr = l4/2 k1 = 1 fn3 = DZI DO l = 2,N3 fn3 = fn3 + 0.25D+0/DSIN((l-1)*pi/N3) END DO ! l C--- Углы первых тел на кольцах IF(Ivr.EQ.1) THEN DO I = 1,N2 fi(I,1) = DZI END DO ! I ELSE I2 = 0 DO I = 1,N2 I2 = I2 +1 fi(I2,1) = DZI I2 = I2 +1 fi(I2,1) = 0.5D+0*dfi0 IF(I2.EQ.N2) GO TO 3 END DO ! I 3 END IF ! Ivr C--- Углы первых тел на кольцах из файла RtStCfi0.dat. IF(Kl3.EQ.1) THEN OPEN(UNIT=6,FILE=FIMA6) DO I = 1,N2 READ(6,*)I1,fi(I,1) END DO ! I CLOSE(6) ELSE OPEN(UNIT=7,FILE=FIMA7) DO I = 1,N2 WRITE(7,*)I,fi(I,1) END DO ! I CLOSE(7) END IF ! C--- Углы всех тел на кольцах DO J = 1,N2 DO l = 2,N3 fi(J,l) = fi(J,l-1) + dfi0 END DO ! l END DO ! J C WRITE(*,*)'*01*** ' C--- Радиусы колец, отнесенные к R0. IF(Kl1.EQ.1) GO TO 27 R(1) = (nZ - fn3)**PND2 OPEN(UNIT=3,FILE=FIMA3) WRITE(3,*) 'Devided on R0 radii of the circles' WRITE(3,*) ' 1 ',R(1) DO I = 2,N2 R(I) = okr*I*R(1) WRITE(3,*) I,R(I) END DO ! I CLOSE(3) GO TO 11 27 OPEN(UNIT=2,FILE=FIMA2) DO J = 1,N2 READ(2,*)J1,R(J) END DO ! J CLOSE(2) C--- Cила Qs воздействия на тело кольца j всех тел всех колец C--- В RtStClb2_0.for Qs - на кольцо j не учитывалось действие нар. колец 11 DO J = 1,N2 As = DZI C DO I = 1,J ! Было в RtStClb2_0.for DO I = 1,N2 Qs = DZI RIJ = R(I)/R(J) DO l = 1,N3 dfi = fi(I,l) - fi(J,1) C--- Предотвращение деления на нуль в случаях C--- сил воздействия тела на себя. IF(l.EQ.1.AND.I.EQ.J) dfi = pi C--- Расчет сил воздействия на тело данного кольца тела другого кольца TMP = RIJ*DCOS(dfi) Q1 = (1-TMP)/(1+RIJ*RIJ-2*TMP)**PND C--- Обнуление сил воздействия тела кольца на себя. IF(l.EQ.1.AND.I.EQ.J) Q1 = DZI Qs = Qs + Q1 C WRITE(*,*)'J = ',J,' I = ',I,Q1 END DO ! l A1(J,I) = Qs As = As + A1(J,I) END DO ! I R3 = R(J)*R(J)*R(J) TMP = nZ - As IF(TMP.GT.DZI) GO TO 12 OPEN(UNIT=3,FILE=FIMA3) WRITE(*,*)'*02Rt Error: (nZ - As) < 0 *** ',TMP,' J = ',J, ! 02.02.2015 *' R(J) = ',R(J) OPEN(UNIT=3,FILE=FIMA3) WRITE(3,*)'*02Rt Error: (nZ - As) < 0 *** ',TMP,' J = ',J, ! 02.02.2015 *' R(J) = ',R(J) CLOSE(3) GO TO 150 12 Prd(J) = DSQRT(R3/TMP) wo(J) = p2/(Prd(J)*Prdsc*okt) END DO ! J OPEN(UNIT=3,FILE=FIMA3) WRITE(3,*)'Relative radii (R(I)/R0) of the circles and periods' DO I = 1,N2 WRITE(3,*) I,R(I),Prd(I) END DO ! I WRITE(3,*) 'Coefficients A1(J,I) of the equations' DO J = 1,N2 WRITE(3,15) J,(A1(J,I),I=1,N2) END DO ! J CLOSE(3) 15 FORMAT(I5,20E13.4) C--- Массы тел в программе Galactica и параметры в плоскости орбиты omo(1) = cm0/ssm qo(1) = qo0 Ra(1) = Rn/Am k = 1 DO J = 1,N2 Ro(J) = R(J)*R0/Am DO l = 1,N3 k = k+1 omo(k) = cm(J)/ssm qo(k) = qo1 Ra(k) = Re/Am Xo(k,1) = Ro(J)*DCOS(fi(J,l)) Xo(k,2) = Ro(J)*DSIN(fi(J,l)) Xo(k,3) = DZI Uo(k,1) = -Ro(J)*wo(J)*DSIN(fi(J,l)) Uo(k,2) = Ro(J)*wo(J)*DCOS(fi(J,l)) Uo(k,3) = DZI END DO ! l END DO ! J C--- Координаты и скорости периферийных тел в экваториальной плокости DO k = 2,N X1(k,1) = Xo(k,1) X1(k,2) = Xo(k,2)*DCOS(eps0) - Xo(k,3)*DSIN(eps0) X1(k,3) = Xo(k,2)*DSIN(eps0) + Xo(k,3)*DCOS(eps0) U1(k,1) = Uo(k,1) U1(k,2) = Uo(k,2)*DCOS(eps0) - Uo(k,3)*DSIN(eps0) U1(k,3) = Uo(k,2)*DSIN(eps0) + Uo(k,3)*DCOS(eps0) END DO ! k C--- Координаты и скорости центрального тела в экваториальной плокости DO q = 1,3 X1(1,q) = DZI U1(1,q) = DZI END DO ! q C--- Параметры центра масс. 30 DO q = 1,3 Pns(q) = DZI Xc(q) = DZI AM1(q) = DZI END DO ! q DO k = k1,N IF(omo(k).EQ.DZI) GO TO 32 DO q = 1,3 Xc(q) = Xc(q) + omo(k) * X1(k,q) Pn(q) = omo(k) * U1(k,q) Pns(q) = Pns(q) + Pn(q) END DO ! q 32 END DO ! k C--- Приведение к центру масс. DO k = k1,N IF(omo(k).EQ.DZI) GO TO 16 DO q = 1,3 X(k,q) = X1(k,q) - Xc(q) U(k,q) = U1(k,q) - Pns(q) END DO ! q AM1(1) = AM1(1)-omo(k)*U(k,2)*X(k,3)+omo(k)*U(k,3)*X(k,2) AM1(2) = AM1(2)+omo(k)*U(k,1)*X(k,3)-omo(k)*U(k,3)*X(k,1) AM1(3) = AM1(3)+omo(k)*U(k,2)*X(k,1)-omo(k)*U(k,1)*X(k,2) 16 END DO ! k C WRITE(*,*)' 1. AM1(1) = ',Am1(1),Amn,vkn DO q = 1,3 Pns(q) = DZI END DO ! q omm = DZI DO k = k1,N DO q = 1,3 Pn(q) = omo(k) * U(k,q) Pns(q) = Pns(q) + Pn(q) dUp(k,q) = DZI Sp(k,q) = DZI END DO ! q Et(k) = DZI IF(omo(k).LT.omm) GO TO 13 omm = omo(k) 13 END DO ! k C WRITE(*,*)' 2. Ra(1) = ',Ra(1),Amn,vkn A = 1.0D+0 B = 0.5D+0 C = 0.25D+0 TMP = (N+1)**PND2 Mu = NINT(TMP) C WRITE(*,*)TMP,' MU = ',MU T = DZI Um = DZI Spsx = DZI Spsy = DZI Spsz = DZI E = DZI Em = DZI Ett = DZI dtp = DZI C--- Создание файла начальных условий FIMA5 для программы Galactica. OPEN(UNIT=5,FILE=FIMA5) WRITE(5,*)T,omm,Um,dtp,Pns(1),Pns(2),Pns(3),AM1(1),AM1(2), *AM1(3),Spsx,Spsy,Spsz,E,Em,Ett,dt,1,1,N,A,B,C,Mu,ed DO k = k1,N WRITE(5,*) omo(k),(X(k,q),q=1,3),(U(k,q),q=1,3), *(dUp(k,q),q=1,3),(Sp(k,q),q=1,3),Ra(k),Et(k),qo(k) END DO ! k C--- Information line at end of file WRITE(5,*)Idta,N,ssm,Am,okv,0,0,0,0,0,0,0,0,0,0,0 CLOSE(5) 150 CONTINUE END