-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathplayjet.pro_old
94 lines (79 loc) · 2.91 KB
/
playjet.pro_old
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
pro playjet,fname,var=var,sample=sample,startn=startn,endn=endn,step=step,png=png,dx=dx
device, decomposed=0
if not keyword_set(fname) then fname='JetSet_hdf5_plt_cnt_'
if not keyword_set(var) then var='dens'
;if not keyword_set(sample) then sample=0
if not keyword_set(step) then step=1
fnames = file_search(fname+'*')
nfiles = n_elements(fnames)
if keyword_set(startn) then begin
str_start = string(startn,format='(I4.4)')
start_j_tmp = where(strpos(fnames,str_start) ne -1)
start_j = start_j_tmp[0]
if start_j eq -1 then begin
print,'out of range of startnumber'
stop
endif
endif else start_j = 0
if keyword_set(endn) then begin
str_end = string(endn,format='(I4.4)')
end_j_tmp = where(strpos(fnames,str_end) ne -1)
end_j = end_j_tmp[0]
if end_j eq -1 then begin
print,'out of range of startnumber'
stop
endif
endif else end_j = nfiles-1
nfiles = end_j - start_j + 1
; Making PNG files
if keyword_set(png) then begin
spawn,'mkdir png_'+var
opt=',/pixmap'
endif else opt=''
xc0 = 0. & yc0 = 0. & zc0 = 0.
;maxd = 2.2d-13 & mind = 1.1d-15
;maxd = 2.76d-13 & mind = 2.9d-16
;maxd = 2.76d-13 & mind = 1.d-16
maxd = 2.76d-13 & mind = 1.d-17
m1 = 4.d34
m2 = 2.d34
mtot = m1+m2
peri = 3.d12
starrad = 1.4e12
th0 = !pi
for j=start_j,end_j,step do begin
print,j+1,' of',end_j+1,' ', fnames[j]
if (n_elements(dx) eq 0) then begin
if (n_elements(sample) eq 0) then begin
read_amr,fnames[j],var='dens',tree=tree,/nodata
smp = max(tree.lrefine) - 7
endif else smp = sample
d = dload(j,var=var,xc=xc0,yc=yc0,zc=zc0,x,y,z,sample=smp,time)
endif else begin
if (n_elements(sample) eq 0) then smp=0 else smp=sample
d = loaddata(fnames[j],var,xra=[xc0-dx,xc0+dx],yra=[yc0-dx,yc0+dx] $
,zra=[zc0-dx,zc0+dx],time=time,sample=smp,xcoords=x,ycoords=y,zcoords=z)
endelse
sd = size(d)
;maxd = max(d) & mind = min(d)
loadct,0,/sil
exestr = execute('window,0,xs='+strtrim(sd[1],1)+',ys='+strtrim(sd[2],1)+opt)
plot,x,y,/xst,/yst,/nodata,xr=[x[0],x[sd[1]-1]],yr=[y[0],y[sd[2]-1]],/iso,position=[0,0,sd[1],sd[2]],/dev
tv,bytscl(alog10(d[*,*,256]),max=alog10(maxd),min=alog10(mind))
; overorbit
th = th0 + sqrt(!unit.g * mtot / peri^3.d0) * time
; center of mass, position of BB
cm2xBh = peri*m1/mtot*cos(th)
cm2yBh = peri*m1/mtot*sin(th)
; center of mass, position of Star
cm2xSt = peri*m2/mtot*cos(th+!pi)
cm2ySt = peri*m2/mtot*sin(th+!pi)
plots,ring(cm2xSt,cm2ySt,starrad),/data,color=fsc_color('yellow')
plots,cm2xBh,cm2yBh,psym=1,color=fsc_color('yellow')
; legend,strtrim(time,1)+ ' s',/right,/top,text=fsc_color('yellow'),box=0
legend,string(sqrt(!unit.g * mtot / peri^3.d0)/2./!pi*time,format='(f4.2)')+' Period',/right,/top,box=0,textcolor=fsc_color('yellow')
if keyword_set(png) then $
draw,'png_'+var+'/Jet_plt_'+string(j,format='(I4.4)')
endfor
stop
end