c ssbox_rinit_eh_2007.f c test if anything can be speeded up C***************************************************************** SUBROUTINE RINIT C***************************************************************** INCLUDE 'ssbox.inc' DIMENSION SV(NSTAR),INDX(NSTAR) DOUBLE PRECISION MYYV(NSTAR),MKES(10) REAL XPRO(NSTAR),YPRO(NSTAR) DIMENSION CX1(10),CX2(10),CI(10),CE(10),NC(10),TAU(10) DOUBLE PRECISION INK,INK2,INKB double precision cifac(10),cefac(10) real sran,dummy,ran1 real ampfidl(4000),phasefidl(4000) dimension induse(nstar) TIME=0.D0 time0=time ENES=0.D0 ENEU=0.D0 enes_coll=0.d0 eneu_coll=0.d0 DSIG(0)=0. c ssbox_c1.f: c add initial seed for perturbations in vx c amounting to afidl*cos(m*(xx/xl*!pi-phase)) c phase is randomly chosen c------------------------------------------------------------ c read new parameters (FROM IIOCT 141203) c note: tamp,eamp not used c tamp = relative amplitude of T-perturbation c vamp = amplitude of V-perturbation (units sigma*1m) c namp = number of zones used in discretization of perturbation c defines also number of collection zones c interval of storing is the same as for the vertical c distribution c new-method: c eamp = relative amplitude of eps-perturbation c nov99 c nwave = number of wavelengths fitting to box (previously 1) c _c4 Kprof-versio = 5 c start_amp,dstart_amp -> recurrent perturbation at c ORB= start_amp+dstart_amp*i write(*,*) 'tamp,vamp,eamp,namp,nwave' read(*,*) tamp,vamp,eamp,namp,nwave write(*,*) tamp,vamp,eamp,namp,nwave c eboxnov99_c4.for: write(*,*) 'start_amp,dstart_amp' read(*,*) start_amp,dstart_amp write(*,*) start_amp,dstart_amp c new write(*,*) 'dorbk1,dorbk2' read(*,*) dorbk1,dorbk2 write(*,*) dorbk1,dorbk2 C***************************************************************** C READ INPUT PARAMETERS write(idev6,*) '**********************************' write(idev6,*) 'SSBOX_EM-SIMULATION with kprof-tab' write(idev6,*) '**********************************' write(idev7,*) '**********************************' write(idev7,*) 'SSBOX_EM-SIMULATION with kprof-tab' write(idev7,*) '**********************************' WRITE(IDEV6,*) ' STARTUP (0) OR CONTINUATION (1)' READ(IDEV5,*) IJATKO WRITE(IDEV6,*) ' norbit, orbtest' READ(IDEV5,*) norbit,orbtest WRITE(IDEV6,*) ' MIX (1) OR NOMIX (0)' READ(IDEV5,*) IMIX WRITE(IDEV6,*) ' LKM,ETA,SEED' READ(IDEV5,*) LKM,ETA,ISEED WRITE(IDEV6,*) ' KAPPA/OMEGA, WZ/OMEGA' READ(IDEV5,*) KAPPAFAC,WZFAC W0=DSQRT(GM/ETA/ETA/ETA) OMEGA=W0 W2=2.D0*OMEGA KAPPA=KAPPAFAC*OMEGA OMEGA2=OMEGA*OMEGA KAPPA2=KAPPA*KAPPA WZ=W0*WZFAC WZ2=WZ*WZ PERIOD=2.D0*PII/OMEGA WRITE(IDEV6,*) ' NUMBER OF RAD.GROUPS ' READ(IDEV5,*) LKMR LKMC=0 CI2=0. CE2=0. DO 77 I=1,LKMR WRITE(IDEV6,*) 'RAD.GROUP: NC,CXMIN,CXMAX,CINK,CEKS' WRITE(IDEV7,*) 'RAD.GROUP: NC,CXMIN,CXMAX,CINK,CEKS' READ(IDEV5,*) NC(I),CX1(I),CX2(I),CI(I),CE(I) WRITE(IDEV6,*) NC(I),CX1(I),CX2(I),CI(I),CE(I) WRITE(IDEV7,*) NC(I),CX1(I),CX2(I),CI(I),CE(I) LKMC=LKMC+NC(I) CI2=CI2+NC(I)*CI(I) CE2=CE2+NC(I)*CE(I) 77 CONTINUE IF(LKMC.NE.LKM) STOP CINK=CI2/LKM CEKS=CE2/LKM WRITE(IDEV6,*) 'BOX SIZE: XX,YY ' WRITE(IDEV6,*) 'IF XX=0 TAU=YY' READ(IDEV5,*) XL,YL WRITE(IDEV6,*) 'LEV,SIG4,CYMAX(IF>0)' READ(IDEV5,*) LEV,SIG4,DEF CXMAX=(CX2(LKMR)-CX1(1))/2. CYMAX=YL IF(DEF.GT.0) CYMAX=DEF C INITIALIZE NAG-RANDOM NUMBER GENERATOR dummy=ran1(iseed) iidum=0 C********************************************************************* C MASS DISTRIBUTION C********************************************************************* WRITE(IDEV6,*) 'DISCRETE (1) OR CONTINUOUS DISTRIBUTION (2)' READ(IDEV5,*) VAL C **************************************** C DISKREETTI IF(VAL.EQ.1) THEN MJAK=1 WRITE(IDEV6,*) 'NUMBER OF SIZE GROUPS' READ(IDEV5,*) NSIG WRITE(IDEV6,*) ' SIG, NUMBER, CIFAC, CEFAC' IND=0 DO 5 I=1,NSIG READ(IDEV5,*) SIG,NUM(I),cifac(i),cefac(i) IF(I.EQ.1) SIGA=SIG DSIG(I)=1.0001*SIG DO 6 J=1,NUM(I) IND=IND+1 S(IND)=SIG xmass(ind)=4.*pii/3.*roo*s(ind)**3 IS(IND)=I c------------------------------------- c MYY(IND)=(SIG/SIGA)**3. c 14/01/03 myy(ind)=xmass(ind) if(roo.le.0) + MYY(IND)=4.*pii/3.*900.*s(ind)**3 c------------------------------------- 6 CONTINUE 5 CONTINUE WRITE(IDEV6,*) ' NTOT ',IND LKM=IND ENDIF C **************************************** C CONTINUOUS DISTRIBUTION IF(VAL.EQ.2) THEN MJAK=2 WRITE(IDEV6,*) ' NUMBER OF SIZE-BINS' READ(IDEV5,*) NSIG WRITE(IDEV6,*) ' Q,SIGA,SIGB' READ(IDEV5,*) Q,SIGA,SIGB QQ=Q VAK= DLOG(SIGB/SIGA)/NSIG DO 10 I=1,NSIG DSIG(I)=DEXP(VAK*I)*SIGA*1.00001 10 CONTINUE EXPO=1-Q S1=SIGA**EXPO S2=SIGB**EXPO DO 19 I=1,LKM SS=(1.*(I-1)/(LKM-1)*(S2-S1)+S1)**(1./EXPO) MYY(I)=(SS/SIGA)**3 S(I)=SS xmass(i)=4.*pii/3.*roo*s(i)**3 c------------------------------------- c MYY(IND)=(SIG/SIGA)**3. c 14/01/03 myy(i)=xmass(i) if(roo.le.0) + MYY(I)=4.*pii/3.*900.*s(i)**3 c------------------------------------- 19 continue write(idev6,*) 'smin,smax',s(1),s(lkm) write(idev7,*) 'smin,smax',s(1),s(lkm) ENDIF C********************************************************************* C CALCULATE MEAN VALUES FOR EACH SIZE GROUP C CREATE IS-TABLES C IF IMIX=1 MIX PARTICLE INDEXES BEFORE DOING THIS C********************************************************************* IF(IMIX.EQ.1) THEN DO 290 I=1,LKM-1,2 XPRO(I)=ran1(iidum) XPRO(I+1)=XPRO(I) 290 CONTINUE CALL INDEXR(LKM,XPRO,INDX) DO 291 I=1,LKM SV(I)=S(I) MYYV(I)=MYY(I) 291 CONTINUE DO 292 I=1,LKM S(I)=SV(INDX(I)) MYY(I)=MYYV(INDX(I)) 292 CONTINUE IF(INDLARGE.GT.0) INDLARGE=INDX(INDLARGE) ENDIF WRITE(IDEV6,*) 'DSIG' WRITE(IDEV6,*) (DSIG(I),I=0,NSIG) C ******* IF XL=0 USE YL AS TAU IF(XL.EQ.0) THEN PA=0. DO 4567 I=1,LKM PA=PA+PII*S(I)*S(I) 4567 CONTINUE XL=SQRT(PA/YL/4.) YL=XL CYMAX=0.9999*YL CX1(1)=-CYMAX CX2(1)=CYMAX WRITE(6,*) 'XL,YL -> ' WRITE(6,*) 'TAU=',PA/4./XL/YL WRITE(6,*) 'XL=YL=',XL WRITE(6,*) 'CXMAX=CYMAX=',CYMAX ENDIF TAU0=0. ntest=0 mtot_size(0)=0. DO 29 I=1,NSIG NUM(I)=0 TAU(I)=0. MKES(I)=0. mtot_size(i)=0. 29 continue DO 30 I=1,LKM IS(I)=NSIG DO 31 J=1,NSIG IF((S(I).LE.DSIG(J)).AND.(S(I).GT.DSIG(J-1))) + THEN IS(I)=J NUM(J)=NUM(J)+1 c ******* test particle indexes (norbit/group) if(num(j).le.norbit) then ntest=ntest+1 iorbit(ntest)=i isorbit(ntest)=is(i) endif c********************************************** TAU(J)=TAU(J)+PII*S(I)*S(I)/CXMAX/YL/4. TAU0=TAU0+PII*S(I)*S(I)/CXMAX/YL/4. MKES(J)=MKES(J)+MYY(I) c 170403: total masses in size bins mtot_size(j)=mtot_size(j)+myy(i) mtot_size(0)=mtot_size(0)+myy(i) GO TO 30 ENDIF 31 CONTINUE 30 CONTINUE c-------------------------------------------------------------- c 170403: write to viscosize-file write(47,*) nsig write(47,*) xl,yl,lkm write(47,*) omega,kappa,wz do 4747 i=0,nsig if(i.ne.0) apu=num(i) if(i.eq.0) apu=lkm write(47,*) i,dsig(i),apu,mtot_size(i) 4747 continue c-------------------------------------------------------------- c total number of test particles nnorbit=ntest c alternative way of assigning index for test-particles c norbit=neg -> READ FROM FILE 'orbit.index' (fixed name) if(norbit.lt.0) then ntest=0 write(*,*) 'test-orbit indexes from "orbit.index"' open(unit=34,file='orbit.index',status='old') do 3434 i=1,-norbit read(34,*) iorbit(i) write(*,*) 'i,iorbit',i,iorbit(i) isorbit(i)=-1 3434 continue close(34) nnorbit=-norbit endif c minimum optical thickness if particles fill whole box TAUMIN=TAU0*CXMAX/XL C ******* CONSTANTS RELATED TO IMAGE-PARTICLES C X -> X + NN*AVAK C Y -> Y + NN*BVAK +MM*CVAK C VX -> VX C VY -> VY + NN*VVAK C ******* EQUATIONS OF MOTION C XPP - 2*OMEGA * YP + FACTEMP * X = FX C YPP + 2*OMEGA * XP = FY C ZPP + WZ2 * Z = FZ FACTEMP=KAPPA2-4.D0*OMEGA2 SHEAR=-FACTEMP/2.d0/OMEGA AVAK=2.D0*XL CVAK=2.D0*YL VVAK=FACTEMP/OMEGA*XL c continuation simulation -> modified later BVAK=VVAK*TIME WRITE(IDEV6,*) ' SSBOX_EM-SIMULATION ' WRITE(IDEV6,*) ' DISTANCE FROM PLANET (10E3 KM)',ETA/1.D6 WRITE(IDEV6,*) ' KAPPA, WZ',KAPPAFAC,WZFAC WRITE(IDEV6,*) ' ORIGINAL CXMAX CYMAX',CXMAX,CYMAX WRITE(IDEV6,*) ' ORIGINAL INK2 EKS2',CINK,CEKS WRITE(IDEV6,*) ' BOX-SIZE XL YL ',XL,YL WRITE(IDEV6,*) ' ' WRITE(IDEV6,*) 'GROUP SIG1-SIG2 MA-MB NN MKES TAU' WRITE(IDEV6,*) ' ' WRITE(7,*) ' ' WRITE(7,*) 'GROUP SIG1-SIG2 MA-MB NN MKES TAU' WRITE(7,*) ' ' DO 40 I=1,NSIG SIG1=DSIG(I-1) SIG2=DSIG(I) IF(NUM(I).NE.0) MKES(I)=MKES(I)/NUM(I) MA=(SIG1/SIGA)**3 MB=(SIG2/SIGA)**3 WRITE(7,41) I,SIG1,SIG2,MA,MB,NUM(I),MKES(I),TAU(I) WRITE(IDEV6,41) I,SIG1,SIG2,MA,MB,NUM(I),MKES(I),TAU(I) 40 continue 41 FORMAT(I6,2F8.3,2F10.3,I8,F10.3,F12.6) C--------MARCH 2002/modified JAN 2003 ------------------- c embedded moonlet: PLACE THIS AS PARTICLE 1 C (OVERWRITE THE PREVIOUS PARTICLE) c cen_moon= 1, moonlet fixed - dimpacts subtracts accel c cen_moon= 0 , moonlet free to move c always use the actual mass for collisional calculations c give zero mass for moonlet gravitational xmass, since c the force due to moonlet is added separately c group index of moonlet set to is(1)=0 c in the end of rinit myyspin=0. ind_moon=0 if(n_moon.ne.0) then ind_moon=1 s(1)=rad_moon xrho_moon=rho_moon if(rho_moon.le.0) xrho_moon=900. myy(1)=4.*pii/3.*xrho_moon*rad_moon**3 xmass(1)=0. myyspin=myy(1) endif c--------------------------------------------------- c initialize kprof-collection c added from iicot 141203 c--------------------------------------------------- drzone=2.*xl/namp do 740 i=1,namp rampzones1(i)=-xl+(i-1)*drzone rampzones2(i)=-xl+i*drzone rampzone=-xl+(i-0.5)*drzone ampzones(i)=1.+amp*cos(rampzone/xl*pii*nwave) 740 CONTINUE KVERSIO=5 KVERSIO=6 WRITE(41) KVERSIO write(41) lkm WRITE(41) NAMP,tAMP,vamp,eamp,nwave WRITE(41) XL,YL c kversio=6 -> output also kappa c write(41) w0,wz write(41) w0,wz,kappa WRITE(41) (RAMPZONES1(i),i=1,namp) WRITE(41) (RAMPZONES2(i),i=1,namp) WRITE(41) (AMPZONES(i),i=1,namp) write(41) start_amp,dstart_amp C************************************************************** C CREATE RADIUS AND VELOCITY VECTORS C SKIP IF CONTINUATION RUN C************************************************************** IF(IJATKO.EQ.1) GOTO 101 c------------------------------------------------------------ c 270107 c let's at last change it so that big particles are created first if(n_moon.eq.0) then do 2007 i=1,lkm induse(i)=lkm+1-i 2007 continue endif if(n_moon.ne.0) then induse(1)=1 induse(2)=2 do 2008 i=3,lkm induse(i)=lkm+3-i 2008 continue endif c------------------------------------------------------------ I_induse=0 DO 1111 NR=1,LKMR LASK=0 INK2=0.D0 EKS2=0.D0 CXMI=CX1(NR) CXMA=CX2(NR) c added 13.2.98 to assure uniform initial distribution if(cxmi.lt.-xl) cxmi=-xl if(cxma.gt. xl) cxma= xl CINK=CI(NR) CEKS=CE(NR) write(6,*) nr,cxmi,cxma,cink,ceks DO 111 IR=1,NC(NR)/2 I_induse=I_induse+2 i=induse(i_induse) im1=induse(i_induse-1) cifac0=1. cefac0=1. if(val.eq.1) then cifac0=cifac(is(i)) cefac0=cefac(is(i)) endif 112 CONTINUE C ******* CREATE PARTICLES IN PAIRS SO THAT C UVAK=SUMM (MYY XP) = 0 C WVAK=SUMM (MYY(YP+1.5OMEGAX)) = 0 C THIS IS TRUE IF FOR EACH PAIR C MYY1 EKS1 = MYY2 EKS2 C DELTA1 = DELTA2 + PII c RAYLEIGH DISTRIBUTION OF EKS,INC (see TRULSEN 1972) c first particle of the pair index=i Sran=ran1(iidum) CX=CXMI+Sran*(CXMA-CXMI) Sran=ran1(iidum) CY=-CYMAX+2.*Sran*CYMAX Sran=ran1(iidum) DELTA=2.D0*PII*Sran Sran=ran1(iidum) GAMMA=2.D0*PII*Sran Sran=ran1(iidum) INK=cifac0*CINK*DSQRT(DLOG(1.D0/(1.D0-Sran))) Sran=ran1(iidum) EKS=cefac0*CEKS*DSQRT(DLOG(1.D0/(1.D0-Sran))) c if moonlet is included then set eks, ink of its pair to zero c to have zero velocities c c IF(I-1.EQ.IND_moon) THEN c EKS=0. c INK=0. c ENDIF R(I,1)=CX+EKS*DCOS(DELTA) R(I,2)=CY-2.D0*omega/kappa*EKS*DSIN(DELTA) R(I,3)=INK*DSIN(GAMMA) RP(I,1)=-kappa*EKS*DSIN(DELTA) RP(I,2)=-2.D0*W0*EKS*DCOS(DELTA)+ + fh3*CX*(kappa2-4.d0*omega2)/2.d0/omega RP(I,3)=INK*WZ*DCOS(GAMMA) c second particle of the pair index=i-1 Sran=ran1(iidum) CXB=CXMI+Sran*(CXMA-CXMI) Sran=ran1(iidum) CYB=-CYMAX+2.*Sran*CYMAX EKSB=MYY(I)/MYY(Im1)*EKS INKB=MYY(I)/MYY(Im1)*INK R(Im1,1)=CXB+EKSB*DCOS(DELTA+PII) R(Im1,2)=CYB-2.D0*omega/kappa*EKSB*DSIN(DELTA+PII) R(Im1,3)=INKB*DSIN(GAMMA+PII) RP(Im1,1)=-kappa*EKSB*DSIN(DELTA+PII) RP(Im1,2)=-2.D0*W0*EKSB*DCOS(DELTA+PII)+ + fh3*CXb*(kappa2-4.d0*omega2)/2.d0/omega RP(Im1,3)=INKB*WZ*DCOS(GAMMA+PII) c if it is a moonlet then set to origin with c zero velocities c replace by moonlet epicycle 170605 IF((Im1).EQ.ind_moon) THEN R(Im1,1)=CX_moon+EKS_moon*DCOS(DELTA_moon+PII) R(Im1,2)=CY_moon- + 2.D0*omega/kappa*EKS_moon*DSIN(DELTA_moon+PII) R(Im1,3)=INK_moon*DSIN(GAMMA_moon+PII) RP(Im1,1)=-kappa*EKS_moon*DSIN(DELTA_moon+PII) RP(Im1,2)=-2.D0*W0*EKS_moon*DCOS(DELTA_moon+PII)+ + fh3*CX_moon*(kappa2-4.d0*omega2)/2.d0/omega RP(Im1,3)=INK_moon*WZ*DCOS(GAMMA_moon+PII) write(*,*) 'moon - in rinit:' write(*,*) 'moon', + CX_moon,CY_moon,eks_moon,ink_moon, + delta_moon,gamma_moon write(*,*) 'moon-R: ',R(im1,1),R(im1,2),R(im1,3) write(*,*) 'moon-V: ',Rp(im1,1),Rp(im1,2),Rp(im1,3) ENDIF C ************************************** C MOVE PAIR INTO CALCULATION BOX N=0 M=0 IF(R(I,1).GE.0) N=-INT(R(I,1)/AVAK+0.5) IF(R(I,1).LT.0) N=INT(-R(I,1)/AVAK+0.5) R(I,1)=R(I,1)+N*AVAK IF(R(I,2).GE.0) M=-INT(R(I,2)/CVAK+0.5) IF(R(I,2).LT.0) M=INT(-R(I,2)/CVAK+0.5) R(I,2)=R(I,2)+M*CVAK RP(I,2)=RP(I,2)+N*VVAK N=0 M=0 IF(R(Im1,1).GE.0) N=-INT(R(Im1,1)/AVAK+0.5) IF(R(Im1,1).LT.0) N=INT(-R(Im1,1)/AVAK+0.5) R(Im1,1)=R(Im1,1)+N*AVAK IF(R(Im1,2).GE.0) M=-INT(R(Im1,2)/CVAK+0.5) IF(R(Im1,2).LT.0) M=INT(-R(Im1,2)/CVAK+0.5) R(Im1,2)=R(Im1,2)+M*CVAK RP(Im1,2)=RP(Im1,2)+N*VVAK LASK=LASK+1 IF(LASK.GT.NC(NR)) THEN LASK=1 CINK=CINK*1.2 WRITE(IDEV6,*) 'CINK',CINK ENDIF C ************************************** C CHECK FOR OVERLAPPING PARTICLES C (ALSO OVER BORDERS) DO 800 J_induse=1,I_induse-1 j=induse(j_induse) EROX=R(I,1)-R(J,1) EROY=R(I,2)-R(J,2) EROZ=R(I,3)-R(J,3) IF(DABS(EROY).GT.YL) EROY=2.D0*YL-DABS(EROY) IF(DABS(EROX).GT.XL) EROX=2.D0*XL-DABS(EROX) ERO=EROX*EROX+EROY*EROY+EROZ*EROZ SKER=(S(I)+S(J))**2 IF(ERO.LT.SKER) GOTO 112 IF(J_induse.EQ.(I_induse-1)) GOTO 800 EROX=R(Im1,1)-R(J,1) EROY=R(Im1,2)-R(J,2) EROZ=R(Im1,3)-R(J,3) IF(DABS(EROY).GT.YL) EROY=2.D0*YL-DABS(EROY) IF(DABS(EROX).GT.XL) EROX=2.D0*XL-DABS(EROX) ERO=EROX*EROX+EROY*EROY+EROZ*EROZ SKER=(S(Im1)+S(J))**2 IF(ERO.LT.SKER) GOTO 112 800 CONTINUE LASK=LASK-1 C ******* PAIR WAS ACCEPTED INK2=INK2+INK*INK+INKB*INKB EKS2=EKS2+EKS*EKS+EKSB*EKSB XPRO(I)=CX YPRO(I)=CY XPRO(Im1)=CXB YPRO(Im1)=CYB IF(MOD(I,50).EQ.0) WRITE(IDEV6,*) I 111 CONTINUE C ******* OUTPUT MEAN VALUES FOR EACH RADIAL GROUP WRITE(IDEV6,*) 'SQRT(INK2) =',SQRT(INK2/NC(NR)) WRITE(IDEV6,*) 'SQRT(EKS2) =',SQRT(EKS2/NC(NR)) WRITE(idev7,*) 'SQRT(INK2) =',SQRT(INK2/NC(NR)) WRITE(idev7,*) 'SQRT(EKS2) =',SQRT(EKS2/NC(NR)) 1111 CONTINUE C***************************************************************** C POSITIONS AND VELOCITIES CALCULATED C CHECK THAT THE SYSTEM IS AT REST C***************************************************************** 101 CONTINUE c***************************************************************** c create higher densities for 2D-case c if cink < 0 -> put vertical positions to zero c impacts take care of overlapps (hopefully) idime=3 if(cink.le.0) idime=2 if(cink.lt.0) then do 2222 i=1,lkm r(i,3)=0. rp(i,3)=0. 2222 continue endif c************************************************ c ADD INITIAL PERTURBATION 2/10/99 write(6,*) 'afidl=',afidl do 3010 i=1,i60 ampfidl(i)=0. phasefidl(i)=0. 3010 continue if(afidl.ne.0.) then c equal amplitude to all m-components if(afidl.gt.0) then c do 3000 i=1,i60 do 3000 i=m_afidl1,m_afidl2 ampfidl(i)=afidl phasefidl(i)=ran1(iidum)*2.*pii 3000 continue endif c read amplitudes from file if(afidl.lt.0) then open(unit=33 ,file='ampfidl.dat',status='old') do 3001 i=1,i60 ampfidl(i)=0. if(i.le.-afidl) then read(33,*) apu ampfidl(i)=apu phasefidl(i)=ran1(iidum)*2.*pii endif 3001 continue close(33) endif c impose perturbation do 3002 i=1,lkm farg=r(i,1)/xl*pii vperp=0. do 3003 m=1,i60 vperp=vperp+ampfidl(m)* + cos(m*(farg-phasefidl(m))) 3003 continue rp(i,1)=rp(i,1)+vperp*omega rp(i,2)=rp(i,2)+vperp*omega*afidly 3002 continue c afidl ne 0 endif write(71) i60 write(71) (ampfidl(ii),ii=1,i60) write(71) (phasefidl(ii),ii=1,i60) write(6,*) (ampfidl(ii),ii=1,i60) write(6,*) (phasefidl(ii),ii=1,i60) c*********************************************** c 17/01/03 assign spins do 1234 i=1,lkm rota(i,1)=0. rota(i,2)=0. rota(i,3)=qzinit*s(i)*w0 1234 continue if(rad_moon.gt.0) then rota(1,3)=rot_moon*omega*rad_moon write(6,*) 'wz/omega moon=',rot_moon write(6,*) 'wz/omega part=',rota(2,3)/omega/s(2) is(1)=0 num(1)=num(1)-1 endif c*********************************************** C ******************************************************** C IF IJATKO=1 READ PARTICLE POSITIONS FROM FILE 50 C ******************************************************** IF(IJATKO.EQ.1) THEN open(unit=50,file='inipos',status='old') read(50,*) time,lkm write(6,*) 'READ INITIAL POSITIONS FROM inipos' ISP=0 IF(LKM.LT.0) WRITE(6,*) 'INCLUDES ALSO SPINS' IF(LKM.LT.0) ISP=1 LKM=ABS(LKM) write(6,*) 'time,lkm',time,lkm read(50,*) (r(i,1),r(i,2),r(i,3),i=1,lkm) read(50,*) (rp(i,1),rp(i,2),rp(i,3),i=1,lkm) read(50,*) (myy(i),i=1,lkm) read(50,*) (s(i),i=1,lkm) read(50,*) (xmass(i),i=1,lkm) IF(ISP.EQ.1) read(50,*) + (ROTA(i,1),ROTA(i,2),ROTA(i,3),i=1,lkm) time0=time write(6,*) 'inipos-READ SUCCESFULL' close(50) if(n_moon.ne.0) then xmass(1)=0. myyspin=myy(1) write(6,*) 'after inipos: moonlet:' write(6,*) 'after inipos: xmass=',xmass(1) write(6,*) 'after inipos: myy=',myy(1) write(6,*) 'after inipos: myyspin=',myyspin endif ENDIF C *********************************************** UVAK=0.D0 WVAK=0.D0 ZVAK=0.D0 ZPVAK=0.D0 MSUM=0. EPOT=0. EKIN=0. DO 102 I=1,LKM UVAK=UVAK+RP(I,1)*MYY(I) WVAK=WVAK+(RP(I,2)+shear*R(I,1))*MYY(I) ZVAK=ZVAK+R(I,3)*MYY(I) ZPVAK=ZPVAK+RP(I,3)*MYY(I) MSUM=MSUM+MYY(I) EKIN=EKIN+0.5*MYY(I)*(RP(I,1)**2+RP(I,2)**2+RP(I,3)**2) EPOT=EPOT+0.5*MYY(I)*(factemp*R(I,1)**2+WZ2*R(I,3)**2) 102 CONTINUE UVAK=UVAK/W0/XL/MSUM WVAK=UVAK/W0/XL/MSUM ZVAK=ZVAK/MSUM ZPVAK=ZPVAK/W0/XL/MSUM uvak_all=uvak wvak_all=wvak zvak_all=zvak zpvak_all=zpvak ETOT=EKIN+EPOT ETOT0=ETOT WRITE(IDEV6,*) 'UVAK JA WVAK' WRITE(IDEV6,*) UVAK,WVAK WRITE(IDEV6,*) 'ZVAK JA ZPVAK' WRITE(IDEV6,*) ZVAK,ZPVAK WRITE(IDEV6,*) 'EKIN JA EPOT' WRITE(IDEV6,*) EKIN,EPOT WRITE(idev7,*) 'UVAK JA WVAK' WRITE(idev7,*) UVAK,WVAK WRITE(idev7,*) 'ZVAK JA ZPVAK' WRITE(idev7,*) ZVAK,ZPVAK WRITE(idev7,*) 'EKIN JA EPOT' WRITE(idev7,*) EKIN,EPOT tmass=0.d0 DO 13 I=1,LKM tmass=tmass+xmass(i) 13 CONTINUE write(*,*) 'rinit: tmass',tmass c store the current seed iidum_save=iidum RETURN END