Changes between Version 3 and Version 4 of PolarizationMaps


Ignore:
Timestamp:
05/02/18 09:09:45 (8 years ago)
Author:
Jonathan
Comment:

—

Legend:

Unmodified
Added
Removed
Modified

  • PolarizationMaps

    v3 v4  
    88Each pixel will be given by the integral of the function along the path
    99
    10 $s = x_i, y_i+s \sin(\theta), z_i+s \cos(\theta)$
     10$s = x_i, y_i+s \sin(\theta), s \cos(\theta)$
    1111
    12 where $\delta y = \cos(\theta) \delta x$ and $\delta z =\sin(\theta) \delta x$
     12where $\delta y = \cos(\theta) \delta x$
    1313
    1414
    15 For 3D we just need to create the expressions in visit, create the expressions and the lineouts, and then integrate the resulting query
     15For 3D we just need to create the expressions in visit, create the lineouts, and then integrate the resulting query
    1616
    1717
    1818For 2.5D, we need to transform the integral along $s$ into an integral in the $xy$ plane.
     19
     20the value at $x_i,y_i+s \sin(\theta), s \cos(\theta)$ is the same as the value rotated into the xy plane at
     21
     22$\left [ \sqrt{x_i^2+\left(s \cos(\theta) \right ) ^2}, y_i+s \sin(\theta), 0 \right ]$
     23
     24or at
     25
     26$X'$
     27
     28where
     29
     30$X' = \left [ x', y_i+s \sin(\theta), 0 \right ]$
     31
     32and
     33
     34$x'=\sqrt{x_i^2+(s \cos (\theta))^2}$
     35
     36so we can write the integral as
     37
     38$\int F(x'(s)
     39
     40and we have $ds = \frac{\sqrt{(x_i^2+(s \cos(\theta))^2}}{s \cos(\theta)} dx + \frac{1}{\sin(\theta)} dy$
     41
     42which under the substitution  gives
     43
     44$ds = \frac{x'}{\sqrt{x'^2-x_i^2}} dx'$
     45
    1946
    2047
    … …  
    3360If $\theta /= 0$ then we have to tilt the integral and calculate $B_p$ using the expression above.
    3461
    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 }}}