C********************************************************************* C* *****MLSpStr2.for from MLSpStr1.for ****** C Сферически распределенные многослойные структуры C по 4-му варианту с уточнением расстояния между N-м и 1-м телами C Joseph J. Smulsky C Начато 05.09.2017 г. Последнее изменение 25.10.2021 г. C Памятки. C C********************************************************************* C--------------------------------------------------- C Нужно задать значения параметров! PARAMETER (l2=30002) ! C --------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION XV(2,l2),YV(2,l2),ZV(2,l2),XV0(2,l2),YV0(2,l2), *ZV0(2,l2),U1(l2,3),Pns(3),Xc(3),AM1(3),Pn(3),omo(l2), *Xo(l2,3),Uo(l2,3),X(l2,3),U(l2,3),Sp(l2,3),Et(l2), *Ra(l2),dUp(l2,3),X1(l2,3) INTEGER*4 N,N3,Idta,q,lmnn1,lmnn2,Jl CHARACTER*12 FIMA1,FIMA3,FIMA5,FIMA6,FIMA7 FIMA1 ='MlSpStr2.dat' FIMA3 ='MlSpStr2Err' FIMA6 ='fN3fvout.dat' FIMA7 ='fN3fvinp.dat' C--- Считывание файла MlSpStr2.dat. OPEN(UNIT=1,FILE=FIMA1,STATUS='OLD') READ(1,39,ERR=40)N2 READ(1,39,ERR=40)N30 READ(1,54,ERR=40)smi READ(1,54,ERR=40)pm0 READ(1,54,ERR=40)Asm READ(1,54,ERR=40)ec READ(1,54,ERR=40)oka READ(1,54,ERR=40)okN3 READ(1,54,ERR=40)okf0 READ(1,54,ERR=40)okfv READ(1,54,ERR=40)EPS READ(1,39,ERR=40)Inx READ(1,39,ERR=40)Icm 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) 39 FORMAT(59X,I12) 54 FORMAT(59X,D11.5) 55 FORMAT(I4,3(D13.6)) GO TO 1 40 WRITE(*,*)'*01MS Error reading file MlSpStr2.dat *** ', *N2,N30,smi,pm0,Asm,ec,oka,okN3,okf0,okfv,EPS,Inx,eps0,Rob,dt, *FIMA5,Idta OPEN(UNIT=3,FILE=FIMA3) WRITE(3,*)'*01MS Error reading file MlSpStr2.dat *** ', *N2,N30,smi,pm0,Asm,ec,oka,okN3,okf0,okfv,EPS,Inx,eps0,Rob,dt, *FIMA5,Idta CLOSE(3) GO TO 150 1 pi = 3.1415926535897932385D+00 p2 = 2.0D+00*pi k1 = 1 J2 = 1 DZI = 0.0D+00 G = 6.67259D-11 AU = 1.49597870691D+11 Yrs = 365.25636042D+00 dsc = 24D+00*3600D+00 PND1 = 3.0D+00/2.0D+00 PND2 = 1.0D+00/3.0D+00 Rminn = 1.0D+99 C--- Массы центрального тела и одного периферйного тела cm00 = smi*pm0 sm = smi - cm00 db = p2/N2 N3 = N30 C N = N3 + 1 N = 1 Np = 0 cm1 = sm/N3 C cm1p = cm1 cm0 = cm00 ax0 = Asm*AU C ax = ax0 al1 = -1/(1 + ec) OPEN(UNIT=3,FILE=FIMA3) OPEN(UNIT=6,FILE=FIMA6) DO J = 1,N2 okf = okf0 bj = db*(J-1) cb = DCOS(bj) sb = DSIN(bj) C IF(J.EQ.1) GO TO 12 ax = Asm*AU*(1 + oka*(J - 1)) N3p = N3 cm1p = cm1 N3 = N30*okN3*(ax/ax0) OPEN(UNIT=7,FILE=FIMA7,STATUS='OLD',ERR=20) DO J1 = 1,J READ(7,*,ERR=20)J2,N3,okf,okfv,ax,cm1 END DO ! J1 CLOSE(7) 20 IF(J.NE.1) Np = N N = N + N3 IF(J.NE.1) cm0 = cm0 + cm1p*N3p 12 Rmin = 1.0D+99 okf2 = DZI YNw2 = DZI dfi0 = p2/N3 fn3 = DZI DO l = 2,N3 fn3 = fn3 + 0.25D+0/DSIN((l-1)*pi/N3) END DO ! l C--- Период обращения Rp = ax*(2*al1 + 1)/al1 omu1 = -G*(cm0 + cm1*fn3) Vp = DSQRT(omu1/(al1*Rp)) TMP = Vp*(-2*al1-1)**PND1 Prdsc = -p2*al1*Rp/TMP IF(J.Eq.1) Prdsc1 = Prdsc Prd = Prdsc/(Yrs*dsc) C--- Время в сотнях периодов Prd, который измеряется C--- в сидерических годах DO I1 = 1,20 I1m = I1 dfi1 = okf*dfi0 dfi2 = okfv*dfi0 dfi3 = dfi1*2.0D+0 C1 = DCOS(dfi1) S1 = DSIN(dfi1) C2 = DCOS(dfi2) S2= DSIN(dfi2) C3 = DCOS(dfi3) S3 = DSIN(dfi3) C Трехмерные структуры C Cтруктура по варианту 4 C l, I - номера тел кольца C--- Параметры тел на своих орбитах 6 DO I = 1,N3 + 1 fi = dfi0*(I-1) CI = DCOS(fi) SI = DSIN(fi) ri = Rp/((al1+1)*CI - al1) Vr = Vp*DSQRT((al1+1)*(al1+1)-(al1+Rp/ri)*(al1+Rp/ri)) IF(SI.LT.DZI) vr = -vr Vt = Vp*Rp/ri XV(1,I) = ri*CI YV(1,I) = ri*SI ZV(1,I) = DZI XV(2,I) = Vr*CI - Vt*SI YV(2,I) = Vr*SI + Vt*CI ZV(2,I) = DZI END DO ! I C Поворот векторов скорости тел кольца DO l3 = 1,2 DO l1 = 2,N3+1 l22 = N3 + 3 - l1 DO l = 2,l22 Xt = XV(l3,l)*C1 - YV(l3,l)*C2*S1 + ZV(l3,l)*S1*S2 Yt = XV(l3,l)*S1 + YV(l3,l)*C1*C2 - ZV(l3,l)*C1*S2 Zt = YV(l3,l)*S2 + ZV(l3,l)*C2 XV(l3,l-1) = Xt YV(l3,l-1) = Yt ZV(l3,l-1) = Zt END DO ! l XV0(l3,l1+1) = XV(l3,1) YV0(l3,l1+1) = YV(l3,1) ZV0(l3,l1+1) = ZV(l3,1) END DO ! l22 END DO ! l3 C l- номера цетрального тела и тел кольца XV0(1,1) = DZI YV0(1,1) = DZI ZV0(1,1) = DZI XV0(2,1) = DZI YV0(2,1) = DZI ZV0(2,1) = DZI XV0(1,2) = Rp YV0(1,2) = DZI ZV0(1,2) = DZI XV0(2,2) = DZI YV0(2,2) = Vp ZV0(2,2) = DZI C Уточнение расстояния между N-м и 1-м телами okf1 = okf2 okf2 = okf YNw1 = YNw2 YNw2 = YV0(1,N3+2) IF(I1.NE.1) GO TO 9 okf = 1.05*okf GO TO 10 9 IF(okf2.NE.DZI) THEN EPSI = DABS(okf2-okf1)/okf2 ELSE WRITE(*,*)'Error okf2 = ',okf2 GO TO 150 END IF IF(EPSI.LE.EPS) GO TO 8 okf = okf2 - YNw2*(okf2 - okf1)/(YNw2 - YNw1) 10 CONTINUE END DO ! I1 8 DO k = 1, N3+1 IF(k.EQ.N3+1.AND.J.NE.1) GO TO 14 k2 = Np + k IF(J.EQ.1) THEN k3 = k ELSE k3 = k + 1 END IF omo(k2) = cm1 Xo(k2,1) = XV0(1,k3) Xo(k2,2) = YV0(1,k3) Xo(k2,3) = ZV0(1,k3) Uo(k2,1) = XV0(2,k3) Uo(k2,2) = YV0(2,k3) Uo(k2,3) = ZV0(2,k3) IF(Inx.EQ.0) GO TO 14 IF(J.EQ.1) GO TO 14 Tmpy = Xo(k2,2)*cb + Xo(k2,3)*sb Tmpz = -Xo(k2,2)*sb + Xo(k2,3)*cb Xo(k2,2) = Tmpy Xo(k2,3) = Tmpz Tmpy = Uo(k2,2)*cb + Uo(k2,3)*sb Tmpz = -Uo(k2,2)*sb + Uo(k2,3)*cb Uo(k2,2) = Tmpy Uo(k2,3) = Tmpz 14 CONTINUE END DO ! k C--- Minimal distance between bodies lmn1 and lmn2 DO l = 2,N3 + 1 DO l1 = 2,N3 + 1 IF(l1.EQ.l) GO TO 4 XX = XV0(1,l1) - XV0(1,l) XX = XX*XX YY = YV0(1,l1) - YV0(1,l) YY = YY*YY ZZ = ZV0(1,l1) - ZV0(1,l) ZZ = ZZ*ZZ RIJ = DSQRT(XX + YY + ZZ) IF(l.EQ.2.AND.l1.EQ.3) R23 = RIJ IF(RIJ.GT.Rmin) GO TO 4 Rmin = RIJ lmn1 = l lmn2 = l1 4 END DO ! l1 END DO ! l IF(Rmin.GT.Rminn) GO TO 5 Rminn = Rmin lmnn1 = lmn1 lmnn2 = lmn2 Jl = J 5 WRITE(3,*)'****** Parameters of ',J,' ring ******' IF(I1m.EQ.20) WRITE(3,*)'On 20th step error EPS is no received' IF(Icm.EQ.0) GO TO 18 DO l = 1,N3 + 2 DO l3 = 1,2 WRITE(3,55)l,XV0(l3,l),YV0(l3,l),ZV0(l3,l) END DO ! l3 END DO ! l C WRITE(3,*)l,XV0(1,2),YV0(1,2),ZV0(1,2) C WRITE(3,*)l,XV0(1,100),YV0(1,100),ZV0(1,100) 18 WRITE(3,*)'Period in sideral years = ',Prd,' and N3 = ',N3 WRITE(3,*)'Distance between bodies 2 and 3',R23 WRITE(3,*)'Minimal distance between bodies ',lmn1, *' and ',lmn2,' ',Rmin WRITE(3,*)'Last step of itrrtn I1m',I1m,' and coef-ts okf,okf1', *okf,okf1 WRITE(6,*)J,N3,okf,okfv,ax,cm1 END DO ! J CLOSE(6) C ssm = cm00 + cm1*(N - 1) ssm = cm0 + cm1*N3 okt = 1/(100*Prdsc1) C--- Массы тел в программе Galactica и начальные параметры omo(1) = cm00/ssm TMP = G*ssm/(okt*okt) Am = TMP**PND2 okv = 1/(okt*Am) Robo = Rob*Am*Am*Am/ssm DO k = 1,N IF(k.NE.1) omo(k) = omo(k)/ssm Xo(k,1) = Xo(k,1)/Am Xo(k,2) = Xo(k,2)/Am Xo(k,3) = Xo(k,3)/Am Uo(k,1) = Uo(k,1)*okv Uo(k,2) = Uo(k,2)*okv Uo(k,3) = Uo(k,3)*okv END DO ! k 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) 19 Ra(k) = TMP**PND2 IF(omo(k).LT.omm) GO TO 13 omm = omo(k) 13 END DO ! k A = 1.0D+0 B = 0.5D+0 C = 0.25D+0 TMP = (N+1)**PND2 Mu = NINT(TMP) T = DZI Um = DZI Spsx = DZI Spsy = DZI Spsz = DZI E = DZI Em = DZI Ett = DZI dtp = DZI WRITE(3,*)'****** Parameters of all structure ******' WRITE(3,*)'The minimal distance between bodies ',lmnn1, *' and ',lmnn2,' Rminn = ',Rminn, ' is in ',Jl,' ring.' WRITE(3,*)'Coordinates and velocities of the system baricentre' WRITE(3,*)(Xc(q),q=1,3) WRITE(3,*)(Pns(q),q=1,3) WRITE(3,*)'Initial undimensional coordinates and velocities' WRITE(3,*)' of the 1-st body' WRITE(3,*)(X1(2,q),q=1,3) WRITE(3,*)(U1(2,q),q=1,3) WRITE(3,*)'Initial data in file ',FIMA1 WRITE(3,*)N2,N30,smi,pm0,Asm,ec,oka,okN3,okf0,okfv,EPS,Inx,Icm, *eps0,Rob,dt,FIMA5,Idta CLOSE(3) 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