C ______________________________ Profiles for Galaxy Simulations C C INCLUDE 'PMparameters.h' INCLUDE 'PMgalaxy.h' Include 'leetipsy.h' REAL INPUT,moment Character FileASCII*100,Analysis*3 Character File1*100,File2*100 EQUIVALENCE (Ndisk,extras(90)), & (Nbulge,extras(97)),(aMassDisk,extras(93)), & (Box,extras(100)) Analysis ='PKDG' ScaleT =42.63!9.8e8 Box = 0.7 C................................................................... C Read data and open files write(*,*)' enter 5 digits of the snapshot' read(*,*) imoment c write(File1,'("FILES/dcs3.100.bs.pkd.",i5.5)')imoment c write(File1,'("FILES/dcs3.100.ssc.pkd.",i5.5)')imoment write(File1,'("FILES/KCB.pkd1.",i5.5)')imoment c write(*,*) File1 c File1 = 'FILES/Dcs3.std' write(*,*) File1 CALL leetipsy(File1) ScaleM =amscale*aminmax !1.e6/2. ! /2. write (*,*) ' out of tipsy scaling' write (*,*) ' time= ',time write (*,*) ' AEXPN=',AEXPN c imoment = 0 ISTEP =imoment CALL BLOCK_DATA ! Initialize arrays write(FileASCII,'("RESULTS/Result",a,".DM.",f6.4,".dat")') & Analysis,AEXPN OPEN(37,FILE=FileASCII,STATUS='UNKNOWN') write(FileASCII,'("RESULTS/Result",a,".SPH.",f6.4,".dat")') & Analysis,AEXPN OPEN(39,FILE=FileASCII,STATUS='UNKNOWN') write(FileASCII,'("RESULTS/Result",a,".CYL.",f6.4,".dat")') & Analysis,AEXPN OPEN(41,FILE=FileASCII,STATUS='UNKNOWN') write(FileASCII,'("RESULTS/Result",a,".BAR.",f6.4,".dat")') & Analysis,AEXPN OPEN(43,FILE=FileASCII,STATUS='UNKNOWN') OPEN(18,FILE='Stat.dat',STATUS='UNKNOWN',access='append') write(FileASCII,'("RESULTS/Result",a,".PRT.",f6.4,".dat")') & Analysis,AEXPN OPEN(20,FILE=FileASCII,STATUS='UNKNOWN') c OPEN(30,FILE='particlesu.dat',STATUS='UNKNOWN',form='UNFORMATTED') N = lspecies(Nspecies) aMassDisk = Ndisk*ScaleM time = AEXPN*1.e9 Nbulge = 0 write(*,*) ' aminmax= ',aminmax write(*,*)' Nt = ',N c pause write (*,*) ' mass of disk=',aMassDisk write (*,*) ' RDTAPE is done Time=',time,' a=',AEXPN WRITE (*,100) time/1.e9, + Ndisk,Nbulge,aMassDisk, + bar_max*1.e3,dZ_slice*1.e3 Do i=37,43,2 WRITE (i,100) time/1.e9, + Ndisk,Nbulge,aMassDisk, + bar_max*1.e3,dZ_slice*1.e3 EndDo 100 FORMAT(T2,'Time(Gyrs)', & T15,'Ndisk',T25,'Nbulge',T35,'Mdisk(Msun)',/ & T2,f8.5,T15,i8,T25,i7,T35,g11.4,/ + ' Max Bar radius(kpc)=', & f7.3,' Slice in Z for cylindr shells=',f8.3) c Box =INPUT(' Enter box size in comoving Mpc/h =') c Box =0.5 ! NG3109 c Box = 2 ! 'bulge' BoxV = Box*100. ! Box size in km/s ScaleV = BoxV/AEXPN/NGRID ! scale factor for Velocities (km/s) ScaleV = 1.0 ScaleH = 0.! Hubble flow correction if(Ndisk.le.0)write (*,*) ' wrong number of disk particles:',Ndisk if(Ndisk.le.0)Stop Call ChangePnts(N,Box,ScaleH) c Call DumpPnts(N,Box,ScaleM) c Call AngMom(N,30000,3.,Box,ScaleM) c CALL PV(N,Box,ScaleM) c CALL RandV(N,Box,ScaleM) Call Profile(N,Box,ScaleM) Call Statist(N,Box,ScaleM) c Call Projection(N,Box,0.) c Call Projection(N,Box,5.) END c-------------------------------------------------------------------- Subroutine leetipsy(filename_in) c-------------------------------------------------------------------- include 'PMparameters.h' include 'leetipsy.h' COMMON /GENERAL/time double precision time integer ndark c double precision time real xx(Nmaxpart),yy(Nmaxpart),zz(Nmaxpart),mass(Nmaxpart) real vxx(Nmaxpart),vyy(Nmaxpart),vzz(Nmaxpart) real eps(Nmaxpart), phi(Nmaxpart) real soft(Nmaxpart),metals(Nmaxpart),tform(Nmaxpart),rho(Nmaxpart) real temp(Nmaxpart),pot(Nmaxpart) character*100 filename_in,filename_out EQUIVALENCE (Ndisk,extras(90)), & (Nbulge,extras(97)),(aMassDisk,extras(93)), & (Box,extras(100)) c filename_in= 'MW1.1024.00512' c filename_in= 'dwf1.6144.02048' write (*,*) ' in lltipsy' Call READPARAM() write (*,*) ' Open files' call open_file(filename_in,icheck) c call readdmtip(ndark,nstar,nsph,xx(1),yy(1),zz(1),vx(1),vy(1), c & vz(1),mass(1),eps(1),phi(1),time,soft(1),metals(1),tform(1), c & rho(1),temp(1),pot(1)) c OR, FOR STANDARD BINARY: write (*,*) ' read files' call readdmstd(ndark,nstar,nsph,xx(1),yy(1),zz(1),vxx(1),vyy(1), & vzz(1),mass(1),eps(1),phi(1),time,soft(1),metals(1),tform(1), & rho(1),temp(1),pot(1)) nt = ndark + nstar + nsph ndm = ndark nd = nstar aminmax = mass(ndark+1) write(*,*)' nt= ',nt,ndark,nstar write(*,*)' time =',time*42.65 write(*,*)' aminmax= ',aminmax,aminmax*amscale tiempo = time AEXPN = time call Translate(xx,yy,zz,vxx,vyy,vzz,mass,nt,nd,ndm,tiempo) call close_file return end c-------------------------------------------------------------------- Subroutine Translate(xxn,yyn,zzn,vxxn,vyyn,vzzn,mass,nt,nd, * ndm,tiempo) c-------------------------------------------------------------------- Include 'PMparameters.h' Include 'PMgalaxy.h' Include 'leetipsy.h' INTEGER INDX(10),ispecies(10),peso(10),ispecies2(10) integer ndark,nd,ndm c double precision time real xxn(Nmaxpart),yyn(Nmaxpart),zzn(Nmaxpart),mass(Nmaxpart) real vxxn(Nmaxpart),vyyn(Nmaxpart),vzzn(Nmaxpart) real eps(Nmaxpart), phi(Nmaxpart) real soft(Nmaxpart),metals(Nmaxpart),tform(Nmaxpart), & rho(Nmaxpart) real temp(Nmaxpart),pot(Nmaxpart) EQUIVALENCE (Ndisk,extras(90)), & (Nbulge,extras(97)),(aMassDisk,extras(93)), & (Box,extras(100)) hsmall = hubble aexpn = tiempo*42.63 ! time in gyrs Ndisk = nd write(*,*)'time= ',aexpn write(*,*)'Box= ',Box xmin =1.e10 xmax =-1.e10 ymin =1.e10 ymax =-1.e10 zmin =1.e10 zmax =-1.e10 Do i=1,10 peso(i) = 0.0 ispecies(i) = 0 ispecies2(i) = 0 enddo Nspecies = 1 write(*,*)'Nt= ',nt,' Nd= ',nd,' Ndark= ',ndm istar = 0 Do i=1,nt itemp = int(mass(i)/aminmax + 0.5) iWeight(i) = float(itemp) if(itemp.lt.1)then write(*,*)'test lspecies',itemp,mass(i),aminmax,i stop endif j = 1 11 if( peso(j).eq.iWeight(i))then !if not ispecies(j) = ispecies(j) + 1 else if(peso(j).ne.iWeight(i).and.peso(j).ne.0.)then j = j + 1 goto 11 else if(peso(j).ne.iWeight(i).and.peso(j).eq.0.)then peso(j) = iWeight(i) ispecies(j) = ispecies(j) + 1 write(*,*)'species test',j,iWeight(i) endif c Go to proper kpc, km/s xxn(i) = xxn(i)*rscale/1.e3 + Box/2! ! 1.e3 from Kpc to Mpc yyn(i) = yyn(i)*rscale/1.e3 + Box/2 zzn(i) = zzn(i)*rscale/1.e3 + Box/2 vxxn(i) = vxxn(i)*vscale vyyn(i) = vyyn(i)*vscale vzzn(i) = vzzn(i)*vscale iWeight(i) = peso(j) if(i.gt.ndm)then istar = istar +1 xp(istar) = xxn(i) yp(istar) = yyn(i) zp(istar) = zzn(i) vxp(istar) = vxxn(i) vyp(istar) = vyyn(i) vzp(istar) = vzzn(i) iWeight(istar) = iWeight(i) endif enddo write(*,*)'loop over particles in traslate done' c Reorder Itest = 0 Do k =1,10 Itest = Itest + ispecies(k) if(ispecies(k).gt.0.)Nspecies = k enddo Do k=1,Nspecies,1 ij =Nspecies +1 - k lspecies(k) = ispecies(ij) wspecies(k) = peso(ij) write(*,*)'specie # ',k,lspecies(k),wspecies(k) if(k.ge.2)then lspecies(k) = lspecies(k) + lspecies(k-1) endif Enddo Do k=1,Nspecies,1 write(*,*)'lspecies',k,lspecies(k),wspecies(k) Enddo c Particle Reshuffle Do kk=1,nt-nd ifactor = int(mass(kk)/aminmax + 0.5) jj=1 22 if(ifactor.eq.wspecies(jj))then ispecies2(jj) = ispecies2(jj) +1 if(jj.eq.1)then !1st specie indext = ispecies2(jj) + nd xp(indext) = xxn(kk) yp(indext) = yyn(kk) zp(indext) = zzn(kk) vxp(indext) = vxxn(kk) vyp(indext) = vyyn(kk) vzp(indext) = vzzn(kk) iWeight(indext) = float(ifactor) xmin =min(Xp(indext),xmin) xmax =max(Xp(indext),xmax) ymin =min(Yp(indext),ymin) ymax =max(Yp(indext),ymax) zmin =min(Zp(indext),zmin) zmax =max(Zp(indext),zmax) else if(jj.ne.1)then indext = lspecies(jj-1) if(indext.le.0.0)then write(*,*)'stop= ',indext,jj,kk stop endif indext = indext + ispecies2(jj) xp(indext) = xxn(kk) yp(indext) = yyn(kk) zp(indext) = zzn(kk) vxp(indext) = vxxn(kk) vyp(indext) = vyyn(kk) vzp(indext) = vzzn(kk) iWeight(indext) = ifactor xmin =min(Xp(indext),xmin) xmax =max(Xp(indext),xmax) ymin =min(Yp(indext),ymin) ymax =max(Yp(indext),ymax) zmin =min(Zp(indext),zmin) zmax =max(Zp(indext),zmax) endif else jj = jj + 1 goto 22 endif Enddo imark1 = lspecies(1) imark2 = lspecies(2) imarkn= lspecies(nspecies) c write(*,*)'Ngrid, Nrow, w(1)',Ngrid,Nrow,wspecies(1) c write(*,*)'test lspecies1,2, n', lspecies(1),lspecies(2),imarkn c write(*,*)'test weights 1 N',iWeight(1), c * iWeight(imark2), iWeight(imarkn) c write(*,*)'test reordering',iWeight(imark2)/max(iWeight(1),1.), c * iWeight(imarkn)/max(iWeight(imark1),1.) write(*,*)'Nspecies= ',Nspecies write(*,*)'Nparticles= ',Itest Do i=1,Nspecies c wspecies(i) = wspecies(i)*weightSmall c Dscale= 2.746e+11*hsmall**2*(Box/NGRID)**3 *wspecies(1) write(*,*)'test w-mass',i, wspecies(i) enddo write(*,*) ' x:',xmin,xmax write(*,*) ' y:',ymin,ymax write(*,*) ' z:',zmin,zmax return end c-------------------------------------------------------------------- Subroutine READPARAM() c-------------------------------------------------------------------- include 'PMparameters.h' include 'leetipsy.h' HEADER = HEADER2 AEXPN = AEXPN2 AEXPN0 = AEXPN02 ASTEP = ASTEP2 ISTEP = ISTEP2 Om0 = Om02 Oml0 = Oml02 hubble = hubble2 Ocurv = Ocurv2 PARTW = PARTW2 AMPLT = AMPLT2 EKIN = EKINB EKIN1 = EKIN1B EKIN2 = EKIN2B BOX = BOX2 c extras(100) = BOX Return End