C+
C NAME:
C	ExpandSW
C PURPOSE:
C CATEGORY:
C CALLING SEQUENCE:
	program ExpandSW
C INPUTS:
C OUTPUTS:
C CALLS:
C SEE ALSO:
C EXTERNAL:
	external	fncWZ
C INCLUDE:
	include		'openfile.h'
C PROCEDURE:
C MODIFICATION HISTORY:
C	???-????, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	character	cSay*8		/'ExpandSW'/

	parameter	(nMax2D=360*180)
	parameter	(nMax3D=nMax2D*31)

	real		VMap(nMax2D)
	real		DMap(nMax2D)
 	real		GMap(nMax2D)

	real		V3D(nMax3D)
	real		D3D(nMax3D)

	real		Info3D(nMax3D*3)

	real		RR		/0./
	real		dRR		/0.1/

	character	cFileD*80	/' '/
	character	cFileV*80	/' '/
	character	cFile3D*80	/' '/
	character	cFileVV*80	/' '/
	character	cFileDD*80	/' '/
	character	cFmt*7		/'(F11.4)'/

	character	cStr(10)*80

	integer		nPad(4)		/4*0/

	real		rVec(3)

	logical		bOpenFile

	include		't3d_grid_fnc.h'

	call Say('ExpandSW','I','Specify','file with 3D information')

	if (bOpenFile(OPN__STOP .or. OPN__TEXT,iU,cFile3D,iRecl)) then
	    iStr = 1
	    read (iU, '(A)') cStr(iStr)
	    do while (cStr(iStr)(:1) .eq. ';')
		cStr(iStr) = cStr(iStr)(3:)
		iStr = iStr+1
		read (iU, '(A)') cStr(iStr)
	    end do
	    iStr = iStr-1
	    iU = iFreeLun(iU)
	end if

	if (iStr .lt. 4) call Say('ExpandSW','E','Err','Not a 3D info file: '//cFile3D)

	call Str2Flt_Exp(.TRUE.)
	nVec = 2
	call Str2Flt(cStr(2),nVec,rVec)
	CR  = rVec(1)			! Start of rotation
	nCR = rVec(2)			! # rotations

	call Str2Flt_Exp(.TRUE.)
	nVec = 2
	call Str2Flt(cStr(3),nVec,rVec)
	RR  = rVec(1)			! Reference distance
	dRR = rVec(2)			! Radial grid distance

	nVec = 3
	call Str2Flt(cStr(4),nVec,rVec)
	nLng = nint(rVec(1))		! # longitudes
	nLat = nint(rVec(2))		! # latitudes
	nRad = nint(rVec(3))		! # radial distances

	write (*,*) ' Start of rotation: ',CR
	write (*,*) ' Number of rotations: ',nCR
	write (*,*) ' Reference distance: ',RR
	write (*,*) ' Radial grid distance: ', dRR
	write (*,*) ' Dimensions: ',nLng,nLat,nRad

	if (iFltArr(cFile3D,10**6+0,0.,1.,fncWZ,0.,nMax3D*3,nLng,nLat*nRad*3,nPad,Info3D) .eq. 0)
     &		call Say(cSay,'E','Error','reading file '//cFile3D)

	call AskChar('Velocity file at source surface$file',cFileV)

	if (iFltArr(cFileV,10**6+40,0.,1.,fncWZ,BadR4(),nMax2D,nLng,nLat,nPad,VMap) .eq. 0)
     &		call Say(cSay,'E','Error','reading file '//cFileV)

	call AskChar(' Density file at source surface$file',cFileD)

	if (iFltArr(cFileD,10**6+40,0.,1.,fncWZ,BadR4(),nMax2D,nLng,nLat,nPad,DMap) .eq. 0)
     &		call Say(cSay,'E','Error','reading file '//cFileD)

	call Say(cSay,'I','Calculating','3D velocity and density matrices')

	N = nLng*nLat*nRad
	iCshift = 0*N+1
	iVratio = 1*N+1
	iDratio = 2*N+1

	N = nLng*nLat
	do I=1,nRad
	    R = RRhght(I)
	    i3D = (I-1)*N+1
	    call ExtractMaps(R,RR,dRR,nLng,nLat,nRad,VMap,DMap,Info3D(iCshift),Info3D(iVratio),Info3D(iDratio),V3D(i3D),D3D(i3D))

	end do

	P = SUN__RAU*SUN__RAU
	call ArrR4GetMinMax(-N,VMap,vmin,vmax)
	call ArrR4GetMinMax(-N,DMap,dmin,dmax)
	print *, '>',vmin,vmax,dmin*P,dmax*P
	do i=1,nRad
	    i3d = (i-1)*N+1
	    call ArrR4GetMinMax(-N,v3d(i3d),vmin,vmax)
	    call ArrR4GetMinMax(-N,d3d(i3d),dmin,dmax)
	    print *, '<',vmin,vmax,dmin,dmax
	enddo
	print *, ' '

	call AskChar('Output file for velocity matrix$parse',cFileVV)

	cStr(iStr) = 'Velocity (km/s)'

	call WR2DARR(1,nLng,nLat*nRad,V3D,cFileVV,cFmt,.FALSE.,iStr,cStr)

	call AskChar('Output file for density matrix$parse',cFileDD)

	cStr(iStr) = 'Normalized densities (cm^-3*AU^2)'

	call WR2DARR(1,nLng,nLat*nRad,D3D,cFileDD,cFmt,.FALSE.,iStr,cStr)

	end