#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
s=math.tan(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", "Bz*cos(angjet)+Bx*sin(angjet)")
     DefineScalarExpression("I", "fact*(By*By+Bproj*Bproj)")
     DefineScalarExpression("Q", "fact*(By*By-Bproj*Bproj)")
     DefineScalarExpression("U", "2.*fact*Bproj*By")

#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*s*(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("I","2*fact*( (Bz*sqrt(x^2-"+str(x)+"^2)/x)^2+By^2) * x/sqrt(x^2-"+str(x)+"^2)")
        DefineScalarExpression("Q","2*fact*( (Bz*sqrt(x^2-"+str(x)+"^2)/x)^2-By^2) * x/sqrt(x^2-"+str(x)+"^2)")           
        DefineScalarExpression("U","2*2.*fact*By*Bz")    
  
   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"')