Goal is to integrate various functions along lines of sight which is constrained to be in the y-z plane.

\int f(B_\perp(s), B_x(s)) ds

B_\perp(s) = B_y(s) \cos(\theta) + B_z(s) \sin(\theta)

Each pixel will be given by the integral of the function along the path

s = x_i, y_i+s \sin(\theta), z_i+s \cos(\theta)

where \delta y = \cos(\theta) \delta x and \delta z =\sin(\theta) \delta x

For 3D we just need to create the expressions in visit, create the expressions and the lineouts, and then integrate the resulting query

For 2.5D, we need to transform the integral along s into an integral in the xy plane.

Assume jet is oriented along y and that we are integrating along the direction (0, sin(\theta), cos(\theta))

where \theta is the angle of inclination.

We are interested in calculating the integral of B_p^2, B_x^2, and B_xB_p where B_p = cos(\theta)By + sin(\theta)Bz

If \theta = 0 this is just \int B_y^2 dz, \int B_x^2 dz, and \int B_x B_y dz

If \theta /= 0 then we have to tilt the integral and calculate B_p using the expression above.

#run it at command line with: 
#visit -assume_format Chombo -cli -s /Users/mhuartee/Box\ Sync/py/polarisation_xmdv_3.py
#meant to produce DATE.okc file to be opened with visit then.

MODE25D=1
MODE3D=0

mode=MODE25D

SetActiveWindow(1)
DeleteActivePlots()

import math, datetime, sys, time, os
#import numpy as np
SuppressMessages(2)
start = time.time()#measure program time


#  *** params ***
angjet=0.0
tt=math.tan(angjet)
cc=math.cos(angjet)
ss=math.sin(angjet)
DefineScalarExpression("angjet", str(angjet))
xcells=16 #from data
pix=(8,8) #36min nx,ny
ylim=(0.,8.)
zlim=(28,36) #ignorde in MODE25D
xlim=(0,8)
dy=(ylim[1]-ylim[0])/float(pix[1])
dx=(xlim[1]-xlim[0])/float(pix[0])
where='/scratch/bliu17/Laurence/run_dir_3d_mhd/out/'
###############################

hoy=datetime.datetime.now().strftime("%m-%d")

#   z=jet                                 
#   ^                                     
#   |                                     
#   |---> y    x axis is pointing out      
DefineScalarExpression("3dB_strength", "sqrt(Bx*Bx+By*By+Bz*Bz)")
DefineScalarExpression("fact", "(px*px+py*py+pz*pz)/rho/2.")#Ekinetic

if mode == MODE3D #Create expressions in advance
     DefineScalarExpression("Bproj", "By*cos(angjet)+Bz*sin(angjet)")
     DefineScalarExpression("I", "fact*(Bx^2+Bproj^2)")
     DefineScalarExpression("Q", "fact*(Bx^2-Bproj^2)")
     DefineScalarExpression("U", "2.*fact*Bproj*Bx")
else if mode == MODE25D
     DefineScalarExpression("xp", "1.0")
     DefineScalarExpression("Bzz", "Bz*xp/x")
     DefineScalarExpression("Bxx", Bz*sqrt(x^2-xp^2)/x")
     DefineScalarExpression("Bproj", "By*cos(angjet)+Bzz*sin(angjet)")
     DefineScalarExpression("dzdx", "x/sqrt(x^2-xp^2)")   
     DefineScalarExpression("I","2*fact*( Bxx^2+Bproj^2) * dzdx")
     DefineScalarExpression("Q","2*fact*( (Bxx^2-Bproj^2)*dxdz")
     DefineScalarExpression("U","2*2.*fact*Bxx*Bproj")    

#get data extents
AddPlot("Pseudocolor", "rho", 1, 0)
DrawPlots()
SetQueryFloatFormat("%g")
Query("SpatialExtents")
xmax=float(GetQueryOutputValue()[1])
#ymax=GetQueryOutputValue()[3]
#zmax=GetQueryOutputValue()[5]
print xmax

intcells=int(
   float(xcells)*(ylim[1]-ylim[0])/float(xmax)
   )
DeleteActivePlots()

#lineout initialize
AddPlot("Curve", "operators/Lineout/I", 1, 0)
DrawPlots()
LineoutAtts = LineoutAttributes()
LineoutAtts.point1 = (0., 0., 32.)
LineoutAtts.point2 = (xmax, 0., 32.)
LineoutAtts.interactive = 0
LineoutAtts.ignoreGlobal = 0
LineoutAtts.samplingOn = 0
LineoutAtts.numberOfSamplePoints = 10
LineoutAtts.reflineLabels = 0
LineoutAtts.interactive = 1
LineoutAtts.numberOfSamplePoints = intcells


#initial IO
lines=pix[0]*pix[1]
filename='/scratch/bliu17/Laurence/run_dir_3d_mhd/polarMaps/'+hoy+'pol.okc'
f = open(filename, 'w+')
print f.write("%g %g %g\n" % (6,lines,6))
print f.write("x\ny\nz\npolx\npoly\npolz\n")
print f.write("%g %g %g \n" % (ylim[0],float(1+pix[0])*dy, 10))
print f.write("%g %g %g \n" % (zlim[0],float(1+pix[1])*dz, 10))
print f.write("%g %g %g \n" % (0.,0.00, 10))
print f.write("%g %g %g \n" % (0.,0.75, 10))
print f.write("%g %g %g \n" % (0.,0.75, 10))
print f.write("%g %g %g \n" % (0.,0.00, 10))

#compute (integral in z)
yshift=.5*tt*(zlim[1]-zlim[0])

x=xlim[0]
while x<=xlim[1]:
   print str(int(98.*(x-xlim[0])/(xlim[1]-xlim[0])))+'%'

   if mode == MODE25D: #update expressions for integral calculations
        DefineScalarExpression("xp",str(x))
  
   y=ylim[0]
   while y<=ylim[1]:
   
      if mode == MODE25D: #update expressions for integral calculations
            LineoutAtts.point1 = (x,y,0)
            LineoutAtts.point2 = (xlim[1],y+dz,0)
      elsif mode == MODE3D:
           LineoutAtts.point1 = (x,y-yshift,zlim[0])
           LineoutAtts.point2 = (x,y+yshift,zlim[1])

      SetOperatorOptions(LineoutAtts, 0)


      #I
      ChangeActivePlotsVar("I")                                             
      Query("Integrate")                                                     
      I0 = GetQueryOutputValue()/float(intcells)                       
      #Q                                                                     
      ChangeActivePlotsVar("Q")                                              
      Query("Integrate")                                                     
      Q0 = 0.75*GetQueryOutputValue()/float(intcells)                       
      #U                                                                     
      ChangeActivePlotsVar("U")                                              
      Query("Integrate")                                                     
      U0 = 0.75*GetQueryOutputValue()/float(intcells)                                                                             
      
      if I0 != 0.:
         degree=( math.sqrt(math.pow(U0,2.) +
            math.pow(Q0,2.))/I0 
            )

      if Q0 != 0.:
         psi=0.5*(math.atan(U0/Q0)+0.*math.pi)

      f.write("%g %g %g %g %g %g\n" % (x,y,z,degree*math.sin(psi),
         degree*math.cos(psi),0.))
      
      
      y=y+dy
   x=x+dx

f.close()
stop = time.time()
print '\n%3.1f min \n' % float((stop-start)/60.)
print 'Done!, the xmdv file is called '+filename+'\n'
#os.system('say "Your calculation is done"')