35 | | |
36 | | |
37 | | |
38 | | |
39 | | |
40 | | {{{ |
41 | | #run it at command line with: |
42 | | #visit -assume_format Chombo -cli -s /Users/mhuartee/Box\ Sync/py/polarisation_xmdv_3.py |
43 | | #meant to produce DATE.okc file to be opened with visit then. |
44 | | |
45 | | MODE25D=1 |
46 | | MODE3D=0 |
47 | | |
48 | | mode=MODE25D |
49 | | |
50 | | SetActiveWindow(1) |
51 | | DeleteActivePlots() |
52 | | |
53 | | import math, datetime, sys, time, os |
54 | | #import numpy as np |
55 | | SuppressMessages(2) |
56 | | start = time.time()#measure program time |
57 | | |
58 | | |
59 | | # *** params *** |
60 | | angjet=0.0 |
61 | | tt=math.tan(angjet) |
62 | | cc=math.cos(angjet) |
63 | | ss=math.sin(angjet) |
64 | | DefineScalarExpression("angjet", str(angjet)) |
65 | | xcells=16 #from data |
66 | | pix=(8,8) #36min nx,ny |
67 | | ylim=(0.,8.) |
68 | | zlim=(28,36) #ignorde in MODE25D |
69 | | xlim=(0,8) |
70 | | dy=(ylim[1]-ylim[0])/float(pix[1]) |
71 | | dx=(xlim[1]-xlim[0])/float(pix[0]) |
72 | | dz=0 |
73 | | ############################### |
74 | | |
75 | | hoy=datetime.datetime.now().strftime("%m-%d") |
76 | | |
77 | | # z=jet |
78 | | # ^ |
79 | | # | |
80 | | # |---> y x axis is pointing out |
81 | | DefineScalarExpression("3dB_strength", "sqrt(Bx*Bx+By*By+Bz*Bz)") |
82 | | DefineScalarExpression("fact", "(px*px+py*py+pz*pz)/rho/2.")#Ekinetic |
83 | | |
84 | | if mode == MODE3D: #Create expressions in advance |
85 | | DefineScalarExpression("Bproj", "By*cos(angjet)+Bz*sin(angjet)") |
86 | | DefineScalarExpression("I", "fact*(Bx^2+Bproj^2)") |
87 | | DefineScalarExpression("Q", "fact*(Bx^2-Bproj^2)") |
88 | | DefineScalarExpression("U", "2.*fact*Bproj*Bx") |
89 | | elif mode == MODE25D: |
90 | | DefineScalarExpression("xp", "1.0") |
91 | | DefineScalarExpression("Bzz", "Bz*xp/x") |
92 | | DefineScalarExpression("Bxx", "Bz*sqrt(x^2-xp^2)/x") |
93 | | DefineScalarExpression("Bproj", "By*cos(angjet)+Bzz*sin(angjet)") |
94 | | DefineScalarExpression("dzdx", "x/sqrt(x^2-xp^2)") |
95 | | DefineScalarExpression("I","2*fact*( Bxx^2+Bproj^2) * dzdx") |
96 | | DefineScalarExpression("Q","2*fact*( Bxx^2-Bproj^2)*dzdx") |
97 | | DefineScalarExpression("U","2*2.*fact*Bxx*Bproj") |
98 | | |
99 | | #get data extents |
100 | | AddPlot("Pseudocolor", "rho", 1, 0) |
101 | | DrawPlots() |
102 | | SetQueryFloatFormat("%g") |
103 | | Query("SpatialExtents") |
104 | | xmax=float(GetQueryOutputValue()[1]) |
105 | | #ymax=GetQueryOutputValue()[3] |
106 | | #zmax=GetQueryOutputValue()[5] |
107 | | print xmax |
108 | | |
109 | | intcells=int( |
110 | | float(xcells)*(ylim[1]-ylim[0])/float(xmax) |
111 | | ) |
112 | | DeleteActivePlots() |
113 | | |
114 | | #lineout initialize |
115 | | AddPlot("Curve", "operators/Lineout/I", 1, 0) |
116 | | DrawPlots() |
117 | | LineoutAtts = LineoutAttributes() |
118 | | LineoutAtts.point1 = (0., 0., 32.) |
119 | | LineoutAtts.point2 = (xmax, 0., 32.) |
120 | | LineoutAtts.interactive = 0 |
121 | | LineoutAtts.ignoreGlobal = 0 |
122 | | LineoutAtts.samplingOn = 0 |
123 | | LineoutAtts.numberOfSamplePoints = 10 |
124 | | LineoutAtts.reflineLabels = 0 |
125 | | LineoutAtts.interactive = 1 |
126 | | LineoutAtts.numberOfSamplePoints = intcells |
127 | | |
128 | | |
129 | | #initial IO |
130 | | lines=pix[0]*pix[1] |
131 | | filename='/scratch/jcarrol5/'+hoy+'pol.okc' |
132 | | f = open(filename, 'w+') |
133 | | print f.write("%g %g %g\n" % (6,lines,6)) |
134 | | print f.write("x\ny\nz\npolx\npoly\npolz\n") |
135 | | print f.write("%g %g %g \n" % (ylim[0],float(1+pix[0])*dy, 10)) |
136 | | print f.write("%g %g %g \n" % (zlim[0],float(1+pix[1])*dz, 10)) |
137 | | print f.write("%g %g %g \n" % (0.,0.00, 10)) |
138 | | print f.write("%g %g %g \n" % (0.,0.75, 10)) |
139 | | print f.write("%g %g %g \n" % (0.,0.75, 10)) |
140 | | print f.write("%g %g %g \n" % (0.,0.00, 10)) |
141 | | |
142 | | #compute (integral in z) |
143 | | yshift=.5*tt*(zlim[1]-zlim[0]) |
144 | | z=0; |
145 | | x=xlim[0] |
146 | | while x<=xlim[1]: |
147 | | print str(int(98.*(x-xlim[0])/(xlim[1]-xlim[0])))+'%' |
148 | | |
149 | | if mode == MODE25D: #update expressions for integral calculations |
150 | | DefineScalarExpression("xp",str(x)) |
151 | | |
152 | | y=ylim[0] |
153 | | while y<=ylim[1]: |
154 | | |
155 | | if mode == MODE25D: #update expressions for integral calculations |
156 | | LineoutAtts.point1 = (x,y,0) |
157 | | LineoutAtts.point2 = (xlim[1],y+dz,0) |
158 | | elif mode == MODE3D: |
159 | | LineoutAtts.point1 = (x,y-yshift,zlim[0]) |
160 | | LineoutAtts.point2 = (x,y+yshift,zlim[1]) |
161 | | |
162 | | SetOperatorOptions(LineoutAtts, 0) |
163 | | |
164 | | |
165 | | #I |
166 | | ChangeActivePlotsVar("I") |
167 | | Query("Integrate") |
168 | | I0 = GetQueryOutputValue()/float(intcells) |
169 | | #Q |
170 | | ChangeActivePlotsVar("Q") |
171 | | Query("Integrate") |
172 | | Q0 = 0.75*GetQueryOutputValue()/float(intcells) |
173 | | #U |
174 | | ChangeActivePlotsVar("U") |
175 | | Query("Integrate") |
176 | | U0 = 0.75*GetQueryOutputValue()/float(intcells) |
177 | | |
178 | | if I0 != 0.: |
179 | | degree=( math.sqrt(math.pow(U0,2.) + |
180 | | math.pow(Q0,2.))/I0 |
181 | | ) |
182 | | else: |
183 | | degree=0 |
184 | | |
185 | | if Q0 != 0.: |
186 | | psi=0.5*(math.atan(U0/Q0)+0.*math.pi) |
187 | | else: |
188 | | psi=0 |
189 | | |
190 | | f.write("%g %g %g %g %g %g\n" % (x,y,z,degree*math.sin(psi), |
191 | | degree*math.cos(psi),0.)) |
192 | | |
193 | | |
194 | | y=y+dy |
195 | | x=x+dx |
196 | | |
197 | | f.close() |
198 | | stop = time.time() |
199 | | print '\n%3.1f min \n' % float((stop-start)/60.) |
200 | | print 'Done!, the xmdv file is called '+filename+'\n' |
201 | | #os.system('say "Your calculation is done"') |
202 | | }}} |