c
c 	program 3-D proton (2-D electron) velocity distribution function
c
c       some indices and plot variables
c	
	integer i, j, k, l, m, n, o, nx, ny, nz
c
c       anzahl_e             = number of electron data
c       anzahl_p             = number of proton data
c       dim                  = maximal number of data points (protons)
c	earcarrot            = Carrington rotation number from Earth
c       helcarrot            = Carrington rotation number from spacecraft
c       i1a_3_on             = instrument i1a on (=0) or i3 on (=1)
c       index_shift
c       minus                = indicates a HDM file which contained bad
c                              data (frames), but could be handled as NDM 
c                              file 
c	mode                 = HDM or NMD
c	ndim                 = defines scale for electron plot
c       perihelion_shift     = on or of (0 or 1)
c       status               = file handling variable (0 = file exist)
c       win_x, win_y         = define plot windows
c       win_heigth, win_width= define plot windows
c       dummy, yesno         = dummy variables
c	
	integer anzahl_e, anzahl_p, dim, earcarrot,
	1   helcarrot, i1a_3_on, index_shift, minus, mode, ndim,
	1   perihelion_shift,shift, status, dummy, win_x, win_y,
	2   win_width, win_height,yesno
c
c       indexi1.             = position of maximum value of instrument
c	
	integer index_i1a(1), index_i1b(1)
	parameter(dim = 1568)
c
c       bsacle               = scale magnetic vector to fit into plot
c       bbetrag              = modulus of magnetic field
c       earcarlat,earcarlng  = Carrington longitude and latitude from Earth
c       earsdissun           = distance Earth-Sun (always zero, only for
c                              history)
c       helcarlng, helcarlng = Carrington longitude and latitude from Sun
c       helvnorm, helvrad    = tangential (in ecliptic) and radial velocity
c                              in AU / day
c       hseangle             = spacvecraft-Sun-Earth angle
c       nscale               = scaling for electron plot
c       n_p                  = proton density [cm^{-3}]
c       vxs,vys,vzs          = center of mass velocity (1-fluid velocity)
c       x,y, zz, zzz         = dummy parameters used to read in 3D proton
c                              distribiution
c
c       xmin, xmax, xor, xstep  | maximum, minimum origin and stepsize
c	ymax, ymin, yor, ystep  | for plot routines
c       zmax, zmin, zor, zstep  |
c       sinfi,    cosfi       = sine and cosine of azimuth of proton velocity
c       sintheta, costheta    = sine and cosine of elevation of proton velocity
c       background_i1.        = background of instruments (here the first two
c                               and last two channels)
c
	real bscale, bbetrag, earcarlat, earcarlng, eardissun, helcarlat
	1   ,helcarlng, heldissun, helvnorm, helvrad,hseangle
	1   ,nscale, n_p, vxs, vys,vzs, x, y, zz, zzz, sintheta,
	2   costheta,sinfi,cosfi, xmax,xmin, xor, xstep, ymax, ymin, yor
	2   , ystep,zmax, zstep, zmin, zor, zzmax, background_i1a
	3   ,background_i1b
c   
c       b                     = 3 components of magnetic field [nT]
c       bsig                  = stanadrad deviation of the 3 components 
c                               of magnetic field [nT]
c       i1aint                = 32 onboard integrated enegry channels of
c                               protons insrtument i1a
c       i1bint                = 32 onboard integrated enegry channels of
c                               protons insrtument i1b
c       i1av, iabv            = corresponding velocities to i1aint, i1bint
c       i1abv, i1int          = sum and ratio of i1aint, i1bint (must be
c                               improved)
c       maplh                 = 3 dimenional array: density [cm^{-3}], 
c                                velocity [km/s] and temperature [K] of alpha
c	                         particles
c	mi1b                  = 3 dimenional array: density [cm^{-3}], 
c                                velocity [km/s], temperature [K]
c                                of protons
        real b(3), bsig(3), i1aint(32),
	1   i1av(32), i1bint(32), i1bv(32), i1abv(32), i1int(32),
	1   malph(3), mi1a(5), mi1b(3)
c
c       imat, zmat, wmat       = used for contour plots
c	
	integer, allocatable, dimension(:,:) :: imat
	real,    allocatable, dimension(:,:) :: zmat, wmat
c
c
c       vx, vy, vz            = velocity components corresponding
c	                        to respective distribution element
c                               1-dimensional arrays (protons) [km/s]
c       z                     = value of distribution function
c                               1-dimensional arrays (protons)
c                       	[cm^[-6] s^[-3])]
c       vxe, vye, ze          = 2-D electron parameters (see above)
c	
	real vx(dim), vy(dim), z(dim), vz(dim),
	1   vxe(dim), vye(dim), ze(dim)
c
c      
c       opens the file 
c	
 	open(16, file = 'xyz_file')
c
c       now read in the data. The first data set, the correction from
c       Pizzo et al. and for the finite solar diameter will not be used
c       (the first three records).
c
c	All other parameters are read in. For a definition see top.
c
c       A 'read' command without a parameter reads of a comment in the
c       data file 	
c	
	read(16,*)
        read(16,*)
        read(16,*)
	read(16,*,err=100) mode,shift,perihelion_shift, minus, i1a_3_on
	read(16,*,err=101) heldissun,helcarlng, helcarlat, helcarrot
	read(16,*,err=102) eardissun, earcarlng,earcarlat,
	1   hseangle, earcarrot
	read(16,*,err=103) helvrad, helvnorm
	read(16,*,err=104) mi1a
	read(16,*,err=105) malph
	read(16,*,err=106) mi1b
	read(16,*,err=107) b, bsig
	read(16,*) 
	read(16,*,err=108) i1aint
	read(16,*) 
	read(16,*,err=109) i1av
	read(16,*) 
	read(16,*,err=110) i1bint
 	read(16,*) 
	read(16,*,err=111) i1bv
	read(16,*)
	read(16,*)
c
c       Define velocity with respect to center of mass (1-D velocity).
c       Transform from degree to radian. Calculate center of mass
c       velocity 'v.s'
c	
	i1av = i1av - mi1a(2)
	i1bv = i1bv - mi1a(2)
	sinfi = sin(mi1a(4) * pi180)
	cosfi = cos(mi1a(4) * pi180)
	sintheta = sin(mi1a(5) * pi180)
	costheta = cos(mi1a(5) * pi180)
	vxs = mi1a(2) * cosfi * costheta
	vys = mi1a(2) * sinfi * costheta
	vzs = mi1a(2) * sintheta
c
c       set to predefined values (z,vx,vy,vz are arrays)
c	
	anzahl_p = 0 ; zmax = -1.e10 ; z = 0 ; vx = 0. ; vy = 0. ; vz = 0.
c
c       now loop through the distribution fucntion. First read in the
c       the indices i,j,k of distribution function and the velocity
c       components. The dummy variables zz, x,y,zzz contain the value of
c       the distrubution function, and the repective velocity components.
c       The variable 'dummy' contains the count rate, which can be
c       used for statistics. The dummy variables x,y,zzz are then 
c       transformed in the center of mass system and written to the
c       arrays zz,vx,vy,vz. If an error occurs during read, the program
c       will continue at label 1 ('err=1'). The error is forced by a
c       comment line in the data file. The next data will be the
c       electrons, which are handled later.
c	
 2	anzahl_p = anzahl_p + 1
	read(16,*,err=1) i,j,k,zz,dummy,x,y,zzz
	z(anzahl_p) = zz
	vx(anzahl_p) = x - vxs
	vy(anzahl_p) = y - vys
	vz(anzahl_p) =zzz - vzs
	go to 2
c
c       end loop reading proton distribution
c       
c
c============================ electrons ===========================
c
c       Now read electron data. Procedure similar to that of protons.
c       First skip one coment line, then loop, Structure: i,j indices
c       of distribution function, zz value of distribution function,
c       dummy count rate, x,y electron velcoities
c       
c	
	read(16,*)
 3	read(16,*, end=4) i,j,zz,dummy,x,y
	ze(anzahl_e) = zz ; vxe(anzahl_e) = x ; vye(anzahl_e) = y
	anzahl_e = anzahl_e + 1
	go to 3
c
c       end electron loop
c
c       
c       close(16,*)
        stop  
        end