import zeep
import numpy as np

client = zeep.Client('http://turbulence.pha.jhu.edu/service/turbulence.asmx?WSDL')
ArrayOfFloat = client.get_type('ns0:ArrayOfFloat')
ArrayOfArrayOfFloat = client.get_type('ns0:ArrayOfArrayOfFloat')
SpatialInterpolation = client.get_type('ns0:SpatialInterpolation')
TemporalInterpolation = client.get_type('ns0:TemporalInterpolation')

token="edu.jhu.pha.turbulence.testing-201406" #replace with your own token

nnp=5 #number of points
points=np.random.rand(nnp,3)

# convert to JHTDB structures
x_coor=ArrayOfFloat(points[:,0].tolist())
y_coor=ArrayOfFloat(points[:,1].tolist())
z_coor=ArrayOfFloat(points[:,2].tolist())
point=ArrayOfArrayOfFloat([x_coor,y_coor,z_coor]);

print(points)

"""
Explanation: pyJHTDB are failed to compile on windows. One alternative way might be to use zeep package.  
More details can be found at http://turbulence.pha.jhu.edu/service/turbulence.asmx
End of explanation
"""


Function_name="GetVelocityGradient" 
time=0.6
number_of_component=9 # change this based on function_name, see http://turbulence.pha.jhu.edu/webquery/query.aspx

result=client.service.GetData_Python(Function_name, token,"isotropic1024coarse", 0.6, 
                                     SpatialInterpolation("None_Fd4"), TemporalInterpolation("None"), point)
result=np.array(result).reshape((-1, number_of_component))

print(result)

"""
Explanation: In GetData_Python, Function_name could be
GetVelocity,            GetMagneticField,          GetVectorPotential,  
GetVelocityGradient,    GetMagneticFieldGradient,  GetVectorPotentialGradient,  
GetVelocityHessian,     GetMagneticHessian,        GetVectorPotentialHessian,  
GetVelocityLaplacian,   GetMagneticFieldLaplacian, GetVectorPotentialLaplacian,    
GetPressure,            GetTemperature,            GetDensity,  
GetPressureGradient,    GetTemperatureGradient,    GetDensityGradient,  
GetPressureHessian,     GetTemperatureHessian,     GetDensityHessian,     
GetVelocityAndPressure, GetVelocityAndTemperature, GetForce, GetInvariant
End of explanation
"""


Function_name="GetPosition"
startTime=0.1
endTime=0.2
dt=0.02
number_of_component=3 # change this based on function_name, see http://turbulence.pha.jhu.edu/webquery/query.aspx

result=client.service.GetPosition_Python(Function_name, token,"isotropic1024coarse", startTime, endTime, dt, 
                                         SpatialInterpolation("None"), point)
result=np.array(result).reshape((-1, number_of_component))
print(result)

"""
Explanation: In GetPosition_Python, Function_name could be
GetPosition only.
End of explanation
"""


Function_name="GetBoxFilter" #could also be
field="u"
time=0.6
filterwidth=0.05
spacing=0 #spacing is only used in GetBoxFilterGradient, but always provide it.
number_of_component=3 # change this based on function_name, see http://turbulence.pha.jhu.edu/webquery/query.aspx

result=client.service.GetFilter_Python("GetBoxFilter",token,"isotropic1024coarse", field, 
                                       time, filterwidth, SpatialInterpolation("None"), point, spacing) 
result=np.array(result).reshape((-1, number_of_component))
print(result)

import struct
import base64

field="u"
timestep=1
x_start=1
y_start=1
z_start=1
x_end=2
y_end=5
z_end=8
x_step=1
y_step=1
z_step=1
filter_width=0
    
result=client.service.GetAnyCutoutWeb(token,"isotropic1024coarse", field, timestep,
                                       x_start, y_start, z_start, x_end, y_end, z_end,
                                      x_step, y_step, z_step, filter_width, "")  # put empty string for the last parameter

# transfer base64 format to numpy
number_of_component=3 # change this based on the field
nx=len(range(x_start, x_end+1, x_step))
ny=len(range(y_start, y_end+1, y_step))
nz=len(range(z_start, z_end+1, z_step))
base64_len=int(nx*ny*nz*number_of_component)
base64_format='<'+str(base64_len)+'f'

result=struct.unpack(base64_format, result)
result=np.array(result).reshape((nz, ny, nx, number_of_component))
print(result.shape)  # see the shape of the result and compare it with nx, ny, nz and number of component

"""
Explanation: In GetFilter_Python, Function_name could be
GetBoxFilter, GetBoxFilterSGSscalar, GetBoxFilterSGSvector,   
GetBoxFilterSGSsymtensor, GetBoxFilterSGStensor, GetBoxFilterGradient.
End of explanation
"""


