PROGRAM setspiral_gadget C C this program reads GADGET-2 initial conditions C and adds the initial conditions. C C Author: Thorsten Naab C Modified for GADGET-2: PHJ 12/12/2006 C C C C Cc IMPLICIT NONE integer anstot,bnstot,nstot parameter(anstot=4400000,bnstot=4400000,nstot=9000000) integer anbodies,bnbodies,nd,nb,folge,parorb, $ byhand,p,j,k,i,npart,jg REAL amass,apos,avel,au,bmass,bpos,bvel,bu,tnow,nbodies, & mass,pos,vel,u,thetmod1,thetmod2,phimod1,phimod2, & costhet1,costhet2,sinthet1,sinthet2,cosphi1,cosphi2, & sinphi1,sinphi2,one,pi,arot,brot,avrot,bvrot,m1,m2, $ mtotal,rmod1,rmod2,vsep,xmod1,ymod1,zmod1,xmod2, $ ymod2,zmod2,rsep,rp,rseptest,vxmod1,vymod1,vzmod1, $ vxmod2,vymod2,vzmod2,vsep1,vsep2,eta,eta1,mr,vcm, $ pcm,mtot REAL GRAV integer massloop C logical folge dimension amass(anstot),apos(anstot,3),avel(anstot,3),au(anstot), & bmass(bnstot),bpos(bnstot,3),bvel(bnstot,3),bu(bnstot), & mass(nstot),pos(nstot,3),arot(anstot,3),brot(bnstot,3) & ,vel(nstot,3),u(nstot),avrot(anstot,3),bvrot(bnstot,3), & vcm(3),pcm(3) integer ndim,loop parameter(ndim=3) integer*4 npart_a(0:5), nall_a(0:5) real*8 massarr_a(0:5), time_a, redshift_a integer*4 unused_a(24) integer*4 flag_sfr_a,flag_feedback_a,flag_cool_a,numfiles_a real*8 boxsize_a,omega0_a,omegav_a,h0_a integer*4 aid(anstot) integer anwithmass real amassdum(anstot) character*40 afile,bfile integer*4 npart_b(0:5), nall_b(0:5) real*8 massarr_b(0:5), time_b, redshift_b integer*4 unused_b(24) integer*4 flag_sfr_b,flag_feedback_b,flag_cool_b,numfiles_b real*8 boxsize_b,omega0_b,omegav_b,h0_b integer*4 bid(bnstot) integer bnwithmass real bmassdum(anstot) integer*4 npart_comb(0:5), nall_comb(0:5),ntot_comb real*8 massarr_comb(0:5), time_comb, redshift_comb integer*4 unused_comb(24) integer*4 flag_sfr_comb,flag_feedback_comb, & flag_cool_comb,numfiles_comb real*8 boxsize_comb,omega0_comb,omegav_comb,h0_comb real*4 pos_comb(3,nstot),vel_comb(3,nstot), & mass_comb(nstot),u_comb(nstot) integer*4 id_comb(nstot) character*50 combfile ! ! Galaxy a ! write(6,15) 15 FORMAT( ' Input file for disk 1? ',$) read(5,*) afile print*, 'file 1: ', afile open (1, file=afile, ! $ 'ics_SW_Y_mhres.dat', $ form='unformatted') read (1) npart_a, massarr_a,time_a, redshift_a, flag_sfr_a, $ flag_feedback_a, nall_a, flag_cool_a, numfiles_a, $ boxsize_a,omega0_a,omegav_a,h0_a,unused_a write(6,*) npart_a, massarr_a anwithmass=0 anbodies=0 do i=0,5 anbodies=anbodies+npart_a(i) if ((npart_a(i).gt.0).AND.(massarr_a(i).eq.0.)) then anwithmass=anwithmass+npart_a(i) endif enddo write(6,*) anwithmass, anbodies, anstot-anbodies read(1) ((apos(i,k),k=1,ndim),i=1,anbodies) read(1) (((avel(i,k)),k=1,ndim),i=1,anbodies) read(1) (aid(i),i=1,anbodies) if (anwithmass.gt.0) then read(1) (amassdum(i),i=1,anwithmass) endif read(1) (au(i),i=1,npart_a(0)) loop=0 massloop=0 do i=0,5 if (massarr_a(i) .ne. 0.) then do k=0,npart_a(i)-1 loop=loop+1 amass(loop)=massarr_a(i) enddo else do k=0,npart_a(i)-1 loop=loop+1 massloop=massloop+1 amass(loop)=amassdum(massloop) enddo endif enddo close(1) ! ! Galaxy b ! write(6,25) 25 FORMAT( ' Input file for disk 2? ',$) read(5,*) bfile print*, 'file 2: ', bfile open (2,file=bfile, ! $ file='ics_SW_Y_mhres.dat', $ form='unformatted') read (2) npart_b, massarr_b ,time_b, redshift_b, flag_sfr_b, $ flag_feedback_b, nall_b, flag_cool_b, numfiles_b, $ boxsize_b,omega0_b,omegav_b,h0_b,unused_b write(6,*) npart_b, massarr_b bnwithmass=0 bnbodies=0 do i=0,5 bnbodies=bnbodies+npart_b(i) if ((npart_b(i).gt.0).AND.(massarr_b(i).eq.0.)) then bnwithmass=bnwithmass+npart_b(i) endif enddo write(6,*) bnwithmass, bnbodies, bnstot-bnbodies read(2) ((bpos(i,k),k=1,ndim),i=1,bnbodies) read(2) (((bvel(i,k)),k=1,ndim),i=1,bnbodies) read(2) (bid(i),i=1,bnbodies) if (bnwithmass.gt.0) then read(2) (bmassdum(i),i=1,bnwithmass) endif read(2) (bu(i),i=1,npart_b(0)) loop=0 massloop=0 do i=0,5 if (massarr_b(i) .ne. 0.) then do k=0,npart_b(i)-1 loop=loop+1 bmass(loop)=massarr_b(i) enddo else do k=0,npart_b(i)-1 loop=loop+1 massloop=massloop+1 bmass(loop)=bmassdum(massloop) enddo endif enddo close(2) ! ! OUTPUT file ! !!!!!!!!!!!!!!!!!!!!!! C set disk inclination C C C WRITE(6,70) 70 FORMAT( ' Rotation angles (degrees) theta, phi for disk 1? ',$) READ(5,*) thetmod1,phimod1 print* ,thetmod1,phimod1 WRITE(6,60) 60 FORMAT( ' Rotation angles (degrees) theta, phi for disk 2? ',$) READ(5,*) thetmod2,phimod2 print* ,thetmod2,phimod2 GRAV=43007.1D0 one=1.0d0 pi=4.0d0*ATAN(one) thetmod1=2.0d0*pi*thetmod1/360.0d0 phimod1=2.0d0*pi*phimod1/360.0d0 thetmod2=2.0d0*pi*thetmod2/360.0d0 phimod2=2.0d0*pi*phimod2/360.0d0 costhet1=COS(thetmod1) sinthet1=SIN(thetmod1) cosphi1=COS(phimod1) sinphi1=SIN(phimod1) costhet2=COS(thetmod2) sinthet2=SIN(thetmod2) cosphi2=COS(phimod2) sinphi2=SIN(phimod2) DO i=1,anstot arot(i,1) = apos(i,1)*cosphi1 & +apos(i,2)*sinphi1*costhet1+ & apos(i,3)*sinphi1*sinthet1 arot(i,2)= -apos(i,1)*sinphi1+ & apos(i,2)*cosphi1*costhet1+ & apos(i,3)*cosphi1*sinthet1 arot(i,3)= -apos(i,2)*sinthet1+apos(i,3)*costhet1 ENDDO DO i=1,bnstot brot(i,1) = bpos(i,1)*cosphi2 & +bpos(i,2)*sinphi2*costhet2+ & bpos(i,3)*sinphi2*sinthet2 brot(i,2)= -bpos(i,1)*sinphi2+ & bpos(i,2)*cosphi2*costhet2+ & bpos(i,3)*cosphi2*sinthet2 brot(i,3)= -bpos(i,2)*sinthet2+bpos(i,3)*costhet2 ENDDO DO i=1,anstot avrot(i,1) = avel(i,1)*cosphi1 & +avel(i,2)*sinphi1*costhet1+ & avel(i,3)*sinphi1*sinthet1 avrot(i,2)= -avel(i,1)*sinphi1+ & avel(i,2)*cosphi1*costhet1+ & avel(i,3)*cosphi1*sinthet1 avrot(i,3)= -avel(i,2)*sinthet1+avel(i,3)*costhet1 ENDDO DO i=1,bnstot bvrot(i,1) = bvel(i,1)*cosphi2 & +bvel(i,2)*sinphi2*costhet2+ & bvel(i,3)*sinphi2*sinthet2 bvrot(i,2)= -bvel(i,1)*sinphi2+ & bvel(i,2)*cosphi2*costhet2+ & bvel(i,3)*cosphi2*sinthet2 bvrot(i,3)= -bvel(i,2)*sinthet2+bvel(i,3)*costhet2 ENDDO C C set distance C C print*, $ 'Want to set initial seperation and velocities by hand (0/1)?' c read*,byhand byhand=0 IF (byhand .EQ. 1) THEN DO 10 i=1,anbodies mass(i)=amass(i) pos(i,1)=arot(i,1) pos(i,2)=arot(i,2) pos(i,3)=arot(i,3) vel(i,1)=avrot(i,1) vel(i,2)=avrot(i,2) vel(i,3)=avrot(i,3) 10 continue DO 20 i=1,bnbodies mass(i+anbodies)=bmass(i) pos(i+anbodies,1)=brot(i,1)+10.d0 pos(i+anbodies,2)=brot(i,2) pos(i+anbodies,3)=brot(i,3)-50.d0 vel(i+anbodies,1)=bvrot(i,1) vel(i+anbodies,2)=bvrot(i,2) vel(i+anbodies,3)=bvrot(i,3)+2.d0 20 continue ENDIF print*,'Want to choose parabolic orbit (0/1) ?' c read*,parorb parorb=1 IF ((parorb .EQ.1) .and. (byhand .EQ. 1)) THEN print*,'Incompatible input!' stop ENDIF c Units are kpc/h IF (parorb .EQ. 1) THEN print*,' Pericenter distance for parabolic orbit? ' read*,rp ! rp= 2.5d0* 0.5d0 * (3.12277d0+2.05704d0)! 2.d0* 0.5d0* (2.92606d0+4.4571d0) print*, 'rp= ', rp print*,'Initial center of mass separation? ' read*,rsep ! rsep= 127.d0!119. print*, 'rsep= ', rsep C********************************************************************** C C Compute total mass C********************************************************************** m1=0.0 m2=0.0 DO i=1,anbodies m1=m1+amass(i) ENDDO print*,'mass of disk 1',m1 DO i=1,bnbodies m2=m2+bmass(i) ENDDO print*,'mass of disk 2',m2 mtotal = m1+m2 print*,'mtotal = ',mtotal C********************************************************************** C Center of mass C********************************************************************** rmod1=rsep - m1*rsep/(m1+m2) rmod1=abs(rmod1) print*,'rmod1 = ',rmod1 rmod2=rsep -rmod1 rmod2=abs(rmod2) print*,'rmod2 = ',rmod2 xmod1=m2/mtotal*(rsep-2.0*rp) ymod1=2.0*m2/mtotal*sqrt(rsep*rp-rp*rp) zmod1=0.0 print*,'xmod1',xmod1 print*,'ymod1',ymod1 xmod2=-m1/mtotal*(rsep-2.0*rp) ymod2=-2.0*m1/mtotal*sqrt(rsep*rp-rp*rp) zmod2=0.0 print*,'rsept = ',rmod1+rmod2 print*,'xmod2',xmod2 print*,'ymod2',ymod2 C********************************************************************** C reduced mass C********************************************************************** mr=m1*m2/(m1+m2) vxmod1=-m2/rsep*sqrt((2.0*rsep-2.0*rp)/mtotal)*sqrt(GRAV) print*,'vxmod1 = ',vxmod1 vymod1= -m2/rsep*sqrt(2*rp/mtotal)*sqrt(GRAV) print*,'vymod1 = ',vymod1 vzmod1=0.0 vxmod2=m1/rsep*sqrt((2.0*rsep-2.0*rp)/mtotal)*sqrt(GRAV) print*,'vxmod2 = ',vxmod2 vymod2=m1/rsep*sqrt(2*rp/mtotal)*sqrt(GRAV) print*,'vymod2 = ',vymod2 vzmod2=0.0 vsep=SQRT((vxmod1-vxmod2)**2+(vymod1-vymod2)**2+ & (vzmod1-vzmod2)**2) print*,'vsep = ',vsep rseptest= SQRT((xmod1-xmod2)**2+(ymod1-ymod2)**2+ & (zmod1-zmod2)**2) print*,'rseptest',rseptest if(abs(rseptest-rsep).gt.1.e-3)then print*,'Something is wrong with distances' stop endif vsep1=vxmod1**2.+vymod1**2. vsep2=vxmod2**2.+vymod2**2. eta=0.5*m1*vsep1+0.5*m2*vsep2-GRAV*m1*m2/rseptest print*,'eta = ',eta/GRAV if (abs(eta/GRAV).gt. 1.e-4) then print*,'This is not a parabolic orbit!!!' stop endif DO i=1,anbodies mass(i)=amass(i) pos(i,1)=arot(i,1)+xmod1 pos(i,2)=arot(i,2)+ymod1 pos(i,3)=arot(i,3)+zmod1 vel(i,1)=avrot(i,1)+vxmod1 vel(i,2)=avrot(i,2)+vymod1 vel(i,3)=avrot(i,3)+vzmod1 ENDDO DO i=1,bnbodies mass(i+anbodies)=bmass(i) pos(i+anbodies,1)=brot(i,1)+xmod2 pos(i+anbodies,2)=brot(i,2)+ymod2 pos(i+anbodies,3)=brot(i,3)+zmod2 vel(i+anbodies,1)=bvrot(i,1)+vxmod2 vel(i+anbodies,2)=bvrot(i,2)+vymod2 vel(i+anbodies,3)=bvrot(i,3)+vzmod2 ENDDO ENDIF DO i=1,npart_a(0) u(i)=au(i) enddo DO i=1,npart_b(0) u(i+npart_a(0))=bu(i) enddo npart=int(anbodies+bnbodies) nbodies=anbodies+bnbodies C C output to file TREEBI C mtot=0. DO p=1,nbodies mtot=mtot+mass(p) ENDDO print*,'mtot',mtot DO k=1,ndim DO p=1,nbodies pcm(k)=pos(p,k)+mass(p)*pos(p,k) vcm(k)=vel(p,k)+mass(p)*vel(p,k) ENDDO pcm(k)=pcm(k)/mtot vcm(k)=vcm(k)/mtot ENDDO print*,'cmpos',pcm print*,'cmvel',vcm C Construct the combined array C ----------------------------------------------------------------------------- jg=1 ! gas 1 do i=1,npart_a(0) do k=1,3 pos_comb(k,jg)=pos(i,k) vel_comb(k,jg)=vel(i,k) enddo mass_comb(jg)=mass(i) id_comb(jg)=jg u_comb(jg)=u(jg) ! write(6,*) u(i),u(jg),jg ! write(6,*) jg,i jg=jg+1 enddo write(6,*) 'NGas Gal a',npart_a(0) ! gas 2 do i=1+anbodies,anbodies+npart_b(0) do k=1,3 pos_comb(k,jg)=pos(i,k) vel_comb(k,jg)=vel(i,k) enddo mass_comb(jg)=mass(i) id_comb(jg)=jg u_comb(jg)=u(jg) ! write(6,*) u(i),u(jg),jg ! write(6,*) jg,i jg=jg+1 enddo write(6,*) 'NGas Gal b',npart_b(0) ! halo 1 do i=1+npart_a(0),npart_a(0)+npart_a(1) do k=1,3 pos_comb(k,jg)=pos(i,k) vel_comb(k,jg)=vel(i,k) enddo id_comb(jg)=jg mass_comb(jg)=mass(i) ! write(6,*) jg,i jg=jg+1 enddo write(6,*) 'Nhalo Gal a',npart_a(1) ! halo 2 do i=1+anbodies+npart_b(0),anbodies+npart_b(0)+npart_b(1) do k=1,3 pos_comb(k,jg)=pos(i,k) vel_comb(k,jg)=vel(i,k) enddo id_comb(jg)=jg mass_comb(jg)=mass(i) ! write(6,*) jg,i jg=jg+1 enddo write(6,*) 'Nhalo Gal b',npart_b(1) ! disk 1 do i=1+npart_a(0)+npart_a(1),npart_a(0)+npart_a(1)+npart_a(2) do k=1,3 pos_comb(k,jg)=pos(i,k) vel_comb(k,jg)=vel(i,k) enddo id_comb(jg)=jg mass_comb(jg)=mass(i) ! write(6,*) jg,i jg=jg+1 enddo write(6,*) 'Ndisk Gal a',npart_a(2) ! disk 2 do i=1+anbodies+npart_b(0)+npart_b(1), & anbodies+npart_b(0)+npart_b(1)+npart_b(2) do k=1,3 pos_comb(k,jg)=pos(i,k) vel_comb(k,jg)=vel(i,k) enddo id_comb(jg)=jg mass_comb(jg)=mass(i) ! write(6,*) jg,i jg=jg+1 enddo write(6,*) 'Ndisk Gal b',npart_b(2) ! bulge 1 do i=1+npart_a(0)+npart_a(1)+npart_a(2), & npart_a(0)+npart_a(1)+npart_a(2)+npart_a(3) do k=1,3 pos_comb(k,jg)=pos(i,k) vel_comb(k,jg)=vel(i,k) enddo id_comb(jg)=jg mass_comb(jg)=mass(i) ! write(6,*) jg,i jg=jg+1 enddo write(6,*) 'Nbulge Gal a',npart_a(3) ! bulge 2 do i=1+anbodies+npart_b(0)+npart_b(1)+npart_b(2), & anbodies+npart_b(0)+npart_b(1)+npart_b(2)+npart_b(3) do k=1,3 pos_comb(k,jg)=pos(i,k) vel_comb(k,jg)=vel(i,k) enddo id_comb(jg)=jg mass_comb(jg)=mass(i) ! write(6,*) jg,i jg=jg+1 enddo write(6,*) 'Nbulge Gal b',npart_b(3) ! stars 1 do i=1+npart_a(0)+npart_a(1)+npart_a(2)+npart_a(3), & npart_a(0)+npart_a(1)+npart_a(2)+npart_a(3)+npart_a(4) do k=1,3 pos_comb(k,jg)=pos(i,k) vel_comb(k,jg)=vel(i,k) enddo id_comb(jg)=jg mass_comb(jg)=mass(i) ! write(6,*) jg,i jg=jg+1 enddo write(6,*) 'Nstars Gal a',npart_a(4) ! stars 2 do i=1+anbodies+npart_b(0)+npart_b(1)+npart_b(2)+npart_b(3), & anbodies+npart_b(0)+npart_b(1)+npart_b(2)+npart_b(3)+npart_b(4) do k=1,3 pos_comb(k,jg)=pos(i,k) vel_comb(k,jg)=vel(i,k) enddo id_comb(jg)=jg mass_comb(jg)=mass(i) ! write(6,*) jg,i jg=jg+1 enddo write(6,*) 'Nstars Gal b',npart_b(4) ! blackhole 1 do i=1+npart_a(0)+npart_a(1)+npart_a(2)+npart_a(3)+npart_a(4), & npart_a(0)+npart_a(1)+npart_a(2)+npart_a(3)+ & npart_a(4)+npart_a(5) do k=1,3 pos_comb(k,jg)=pos(i,k) vel_comb(k,jg)=vel(i,k) enddo id_comb(jg)=jg mass_comb(jg)=mass(i) ! write(6,*) jg,i jg=jg+1 enddo write(6,*) 'NBH Gal a',npart_a(5) ! blackhole 2 do i=1+anbodies+npart_b(0)+npart_b(1)+npart_b(2)+ & npart_b(3)+npart_b(4), & anbodies+npart_b(0)+npart_b(1)+npart_b(2)+ & npart_b(3)+npart_b(4)+npart_b(5) do k=1,3 pos_comb(k,jg)=pos(i,k) vel_comb(k,jg)=vel(i,k) enddo id_comb(jg)=jg mass_comb(jg)=mass(i) ! write(6,*) jg,i jg=jg+1 enddo write(6,*) 'NBH Gal b',npart_b(5) C Write the combined GADGET file with the mergers C ----------------------------------------------------------------------------- ntot_comb=0 do i=0,5 npart_comb(i)=npart_a(i)+npart_b(i) massarr_comb(i)=0 nall_comb(i)=nall_a(i)+nall_b(i) ntot_comb=ntot_comb+npart_comb(i) enddo write(6,*) 'Ntot Comb Ntot gas',ntot_comb,npart_comb(0) time_comb=time_a redshift_comb=redshift_a flag_sfr_comb=flag_sfr_a flag_feedback_comb=flag_feedback_a flag_cool_comb=flag_cool_a numfiles_comb=numfiles_a boxsize_comb=boxsize_a omega0_comb=omega0_a omegav_comb=omegav_a h0_comb=h0_a write(6,35) 35 FORMAT( ' Output filename? ',$) read(5,*) combfile print*, 'Output file: ', combfile open (3, $ file=combfile, ! $ 'ics_combined_SW_Y_mhres.dat', $ form='unformatted') write(3) npart_comb, massarr_comb ,time_comb, redshift_comb, $ flag_sfr_comb,flag_feedback_comb, nall_comb, flag_cool_comb, $ numfiles_comb,boxsize_comb,omega0_comb,omegav_comb,h0_comb, $ unused_comb write(3) ((pos_comb(k,i),k=1,ndim),i=1,ntot_comb) write(3) (((vel_comb(k,i)),k=1,ndim),i=1,ntot_comb) write(3) (id_comb(i),i=1,ntot_comb) write(3) (mass_comb(i),i=1,ntot_comb) write(3) (u_comb(i),i=1,npart_comb(0)) ! write(6,*) ((pos_comb(k,i),k=1,ndim),i=1,ntot_comb) ! write(6,*) (((vel_comb(k,i)),k=1,ndim),i=1,ntot_comb) ! write(6,*) (id_comb(i),i=1,ntot_comb) ! write(6,*) (mass_comb(i),i=1,100) ! write(6,*) (u_comb(i),i=1,npart_comb(0)) ! write(6,*) npart_comb(0) close(3) print*,'nbodies',nbodies print*,'ndim',ndim print*,'amass(1)',amass(1) print*,'apos(1,1)',apos(1,1) print*,'bmass(1)',bmass(1) print*,'bpos(1,1)',bpos(1,1) print*,'avel(1,1)',avel(1,1) print*,'bvel(1,1)',bvel(1,1) C print*,'bvel(60000,1)',bvel(60000,1) C print*,'bvel(60000,1)',bvel(60000,2) print*,'vel(1,1)',vel(1,1) print*,'vel(1,anbodies+1)',vel(anbodies+1,1) print*,'vel(1,nbodies+1)',vel(nbodies,1) print*,'vel_comb(1,1)',vel_comb(1,1) print*,'vel_comb(1,npart_a(0)+1)',vel_comb(1,npart_a(0)+1) print*,'vel_comb(1,nbodies+1)',vel_comb(1,nbodies) end