; $Id: read_is.pro,v 1.1 2010/06/29 15:28:06 cooper Exp $

;PRO read_is
;read insitu instrument data for output with jmap_in_situ
;INPUT:instrument-name of instrument(ST_A/B,WIND,OMNI,ACE)
;	s_time-earliest date to be retrieved
;	e_time-latest(most recent) date to be retrieved
;KEYWORDS:gsm=use (GSM) geocentric solar magnetospheric coordinate system to calculate phi and 
;		theta (output still in RTN) [This only applies to omni/wind, all other only have GSE]
;		 the default is GSE coordinates(geocentric solar equatorial)
;	getace=uses the get_acedata function to get ace data, rather than reading saved files
;		this is  slower, but has less data gaps
;OUTPUT FORMAT: structure with two tags, data and magdata.
;	data is array:		[time,speed,density,temperature]
;	magdata is an array:	[time,btotal,phi(lon),theta(lat),beta] (phi,theta in RTN)
;	beta has corrections to approximately include alpha particle and electron pressures
;SYNTAX: read_is,instrument,start_time,end_time

; $Log: read_is.pro,v $
; Revision 1.1  2010/06/29 15:28:06  cooper
; Version 1
;

FUNCTION read_is, instrument, s_time, e_time, gsm=gsm, getace=getace
	is_path=getenv('INSITU_PATH')
	IF instrument eq 'ST_A' OR instrument eq 'ST_B' THEN BEGIN
		IF instrument eq 'ST_A' THEN ins='A' else ins='B'
		;read plastic
		s=anytim2utc(s_time,/ccsds) ;make dates same format as the input file,
		e=anytim2utc(e_time,/ccsds) ;for string comparisons, quicker than converting
		syear=strmid(s,0,4) & sy=fix(syear)
		eyear=strmid(s,0,4) & ey=fix(eyear)
		readfile=concat_dir(concat_dir(is_path,instrument,/dir),'plastic'+ins+syear+'.txt')
		OPENR,lun,readfile,/get_lun
		hdr=''
		readf,lun,hdr
		i=0L
		temp={date:'',dens:0.,vel:0.,t:0.,vx:0.,vy:0.,vz:0.}
		data=replicate(temp,550000)
		WHILE ~eof(lun) DO BEGIN
			readf,lun,temp
			data[i]=temp
			i++
		ENDWHILE
		free_lun,lun
		ind=where(data.date GE s AND data.date LE e AND data.dens GT 0)
		data1=data[ind]
		records=[[utc2tai(data1.date)],[data1.vel],[data1.dens],[data1.t]]
		records=transpose(records)
		;read impact/mag
		readfile=concat_dir(concat_dir(is_path,instrument,/dir),'impact'+ins+syear+'.txt')
		IF ~file_test(readfile) THEN mess=dialog_message('You do not have impact data for this time. See reformat_impact\.pro for instructions on getting impact data')
		OPENR,lun,readfile,/get_lun
		hdr=''
		readf,lun,hdr
		i=0L
		temp={date:'',br:0.,bt:0.,bn:0.}
		magdata=replicate(temp,550000)
		WHILE ~eof(lun) DO BEGIN
			readf,lun,temp
			magdata[i]=temp
			i++
		ENDWHILE
		free_lun,lun
		ind=where(magdata.date GE s AND magdata.date LE e AND magdata.br LT 10.^5)
		magdata=magdata[ind]
		;find lon/lan (in RTN)
		;find the indices for each quadrant and calibrate accordingly
		br=magdata.br & bn=magdata.bn & bt=magdata.bt
		magb=sqrt(br^2+bn^2+bt^2)
		theta=atan(bn/sqrt(br^2+bt^2))*180./!pi
		phi0=abs(atan(bt/br))
		phi=dblarr(n_elements(phi0))
		;already in RTN coords
		q1=where(br GE 0 and bt GE 0) & IF q1[0] ne -1 THEN phi[q1]=phi0[q1]
		q2=where(br LT 0 and bt GE 0) & IF q2[0] ne -1 THEN phi[q2]=!pi-phi0[q2]
		q3=where(br LE 0 and BT LT 0) & IF q3[0] ne -1 THEN phi[q3]=!pi+phi0[q3]
		q4=where(br GT 0 and BT LT 0) & IF q4[0] ne -1 THEN phi[q4]=2*!pi-phi0[q4]
		phi*=180./!pi
		;get beta, using both sets of data, need same times, use the mag ind on plastic data
		beta=(data[ind].dens*(1.16*data[ind].t+1.08*130000.)*1.38*10.^(-6))/$
			(magb^2/(8.*!pi))
		magrecords=[[utc2tai(magdata.date)],[magb],[theta],[phi],[beta]]
		magrecords=transpose(magrecords)
	ENDIF
	IF instrument eq 'ACE' THEN BEGIN
		IF keyword_set(getace) THEN BEGIN
			s_time=anytim2utc(s_time)
			e_time=anytim2utc(e_time)
			data=get_acedata(s_time,e_time,/swepam) ;solar wind data
			days=data[*].mjd
			times=data[*].time
			dates=utc2tai({mjd:days,time:times})
			;dates from mag instrument are different, so need to store separately
			magdata=get_acedata(s_time,e_time,/mag) ;magnetic data
			magdates=utc2tai({mjd:magdata[*].mjd,time:magdata[*].time})
			records=[[dates],[data[*].b_speed],[data[*].p_density],[data[*].ion_temp]]
			records=transpose(records)
			;To calculate beta, find times that match most closely and find 
			;plasma pressure over magnetic pressure
			;Need this loop because the times dont agree, but its not very fast. better way?
			n=n_elements(magdates)
			beta=fltarr(n)
			n-=1
			for i=0L,n do begin	;is there a faster way to make them match?
				dummy=min(dates-magdates[i],d,/absolute)
				beta[i]=(data[d].p_density*(1.16*data[d].ion_temp+1.08*130000.)*1.38*10.^(-6))/$
					((magdata[i].bt)^2/(8.*!pi)) ;in cgs units
			endfor
			magrecords=[[magdates],[magdata[*].bt],[magdata[*].lat],[magdata[*].long],[beta]]
			magrecords=transpose(magrecords)
		ENDIF ELSE BEGIN ;use stored data
			s=anytim2utc(s_time,/ccsds)
			e=anytim2utc(e_time,/ccsds)
			syear=strmid(s,0,4)
			readfile=concat_dir(concat_dir(is_path,instrument,/dir),'ace'+syear+'.txt')
			OPENR,lun,readfile,/get_lun
			hdr=''
			readf,lun,hdr
			temp={date:'',dens:0.,vel:0.,t:0.,bx:0.,by:0.,bz:0.,ratio:0.}
			data=replicate(temp,550000)
			i=0L
			WHILE ~eof(lun) DO BEGIN
				readf,lun,temp
				data[i]=temp
				i++
			ENDWHILE
			ind=where(data.date GE s and data.date LE e and data.dens NE 999.99)
			data=data[ind]
			dates=utc2tai(data.date)
			records=[[dates],[data.vel],[data.dens],[data.t]]
			records=transpose(records)
			bx=data.bx & by=data.by & bz=data.bz
			btotal=sqrt(bx^2+by^2+bz^2)
			theta=180./!pi*atan(bz/SQRT(bx^2+by^2))
			phi0=abs(atan(by/bx))
			phi=dblarr(n_elements(phi0))
			;already in rtn coords
			q1=where(bx LE 0 AND by GE 0) & phi[q1]=phi0[q1]
			q2=where(bx GT 0 AND by GE 0) & phi[q2]=!pi-phi0[q2]
			q3=where(bx GE 0 AND by LT 0) & phi[q3]=!pi+phi0[q3]
			q4=where(bx LT 0 AND by LT 0) & phi[q4]=2*!pi-phi0[q4]
			phi*=180./!pi
			beta=(data.dens*(1.16*data.t+1.08*130000.)*1.38*10.^(-6))/$
				(btotal^2/(8.*!pi))
			magrecords=[[dates],[btotal],[theta],[phi],[beta]]
			magrecords=transpose(magrecords)
		ENDELSE
	ENDIF
	IF instrument eq 'OMNI' OR instrument eq 'WIND' THEN BEGIN
		s=anytim2utc(s_time,/ccsds)
		e=anytim2utc(e_time,/ccsds)
		syear=strmid(s,0,4)
		temp={date:'',dens:0.,vel:0.,t:0.,bx:0.,bygse:0.,bzgse:0.,bygsm:0.,bzgsm:0,vx:0.,vy:0.,vz:0.}
		data=replicate(temp,550000)
		readfile=concat_dir(concat_dir(is_path,instrument,/dir),strlowcase(instrument)+syear+'.txt')
		OPENR,lun,readfile,/get_lun
		i=0L
		hdr=''
		readf,lun,hdr
		WHILE ~eof(lun) DO BEGIN
			readf,lun,temp
			data[i]=temp
			i++
		ENDWHILE
		free_lun,lun
		ind=where(data.date GE s and data.date LE e AND data.dens NE 999.99 AND data.bx NE 9999.99)
		data=data[ind]
		dates=utc2tai(data.date)
		records=[[dates],[data.vel],[data.dens],[data.t]]
		records=transpose(records)
		bx=data.bx
		IF keyword_set(gsm) THEN BEGIN
			print,'USING GSM'
			by=data.bygsm & bz=data.bzgsm
		ENDIF ELSE BEGIN
			by=data.bygse & bz=data.bzgse
		ENDELSE
		btotal=SQRT(bx^2+by^2+bz^2)
		theta=180./!pi*atan(bz/sqrt(bx^2+by^2))
		phi0=abs(atan(by/bx))
		phi=dblarr(n_elements(phi0))
		;this is in gse/gsm, want to plot it in rtn though
		q1=where(bx LE 0 AND by GE 0) & phi[q1]=phi0[q1]
		q2=where(bx GT 0 AND by GE 0) & phi[q2]=!pi-phi0[q2]
		q3=where(bx GE 0 AND by LT 0) & phi[q3]=!pi+phi0[q3]
		q4=where(bx LT 0 AND by LT 0) & phi[q4]=2*!pi-phi0[q4]
		phi*=180./!pi
		beta=(data.dens*(1.16*data.t+1.08*130000.)*1.38*10.^(-6))/$
			(btotal^2/(8.*!pi))
		magrecords=[[dates],[btotal],[theta],[phi],[beta]]
		magrecords=transpose(magrecords)
	ENDIF
return,{data:records,magdata:magrecords}
END