PRO DATACUBE,r,theta,phi,ModPramsStruct,ModSolStruct,nreinit=nreinit

;+
; Name:
;      DATACUBE
;
; Purpose:
;      Main wrapper program to call the data cube and get back out the
;      parameters of the cube: pressure, density, temp, br, bth, bphi
;
; Calling sequence:
;      DATACUBE,r,theta,phi,ModPramsStruct, ModSolStruct 
;
; Input:
;
;      r : scalar x1 (r) position of the desired interpolation point
;      in units of solar radii
;
;      theta : scalar x2 (theta) position of the desired interpolation
;      point in radians
;
;      phi : scalar x3 (phi) position of the desired interpolation
;      point in radians
;
;      ModPramsStruct : Structure created further up the chain by
;      datacubeprams.pro that includes all the necessary modifications to
;      the cube like rotating the central meridian and latitude, and
;      what to do outside the cube.
;
; Output:
;      ModSolStruct : interpolated value of the model parameters at
;      the given r,theta,phi
;
; Author and history:
;      Laurel Rachmeler: original version
;       modified for hydrostatic model consistency SEG March 2011
;		Vectorized, TAK 8-Feb-2013
;-


COMPILE_OPT IDL2 ;default long and square brackets for array subscripts

;**read in the values from ModPramStruct**

        cuberot=ModPramsStruct.Cuberot

;deal with outside of cube

;***UNWRAP THE POSITION DATA***
;The position input might be a point, a line, or a
;cube. I want them to be a long line, so unwrap if it isn't

;If input arrays have mare than 1D, reform them into 1D arrays

	CoordSize=SIZE(r)	
	NDim=CoordSize[0]
	nunwrapall=PRODUCT(CoordSize[1:NDim]) ;map input onto 1-D 
                                              ;arrays of dimension 
                                              ;nunwrap
        nunwrapall=LONG(nunwrapall)
	p1=REFORM(r,nunwrapall)
	p2=REFORM(theta,nunwrapall)
	p3=REFORM(phi,nunwrapall)
        ;at this point
        ;p1=r
        ;p2=theta
        ;p3=phi
        ;in original passed in coordinates
 
        IF cuberot THEN BEGIN
;***TRANSFORM TO CARTESIAN, ORIGINAL COORDINATES, POSITION ONLY***
            FOR_SPHERETOCART, p1, p2, p3
            ;p1=x
            ;p2=y
            ;p3=z
            ;in original passed in coordinates
        
;***DO THE CUBE ROTATION***
            FOR_CUBEROT, cuberot, p1, p2, p3
            ;p1=x
            ;p2=y
            ;p3=z
            ;in x-rotated coordinates

;***TRANSFORM BACK TO SPHERICAL, ROTATED COORDINATES***
            FOR_CARTTOSPHERE, p1, p2, p3
            ;p1=r
            ;p2=theta
            ;p3=phi
            ;in x-rotated coordinates

        ENDIF
 
;***DO THE INTERPOLATION FROM THE DATA CUBE***
;the actual cube is passed in with ModPramsStruct.cube
        
;only restore the cube the first time. Otherwise, use the common block
;where the cube is stored

IF nreinit EQ 1 THEN BEGIN
   print, 'datacube.pro: starting cube restore'
   RESTORE, ModPramsStruct.cubename
   print, 'datacube.pro: done with cube restore'
   COMMON THECUBE, cube
ENDIF 

;I am removing the main look, but we should still have a limit so we don't overload memory
;I think this program my be called in side a loop, so that would be really unneeded.
BiteSize=30000
npts=n_elements(p1)
NLoops=fix(npts/BiteSize)-1
if (npts mod BiteSize) gt 0 then NLoops=NLoops+1
a_pres=p1*0d0
a_dens=p1*0d0
a_temp=p1*0d0
a_b1=p1*0d0
a_b2=p1*0d0
a_b3=p1*0d0

for i=0,NLoops do begin
   FOR_INTERP_CUBE, ModPramsStruct,cube,$		;this might still have to be a loop, but in bigger chuncks
                p1[i*BiteSize:((i+1)*BiteSize-1)<(Npts-1)],p2[i*BiteSize:((i+1)*BiteSize-1)<(Npts-1)],p3[i*BiteSize:((i+1)*BiteSize-1)<(Npts-1)],$
                pres,dens,temp,br,bth,bph
    a_pres[i*BiteSize:((i+1)*BiteSize-1)<(Npts-1)]=pres
    a_dens[i*BiteSize:((i+1)*BiteSize-1)<(Npts-1)]=dens
    a_temp[i*BiteSize:((i+1)*BiteSize-1)<(Npts-1)]=temp
    a_b1[i*BiteSize:((i+1)*BiteSize-1)<(Npts-1)]=br
    a_b2[i*BiteSize:((i+1)*BiteSize-1)<(Npts-1)]=bth
    a_b3[i*BiteSize:((i+1)*BiteSize-1)<(Npts-1)]=bph
ENDFOR
    ModSolStruct={Pres:a_pres,Dens:a_dens,Temp:a_temp,$
                      Br:a_b1,Bth:a_b2,Bph:a_b3}

        ;p1=r
        ;p2=theta
        ;p3=phi
        ;b1=br
        ;b2=btheta
        ;b3=bphi
        ;in x-rotated coordinates
        
;***TRANSFORM POSITION AND VECTOR TO CARTESIAN, ROTATED COORDINATES***
        IF cuberot THEN BEGIN
            FOR_SPHERETOCART, p1, p2, p3, a_b1, a_b2, a_b3
            ;p1=x
            ;p2=y
            ;p3=z
            ;b1=bx
            ;b2=by
            ;b3=bz
            ;in x-rotated coordinates
        
;***DO THE DE-ROTATION ABOUT THE X-AXIS***

            FOR_CUBEROT, -1*cuberot, p1, p2, p3, a_b1, a_b2, a_b3
            ;p1=x
            ;p2=y
            ;p3=z
            ;b1=bx
            ;b2=by
            ;b3=bz
            ;in original passed in coordinates

;***TRANSFORM POSITION AND VECTOR TO SPHERICAL, ORIGINAL COORDINATES***

            FOR_CARTTOSPHERE, p1, p2, p3, a_b1, a_b2, a_b3
            ;p1=r
            ;p2=theta
            ;p3=phi
            ;b1=br
            ;b2=btheta
            ;b3=bphi
            ;in original passed in coordinates
        ENDIF

;***RE-WRAP THE POSITION AND VECTOR DATA***
   	a_pres=REFORM(a_pres,CoordSize[1:NDim],/overwrite)
   	a_dens=REFORM(a_dens,CoordSize[1:NDim],/overwrite)
   	a_temp=REFORM(a_temp,CoordSize[1:NDim],/overwrite)
   	a_b1=REFORM(a_b1,CoordSize[1:NDim],/overwrite)
  	a_b2=REFORM(a_b2,CoordSize[1:NDim],/overwrite)
  	a_b3=REFORM(a_b3,CoordSize[1:NDim],/overwrite)

        ;don't need to reform positions because I copied the 
        ;original positions from r, theta, phi to p1,p2,p3
        ;r=REFORM(r,CoordSize[1:NDim],/overwrite)
   	;theta=REFORM(theta,CoordSize[1:NDim],/overwrite)
   	;phi=REFORM(phi,CoordSize[1:NDim],/overwrite)

	Vel=a_pres*0.d0

        ModSolStruct={Pres:a_pres,Dens:a_dens,Temp:a_temp,Vel:vel,$
                      Br:a_b1,Bth:a_b2,Bph:a_b3}

END