C********************************************************************* C* ***** RtCrcSt2.for ****** C Вращающиеся кольцевые структуры C Joseph J. Smulsky C Последнее изменение 18 ч. 00 м. 04.06.2014 г. !!-1str C C********************************************************************* C--------------------------------------------------- C Нужно задать значения параметров! PARAMETER (l2=1000,l4=6000) ! C --------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A1(l2,l2+1),A2(l2,l2+1),fi(l2,l2),cm(l2), *R(l2),Ro(l2),cmo(l2),dmo(l2),dmm(l4), *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,Lt,Ltk,Lt4,Lt5,Idta,q CHARACTER*12 FIMA1,FIMA2,FIMA3,FIMA4,FIMA5,FIMA6,FIMA7 CHARACTER*3 Lt5Ima CHARACTER*4 ITER FIMA1 ='RtCrcSt2.dat' FIMA3 ='ErRtCrSt.dat' FIMA4 ='RtNcRMEr.dat' ! Номер кольца, радиус, масса тела, ошибка FIMA6 ='RtStrfi0.dat' FIMA7 ='RtStrfi1.dat' C--- Считывание файла RtCrcSt2.dat. OPEN(UNIT=1,FILE=FIMA1,STATUS='OLD') READ(1,39,ERR=40)N2 READ(1,39,ERR=40)N3 READ(1,54,ERR=40)okr READ(1,38,ERR=40)Kl1 READ(1,'(59X,A12)',ERR=40)FIMA2 READ(1,54,ERR=40)smi READ(1,54,ERR=40)pm0 READ(1,38,ERR=40)Ivr READ(1,38,ERR=40)Kl3 READ(1,54,ERR=40)Prd READ(1,54,ERR=40)EPS READ(1,39,ERR=40)Ltk READ(1,54,ERR=40)eps0 READ(1,54,ERR=40)Rob 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) GO TO 1 40 WRITE(*,*)'*01Rt Error reading file RtCrcSt2.dat *** ', *N2,N3,okr,Kl1,FIMA2,smi,pm0,Ivr,Kl3,Prd,EPS,Ltk,eps0, *Rob,dt,FIMA5,Idta OPEN(UNIT=3,FILE=FIMA3) WRITE(3,*)'*01Rt Error reading file RtCrcSt2.dat *** ', *N2,N3,okr,Kl1,FIMA2,smi,pm0,Ivr,Kl3,Prd,EPS,Ltk,eps0, *Rob,dt,FIMA5,Idta CLOSE(3) GO TO 150 1 pi = 3.1415926535897932385D+00 p2 = 2.0D+00*pi G = 6.67259D-11 Yrs = 365.25636042D+00 dsc = 24D+00*3600D+00 DZI = 0.0D+00 PND = 1.5D+00 PND2 = 1.0D+00/3.0D+00 dfi0 = p2/N3 Prdsc = Prd*Yrs*dsc w = p2/Prdsc w2 = w*w fn3 = DZI Vs = DZI N = N2*N3 + 1 Nstr = l4/2 Lt4 = 1 Lt5 = 1 k1 = 1 C okm = 0.05D+00 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--- Углы первых тел на кольцах из файла RtStrfi0.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--- Масса центрального тела и предварительная масса тела первого кольца cm0 = smi*pm0 cmo0 = pm0 sm = smi - cm0 cm(1) = sm/N3 C--- Радиусы колец IF(Kl1.EQ.1) GO TO 27 R(1) = (G*(cm0 + cm(1)*fn3)/w2)**PND2 OPEN(UNIT=3,FILE=FIMA3) WRITE(3,*) 'Initial masses of the central and the body of the', *' 1-st circle ',cm0,cm(1) WRITE(3,*) '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) OPEN(UNIT=3,FILE=FIMA3) WRITE(3,*) 'Initial masses of the central and the body of the', *' 1-st circle ',cm0,cm(1) WRITE(3,*) 'Radii of the circles' WRITE(3,*) ' 1 ',R(1) DO I = 2,N2 WRITE(3,*) I,R(I) END DO ! I CLOSE(3) 11 okt = 1/(100*Prdsc) wo = w/okt Am = R(1) C--- Относительные радиусы тел DO I = 1,N2 Ro(I) = R(I)/Am END DO ! I TMP = Am*Am*w2 TMP = TMP*Am Ci = TMP/(G*smi) C--- Углы всех тел на кольцах DO J = 1,N2 DO l = 2,N3 fi(J,l) = fi(J,l-1) + dfi0 END DO ! l END DO ! J C--- Cила Qs воздействия на тело кольца j всех тел всех колец DO J = 1,N2 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 END DO ! l A1(J,I) = Qs A2(J,I) = Qs END DO ! I TMP = Ro(J)*Ro(J) A1(J,N2+1) = Ci*TMP*Ro(J) - cmo0 A2(J,N2+1) = A1(J,N2+1) END DO ! J CALL GAUSS(N2,A1,cmo,S,l2) C--- Уточнение масс DO Lt = 1,Ltk DO J = 1,N2 Qs = DZI DO I = 1,N2 Qs = Qs + cmo(I)*A2(J,I) A1(J,I) = A2(J,I) END DO ! I A1(J,N2+1) = A2(J,N2+1) - Qs Vs = Vs + DABS(A1(J,N2+1)) C WRITE(*,*) '3',' Lt = ',Lt,J,cmo(J),A1(J,N2+1) END DO ! J Vs = Vs/N2 IF(Vs.LT.EPS) GO TO 5 CALL GAUSS(N2,A1,dmo,S,l2) dmm0 = DZI DO J = 1,N2 dmoo = dmo(J)/cmo(J) dmm0 = dmm0 + DABS(dmoo) END DO ! J dmm(Lt4) = dmm0/N2 IF(dmm(Lt4).LT.EPS) GO TO 5 DO J = 1,N2 C Ослабление приращения dmo(J) cmo(J) = cmo(J) + 0.1D+00*dmo(J) END DO ! J C--- Запись разностей итерационного процесса Lt4 = Lt4 + 1 IF(Lt.EQ.Ltk) GO TO 7 IF(Lt4.LE.Nstr) GO TO 9 7 WRITE(Lt5Ima,'(I3)') Lt5 ITER = Lt5Ima//'i' OPEN (UNIT=2,FILE=ITER) DO in = 1,Lt4-1 inLt = (Lt5-1)*Nstr + in WRITE(2,*)inLt,dmm(in) END DO ! in CLOSE(2) Lt5 = Lt5 + 1 Lt4 = 1 9 END DO ! Lt GO TO 15 5 WRITE(Lt5Ima,'(I3)') Lt5 ITER = Lt5Ima//'i' OPEN (UNIT=2,FILE=ITER) DO in = 1,Lt4 inLt = (Lt5-1)*Nstr + in WRITE(2,*)inLt,dmm(in) END DO ! in CLOSE(2) 15 ssm = cm0 DO J = 1,N2 cm(J) = cmo(J)*smi ssm = ssm + cm(J)*N3 !!-1str END DO ! J C--- Запись номера кольца, радиуса, массы тела C--- и ошибка в файл RtNcRMEr.dat. 17 OPEN(UNIT=4,FILE=FIMA4) DO I = 1,N2 WRITE(4,*)I,R(I),cm(I),dmo(I) END DO ! I C--- Массы тел в программе Galactica и параметры в плоскости орбиты 23 omo(1) = cm0/ssm TMP = G*ssm/(okt*okt) Am = TMP**PND2 okv = 1/(okt*Am) Robo = Rob*Am*Am*Am/ssm k = 1 DO J = 1,N2 Ro(J) = Ro(J)*R(1)/Am DO l = 1,N3 k = k+1 omo(k) = cm(J)/ssm 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*DSIN(fi(J,l)) Uo(k,2) = Ro(J)*wo*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 TMP = 3*omo(k)/(4*pi*Robo) IF(TMP.GE.DZI) GO TO 19 OPEN(UNIT=3,FILE=FIMA3) WRITE(3,*)'*05Rt Computing Error of Programm RtCrcSt2 *** ', *' after terminal step the mass of body k = ',k,' omo(k) = ', *omo(k),' is negative.' WRITE(*,*)'*05Rt Computing Error of Programm RtCrcSt2 *** ', *' after terminal step the mass of body k = ',k,' omo(k) = ', *omo(k),' is negative.' CLOSE(3) GO TO 150 19 Ra(k) = TMP**PND2 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 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) 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 CLOSE(5) 150 CONTINUE END C --- Подпрограмма решения систем линейных уравнений методом Гаусса GAUSS SUBROUTINE GAUSS(N,A,X,S,l2) DIMENSION A(l2,l2+1),X(l2) DOUBLE PRECISION A,X,S,S1,R N1 = N+1 DO 25 K = 1,N K1 = K+1 S = A(K,K) J = K DO 21 I = K1,N R = A(I,K) IF(DABS(R).LE.DABS(S)) GOTO 21 S = R J = I 21 CONTINUE IF(S.EQ.0) RETURN IF(J.EQ.K) GOTO 23 DO 22 I = K,N1 R = A(K,I) A(K,I) = A(J,I) 22 A(J,I) = R 23 DO 24 J = K1,N1 24 A(K,J) = A(K,J)/S DO 25 I = K1,N R=A(I,K) DO 25 J = K1,N1 25 A(I,J) = A(I,J) - A(K,J)*R X(N) = A(N,N1) DO 27 I = N-1,1,-1 S1 = A(I,N1) DO 26 J= I+1,N 26 S1 = S1 - A(I,J)*X(J) 27 X(I) = S1 RETURN END