C********************************************************************* C* *****SphDsSt4.for from SphDsStr3.for ****** C Сферически распределенные однослойные структуры C Joseph J. Smulsky C Начато 28.11.2015 г. Последнее изменение 22.07.2016 г. C********************************************************************* C--------------------------------------------------- C Нужно задать значения параметров! PARAMETER (l2=1002,l4=6000) ! 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 CHARACTER*12 FIMA1,FIMA3,FIMA5,FIMA6 C CHARACTER*59 Space59 FIMA1 ='SphDsSt4.dat' FIMA3 ='SphDsSt4Err' FIMA6 ='file_okf.dat' C--- Считывание файла SphDsSt3.dat. OPEN(UNIT=1,FILE=FIMA1,STATUS='OLD') READ(1,39,ERR=40)N3 READ(1,54,ERR=40)smi READ(1,54,ERR=40)pm0 READ(1,54,ERR=40)Asm READ(1,54,ERR=40)ec C READ(1,54,ERR=40)Space59,okf READ(1,54,ERR=40)okf READ(1,54,ERR=40)okfv READ(1,54,ERR=40)eps0 READ(1,54,ERR=40)Rob READ(1,39,ERR=40)Ivr READ(1,54,ERR=40)dt READ(1,'(59X,A12)',ERR=40)FIMA5 READ(1,39,ERR=40)Idta CLOSE(1) OPEN(UNIT=1,FILE=FIMA6,STATUS='OLD',ERR=1) READ(1,*,ERR=1)okf CLOSE(1) 39 FORMAT(59X,I12) 54 FORMAT(59X,D11.5) 55 FORMAT(I4,3(D13.6)) GO TO 1 40 WRITE(*,*)'*01Rt Error reading file SphDsSt3.dat *** ', *N3,smi,pm0,Asm,Prd,ec,okf,okfv,eps0,Rob,dt,FIMA5,Idta OPEN(UNIT=3,FILE=FIMA3) WRITE(3,*)'*01Rt Error reading file SphDsSt3.dat *** ', *N3,smi,pm0,Asm,Prd,ec,okf,okfv,eps0,Rob,dt,FIMA5,Idta CLOSE(3) GO TO 150 1 pi = 3.1415926535897932385D+00 p2 = 2.0D+00*pi G = 6.67259D-11 AU = 1.49597870691D+11 Yrs = 365.25636042D+00 dsc = 24D+00*3600D+00 Rmin = 1.0D+99 DZI = 0.0D+00 PND1 = 3.0D+00/2.0D+00 PND2 = 1.0D+00/3.0D+00 dfi0 = p2/N3 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) fn3 = DZI DO l = 2,N3 fn3 = fn3 + 0.25D+0/DSIN((l-1)*pi/N3) END DO ! l N = N3 + 1 k1 = 1 C--- Массы центрального тела и одного периферйного тела cm0 = smi*pm0 sm = smi - cm0 cm1 = sm/N3 ax = Asm*AU al1 = -1/(1 + ec) 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 Prd = Prdsc/(Yrs*dsc) okt = 1/(100*Prdsc) C--- Время в сотнях периодов Prd, который измеряется C--- в сидерических годах C Трехмерные структуры IF(Ivr.EQ.4) GO TO 6 C Cтруктура по варианту 3 C l- номера тел кольца IF(N3.LT.3) GO TO 2 DO l = 2,N3 l5 = (-1)**(l+1) XV(1,l) = Rp*C3 YV(1,l) = Rp*S3 ZV(1,l) = DZI XV(2,l) = -Vp*S3*C2 YV(2,l) = Vp*C3*C2 ZV(2,l) = Vp*S2*l5 C WRITE(*,*)'начальн',l5,l,ZV(2,l) END DO ! l C Поворот векторов скорости нечетных тел кольца DO l3 = 1,2 DO l1 = 3,N3+1,2 l22 = N3 + 3 - l1 DO l = 3,l22,2 l5 = (-1)**l Xt = XV(l3,l)*C3 - YV(l3,l)*C2*S3 + ZV(l3,l)*S3*S2*l5 Yt = XV(l3,l)*S3 + YV(l3,l)*C3*C2 - ZV(l3,l)*C3*S2*l5 Zt = YV(l3,l)*S2*l5 + ZV(l3,l)*C2 XV(l3,l-2) = Xt YV(l3,l-2) = Yt ZV(l3,l-2) = 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 ! l1 END DO ! l3 IF(N3.LT.4) GO TO 2 C Поворот векторов скорости четных тел кольца Vp1 = -Vp DO l = 2,N3 l5 = (-1)**(l+1) XV(1,l) = Rp*C3 YV(1,l) = Rp*S3 ZV(1,l) = DZI XV(2,l) = -Vp1*S3*C2 YV(2,l) = Vp1*C3*C2 ZV(2,l) = Vp1*S2*l5 END DO ! l DO l3 = 1,2 DO l1 = 4,N3+1,2 l22 = N3 + 4 - l1 DO l = 4,l22,2 Xt = XV(l3,l)*C3 - YV(l3,l)*C2*S3 + ZV(l3,l)*S3*S2 Yt = XV(l3,l)*S3 + YV(l3,l)*C3*C2 - ZV(l3,l)*C3*S2 Zt = YV(l3,l)*S2 + ZV(l3,l)*C2 XV(l3,l-2) = Xt YV(l3,l-2) = Yt ZV(l3,l-2) = Zt END DO ! l XV0(l3,l1+1) = XV(l3,2) YV0(l3,l1+1) = YV(l3,2) ZV0(l3,l1+1) = ZV(l3,2) C WRITE(*,*)'четн',l3,l22,l1,ZV0(l3,l1+1) END DO ! l1 END DO ! l3 C l- номера цетрального тела и тел кольца 2 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 XV0(1,3) = Rp*C3 YV0(1,3) = Rp*S3 ZV0(1,3) = DZI XV0(2,3) = -Vp*S3*C2 YV0(2,3) = Vp*C3*C2 ZV0(2,3) = Vp*S2 GO TO 8 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 8 ssm = smi C--- Массы тел в программе Galactica и начальные параметры omo(1) = cm0/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) = cm1/ssm Xo(k,1) = XV0(1,k)/Am Xo(k,2) = YV0(1,k)/Am Xo(k,3) = ZV0(1,k)/Am Uo(k,1) = XV0(2,k)*okv Uo(k,2) = YV0(2,k)*okv Uo(k,3) = ZV0(2,k)*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 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) 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 OPEN(UNIT=3,FILE=FIMA3) 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 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 period in sideral years = ',Prd WRITE(3,*)'Distance between bodies 2 and 3',R23 WRITE(3,*)'Minimal distance between bodies ',lmn1, *' and ',lmn2,' ',Rmin WRITE(3,*)'Initial data in file ',FIMA1 WRITE(3,*)N3,smi,pm0,Asm,ec,okf,okfv,eps0,Rob,Ivr,dt,' ',FIMA5, *Idta CLOSE(3) 150 CONTINUE END