ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                 Fri Aug  1 16:40:59 CEST 2003                         c
c        Michael Heuer, Max Planck Institut fuer Aeronomie              c
c                Katlenburg-Lindau 37191, Germany                       c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      program helios
c
      USE PORTLIB
      include 'types.inc'
c     
c     CEM calibration data see blue book (BMFT-FB-W 81-015)
c     
	logical sq
	integer dday(1),idog,kdog,iday,ilast
      include 'calibration.inc' 
c
c     
c     			
c=========================================================c
c     Initialize program                                  c
c=========================================================c
      write(*,*)'From which day do you want to start?'
	read(5,*)iday
	write(*,*)'To which day do you want to end?'
 	read(5,*)ilast

	azimuth    = azimuth   *pi180 
      elevation  = elevation *pi180 
      azimuth3   = azimuth3  *pi180 
      elevation3 = elevation3*pi180 
      azimuthe   = azimuthe  *pi180 
c	do 150,idog=1,1
c	     dday(idog)=idog+51
c150	 continue
c       kdog=1
 1000 format(1x,i2,1x,i1,1x,i2,1x,1pe15.8,1x,i5,6(1x,1pe15.8))
 1001 format(1x,i2,1x,i2,1x,1pe15.8,1x,i5,2(1x,1pe15.8))
c     
c     define letters for output file s = sec, m= min, h = hour
      out(1:1) = 'h'
      out(4:4) = 'm'
      out(7:7) = 's'
c           
c     default values for input file and directory	
      file = 'h1  d   .txt'	
      directory = 'Y  D   '
c
c     reads raw data files, start in directory where
C     these files are located. The raw data files
c     are provided by the MPAe compressed on 3 CD's
c
c==========================================================c
c     Input parameter: spacecraft id, year day of year     c 
c==========================================================c     
c     
      write(6,*)                
      write(6,*)'sc'                
c     $     'spacecraft ? (= 1,2 other number will cancel the program)'
      read(5,*) sc
c	sc=2
c===========================================================================c
c     If wrong spacecraft id stop.                                          c
c     May also be used to halt the program after one day was processed,     c
c     otherwise enter new input parameter (if spacecraft id changes,        c
c     This must be done in the previous statement.) 	                    c
c===========================================================================c     
      if ((sc .lt. 1) .or. (sc .gt. 2)) stop
c
c     ------ get year and day of year ------
      write(6,*) 'year: yy?'
      read(5,'(a2)') year 
c 11    	year='77'
c     
c     ------ check if correct year 1974 to 1986 ------
4      write(*,*)  
11      read(year,'(i2)') idummy  
      if ((idummy .lt. 74) .or. (idummy .gt. 86)) go to 11 
c     
c 12   write(6,*) 'day of the year: ddd (with 0: e.g., 002 for day 2)'
c      read(5,'(a3)') doy
12	write(doy,'(i3)') iday 
	if (iday .lt. 100) doy(1:1)='0'
 	if (iday .lt. 10) doy(2:2)='0'
 	write(*,*)iday
	iday=iday+1
	if (iday .gt.(ilast+1)) stop	
      read(doy,'(i3)') idummy   
c     
c     ------ check day of year ------
      if ((idummy .lt. 1) .or. (idummy .gt. 366)) go to 12 
      write(file(2:2),'(i1)') sc 
      file(3:4) = year
      file(6:8) = doy 
c     
c=========================================================================c
c     Make directory                                                      c
c     At this stage the directories are named with a capital letter Y     c
c     for the year followed by a 2 digit value (e.g., 75 for 1975),       c
c     next comes a capital D and a 3 digit value for the day of the       c
c     year. The directory is made from the input parameters.              c
c=========================================================================c    
      directory(2:3) = year
      directory(5:7) = doy     
c      call system('mkdir '//directory,status) 
	status = SYSTEM('mkdir '//directory )
      if (status .ne. 0) then  
         write(6,*)
     $        'The directory exists, the existing files are overwritten'
         write(6,*) 'Do you want to cancel (yes / no)'
         read(5,'(a3)') yesno
         if ((yesno(1:1) .eq. 'y') .or. (yesno(1:1) .eq. 'Y')) stop
      end if
      open (15, file=file)
	open(17, file = directory//'/'//'baddata.txt')
	open(18, file = directory//'/'//'gooddata.txt')
c     
c=======================================================================c
c     Begin read raw data                                               c
c     See outpout_file_documentation.txt for the different parameter    c
c=======================================================================c     
 3    read(15,*,end=4) yr_doy, msec 
      read(15,*,err=4) bitrate, mode, fformat, dismod, shift,
     $     perihelion_shift,i1a_3_on, i2a_b 
c     ------ skipped: spinrate, pitch, aspect ------
      read(15,*)
      read(15,*,err=4) hellngasc, heldissun, helvrad, helvnorm,
     $ helcarlng,
     $     helcarlat, helcarrot
c dog      read(15,*) eardissun, earcarlng,earcarlat, hseangle, earcarrot
      read(15,10,err=4) eardissun, earcarlng,earcarlat, hseangle,
     $	 earcarrot
10	format(e15.8,e15.8,e15.8,e15.8,i5)
	read(15,*,err=4) zerorate_i1a, zerorate_i1b, zerovar_i1b
      read(15,*,err=4) init_words
      read(15,*,err=4) q_words
      read(15,*,err=4) gmaxel, gmaxaz, gmaxen, gmass
      read(15,*,err=4) mi1a
      read(15,*,err=4) malph
      read(15,*,err=4) mi1b
      read(15,*,err=4) b, bsig
      read(15,*,err=4) i1b
      read(15,*,err=4) i1a_3_int
      read(15,*,err=4) i2b
      read(15,*,err=4) i1a_3
c     
c     ------ subtract zero count rates ------
      do i = 1, 32
         if (i1a_3_int(i) .lt. zerorate_i1a) then
            i1a_3_int(i) = 0
         else
            i1a_3_int(i) = i1a_3_int(i) - NINT(zerorate_i1a)
         endif
         if (i1b(i) .lt. zerorate_i1b) then
            i1b(i) = 0
         else
            i1b(i) = i1b(i) - NINT(zerorate_i1b)
         endif
      end do
c
c     ------ try to recover the initial and quality bits           ------
c     ------ here a comment is at place:                           ------
c     ------ the reconstruction of the binary words is principally ------
c     ------ possible, but the init and qw have only been read     ------
c     ------ properly for NDM. in case of HDM both fields contain  ------
c     ------ corrupt data (cf. counts_ks.for)                      ------
c     
c     ------ clear all data ------
      i1a_3_all = 0             !  data array
      fd1d1a_3 = 0.d0           !  integrated 1d i1a/i3 distribution
      fd1d1b = 0.d0             !  integrated 1d i1b distribution
      fd3d  = 0.d0              !  3d i1a/i3 distribution
      vx = 0.d0                 !  - self explaining -
      vy = 0.d0                 !  - self explaining -
      vz = 0.d0                 !  - self explaining -
      electron_distribution = 0.d0
      hdm_ndm = 0
      eminus = 0 
c==================================c
c     Start filename procedure     c
c==================================c
c
c     ------ write hour, minutes, seconds in the filename and open ------
      th = int(msec * millisec2hour) 
      tm = int((msec * millisec2hour - th)*60)
      ts = int(((msec * millisec2hour - th)*60 - tm)*60)
      out(2:3) = out_0(th)
      out(5:6) = out_0(tm)
      out(8:9) = out_0(ts)        
c 
c     ------ include shift in file name ------
      write(shiftc,'(i1)') shift 
c
c     ------ include data mode in filename: ndm or hdm ------
      select case (mode)
      case (0)
         endung = 'ndm'
      case (1)
         endung = 'hdm'
      end select
c
c     ------ open output file ------
      open(16, file = directory//'/'//file(1:2)//'y'//file(3:8)/
     $     /out//'_'//endung//'.'//shiftc)
c     
c      write(6,*) directory//'/'//file(1:2)//'y'//file(3:8)//out/
c     $     /'_'//endung//'.'//shiftc
c
c dog
      vrad  = AU_day2km_s * helvrad ! au /day --> km /s
 	vnorm = AU_day2km_s * helvnorm	! au /day --> km /s
	waz = mi1a(4) * pi180	! degree to radian
c  correct modulus of velocity
	vxcor = mi1a(2) * sin(waz) + vnorm ! get velocity components 
	vycor = mi1a(2) * cos(waz) + vrad ! get velocity components
	vcor = sqrt( vxcor**2 + vycor**2 ) ! correct modulus of velocity

      select case (sc) 
      case (1)
         write(16,*) '-1.2 Degree, Pizzo correction'
	azimuth_cor = atan2(vxcor,vycor) / pi180 - 0.25646 /
	1	  heldissun ! azimuth correction depending on spacecraft distance
      case (2)
         write(16,*) ' 1.2 Degree, Pizzo correction' 
	   azimuth_cor = atan2(vxcor,vycor) / pi180 + 0.25646 /
	1	  heldissun ! different sign top Helios 1  
      end select
      write(16,*)vcor, azimuth_cor 
      write(16,*)
     $     'correction for effective solar diameter: |v|, and azimuth'
c
c============================================================================c
c     Start 3d analysis                                                      c
c============================================================================c     
      select case(i1a_3_on)
c     ------ i1a = off, i3 = off ------
c
      case (0)
c============================================================c
c     Start write out (I1A)                                  c
c============================================================c      
         write(16,*) mode,shift,perihelion_shift,hdm_ndm,i1a_3_on 
c     write(16,*) bitrate,mode,dismod,fformat,shift,perihelion_shift,
c     $        i1a_3_on,i2a_b
         write(16,*) heldissun,helcarlng,helcarlat,helcarrot
         write(16,*) eardissun,earcarlng,earcarlat,hseangle,earcarrot
         write(16,*) helvrad,helvnorm
         write(16,*) mi1a
         write(16,*) malph
         write(16,*) mi1b
         write(16,*) b, bsig
         write(16,*) '1-D i1a integrated'
         write(16,'(32(1x,e15.8))') fd1d1a_3*1e-20 !  -- km to cm --
         write(16,*) '1-D i1a velocities (km / s)'
         write(16,'(32(1x,e15.8))') volt1a(:,shift+1,sc)
         write(16,*) '1-D i1b'
         write(16,'(32(1x,e15.8))') fd1d1b*1e-20 !  -- km to cm --
         write(16,*) '1-D i1b velocities (km / s)'
         write(16,'(32(1x,e15.8))') volt1b(:,shift+1,sc)
         write(16,*)
     $   '3-D protonen+alpha distribution function (cm^-6 s^-3) and x,y,
     $   z velocity components (km /s)'
         write(16,33)'Maximum of distribution',maxloc(FD3D),maxval(fd3d) 
33	   format(a23,1x,3(1x,i2),1pe15.8)
c
c     ------ end i1a = off, i3 = off ------
c     ------ i1a = on, i3 = off  ------
c
      case (1)
c     
c     ------ velocity components ------
         do i = 1, 2                 ! sc
            do k = 1, 2              ! perishift
               do l = 1, 2           ! shift
                  do m = 1, 32       ! energy
                     do n = 1, 9     ! elevation
                        do o = 1, 16 ! azimuth
                           vx(n,o,m,l,k,i) = volt1a(m,l,i)
     $                          *cos(azimuth(o,l,k,i))
     $                          *cos(elevation(n,i))
                           vy(n,o,m,l,k,i) = volt1a(m,l,i)
     $                          *sin(azimuth(o,l,k,i))
     $                          *cos(elevation(n,i))
                           vz(n,o,m,l,k,i) = volt1a(m,l,i)
     $                          *sin(elevation(n,i))
                        end do
                     end do
                  end do
               end do
            end do
         end do
c
         select case(mode)
c     ------ NDM ------
c     
         case (0)
c     
c     ------ check for bad data ------	           
            if (minval(i1a_3(1:5,1:5,1:9)) .eq. -1) then
               i1a_3 = 0 
               go to 101
            end if
c     
c     ------ find maximum in 5 x 5 x 9 data array ------
            index = maxloc(i1a_3(1:5,1:5,1:9))
c     
c     ------ place ndm data in 9, 16, 32 array correctly ------
c     ------ elevation (always first!!) ------
            elmin = max(1,gmaxel(1)-2)
            elmax = min(9,gmaxel(1)+2)
c     ------ azimuth ------
            azmin = max(1,gmaxaz(1)-2)
            azmax = min(16,gmaxaz(1)+2)
c     ------ energy(i1a)/frequency(i3) ------
            enmin = max(1,gmaxen(1)-2)
            enmax = min(32,gmaxen(1)+6)
c     
            if ( (elmax-elmin) .ne. 4 .or. 
     $           (azmax-azmin) .ne. 4 .or.
     $           (enmax-enmin) .ne. 8 ) then
               i1a_3 = 0 
               go to 101
            endif
            
            i1a_3_all(elmin:elmax,azmin:azmax,enmin:enmax) =
     $           i1a_3(1:5,1:5,1:9)
c     ------ coordinate transformation ------
c     if (sc .eq. 2) then
            if (.False.) then
               allocate (i1a_dummy(9,16,32))
               do i = 1, 9
                  do j = 1, 16
                     i1a_dummy(i,j,:) = i1a_3_all(10-i,17-j,:)
                  end do 
               end do
               i1a_3_all = i1a_dummy(:,:,:)
               deallocate (i1a_dummy)  
            endif
c     
c     ====== end NDM ======
c     ------ HDM ------
c     
         case(1)         
c     ------ bad data ------
c     
            if (minval(i1a_3) .eq. -1) then 
c     ------ select quasi NDM array in 7 x 7 x 32 HDM data array ------
c     
               index = maxloc(i1a_3) 
c     
               elmin = max(1,index(1)-2) 
               azmin = max(1,index(2)-2)
               enmin = max(1,index(3)-2)
c     
               elmax = min(7,index(1)+2)
               azmax = min(7,index(2)+2)
               enmax = min(32,index(3)+6) 
c     ------ bad bad data ------ 
               imin = minval(i1a_3(elmin:elmax,azmin:azmax,enmin:enmax))
               if (imin .eq. -1) then
                  i1a_3 = 0
                  go to 101
               end if
c     ------ eliminate bad data fragments ------
               where (i1a_3 .eq. -1) i1a_3 = 0 
               hdm_ndm = -1
c     
c     ------ end bad data ------
            endif
c     ------ place data correctly in 9 x 16 x 32 data array ------
c     
            i1a_3_all(2:8,5:11,:) = i1a_3
c     ------ coordinate transformation ------
c     if (sc .eq. 2) then
            if (.False.) then
               allocate (i1a_dummy(9,16,32))
               do i = 1, 9
                  do j = 1, 16
                     i1a_dummy(i,j,:) = i1a_3_all(10-i,17-j,:)
                  end do 
               end do
               i1a_3_all = i1a_dummy(:,:,:)
               deallocate (i1a_dummy)  
            endif
c     
c     ====== end HDM ======
         end select               
c==============================================================c
c     Either NDM(5,5,9) or HDM(7,7,32) has been placed         c
c     correctly within the full (9,16,32) data array           c
c     Thus there is no longer any difference between NDM       c
c     and HDM mode                                             c     
c==============================================================c 
c
c     ------ instrument i1aint ------
 101     fd1d1a_3 = i1a_3_int*k1aint/volt1a(:,shift+1,sc)**2
c     
c     ------ instrument i1b ------
         fd1d1b = i1b*k1b/volt1b(:,shift+1,sc)**2
c     
c     ------ 3d distribution ------
         do i = 1, 9         ! elevation
            do j = 1, 16     ! azimuth
               do k = 1, 32  ! energy   
                  fd3d(i,j,k) = (i1a_3_all(i,j,k)*k1a(sc,i))/
     $                 (volt1a(k,shift+1,sc)**4*cos(elevation(i,sc)))
               end do
            end do
         end do
c     ------ cgs units ------
         fd3d = fd3d*1e-20  
c============================================================c
c     Start write out (I1A)                                  c
c============================================================c      
         write(16,*) mode,shift,perihelion_shift,hdm_ndm,i1a_3_on 
c     write(16,*) bitrate,mode,dismod,fformat,shift,perihelion_shift,
c     $        i1a_3_on,i2a_b
         write(16,*) heldissun,helcarlng,helcarlat,helcarrot
         write(16,*) eardissun,earcarlng,earcarlat,hseangle,earcarrot
         write(16,*) helvrad,helvnorm
         write(16,*) mi1a
         write(16,*) malph
         write(16,*) mi1b
         write(16,*) b, bsig
         write(16,*) '1-D i1a integrated'
         write(16,'(32(1x,e15.8))') fd1d1a_3*1e-20        !  -- km to cm --
         write(16,*) '1-D i1a velocities (km / s)'
         write(16,'(32(1x,e15.8))') volt1a(:,shift+1,sc)
         write(16,*) '1-D i1b'
         write(16,'(32(1x,e15.8))') fd1d1b*1e-20          !  -- km to cm --
         write(16,*) '1-D i1b velocities (km / s)'
         write(16,'(32(1x,e15.8))') volt1b(:,shift+1,sc)
         write(16,*)
     $   '3-D protonen+alpha distribution function (cm^-6 s^-3) and x,y,
     $   z velocity components (km /s)'
         write(16,33)'Maximum of distribution',maxloc(FD3D),maxval(fd3d)
         do k = 1, 32           ! energy 
            do j = 1, 9         ! elevation
               do i = 1, 16     ! azimuth
                  if (fd3d(j,i,k) .ne. 0.) then
                     write(16,1000) i,j,k,fd3d(j,i,k),i1a_3_all(j,i,k),
     $                    vx(j,i,k,shift+1,perihelion_shift+1,sc),
     $                    vy(j,i,k,shift+1,perihelion_shift+1,sc),
     $                    vz(j,i,k,shift+1,perihelion_shift+1,sc),
c  dog    
     $	                volt1a(k,shift+1,sc),
     $                    azimuth(i,shift+1,perihelion_shift+1,sc),
     $                    elevation(j,sc)
c  dog

                  endif
               end do
            end do 
         end do
c     
c     ====== end i1a ======
c     ------ i1a = off, i3 = on ------
c     
      case(2)
c     ------ velocity components ------
         do i = 1, 2                 ! sc
            do k = 1, 2              ! perishift
               do l = 1, 2           ! shift
                  do m = 1, 16       ! frequency(velocity)
                     do n = 1, 9     ! elevation
                        do o = 1, 16 ! azimuth
                           vx(n,o,m,l,k,i) = volt3(m,i)
     $                          *cos(azimuth3(o,l,k,i))
     $                          *cos(elevation3(n,i))
                           vy(n,o,m,l,k,i) = volt3(m,i)
     $                          *sin(azimuth3(o,l,k,i))
     $                          *cos(elevation3(n,i))
                           vz(n,o,m,l,k,i) = volt3(m,i)
     $                          *sin(elevation3(n,i))
                        end do
                     end do
                  end do
               end do
            end do
         end do
c     
         select case(mode)
c     ------ NDM ------
c     
         case (0)
c     
c     ------ check for bad data ------	           
            if (minval(i1a_3(1:5,1:5,1:9)) .eq. -1) then
               i1a_3 = 0 
               go to 201
            end if
c     
c     ------ find maximum in 5 x 5 x 9 data array ------
            index = maxloc(i1a_3(1:5,1:5,1:9))
c     
c     ------ place ndm data in 9, 16, 32 array correctly ------
c     ------ elevation (always first!!) ------
            elmin = max(1,gmaxel(1)-2)
            elmax = min(9,gmaxel(1)+2)
c     ------ azimuth ------
            azmin = max(1,gmaxaz(1)-2)
            azmax = min(16,gmaxaz(1)+2)
c     ------ energy(i1a)/frequency(i3) ------
            enmin = max(1,gmaxen(1)-2)
            enmax = min(32,gmaxen(1)+6)
c     
            if ( (elmax-elmin) .ne. 4 .or. 
     $           (azmax-azmin) .ne. 4 .or.
     $           (enmax-enmin) .ne. 8 ) then
               i1a_3 = 0 
               go to 201
            endif
            
            i1a_3_all(elmin:elmax,azmin:azmax,enmin:enmax) =
     $           i1a_3(1:5,1:5,1:9)
c     ------ coordinate transformation ------
c     if (sc .eq. 2) then
            if (.False.) then
               allocate (i1a_dummy(9,16,32))
               do i = 1, 9
                  do j = 1, 16
                     i1a_dummy(i,j,:) = i1a_3_all(10-i,17-j,:)
                  end do 
               end do
               i1a_3_all = i1a_dummy(:,:,:)
               deallocate (i1a_dummy) 
            endif
c     
c     ====== end NDM ======
c     ------ HDM ------
c     
c     ----- n.b.: the part i1a_3(7 x 7 x (17-32)) contains the mass spectrum ------
c     
         case(1)              
c     ------ bad data ------
c     
            if (minval(i1a_3(:,:,1:16)) .eq. -1) then     
c     ------ find maximum in 7 x 7 x 16 data array ------
c     
               index = maxloc(i1a_3(:,:,1:16)) 
c     
               elmin = max(1,index(1)-2) 
               elmax = min(7,index(1)+2)
c     
               azmin = max(1,index(2)-2)
               azmax = min(7,index(2)+2)
c     
               enmin = max(1,index(3)-2)                  
               enmax = min(16,index(3)+6)         
c     ------ bad bad data ------    
               imin = minval(i1a_3(elmin:elmax,azmin:azmax,enmin:enmax))
               if (imin .eq. -1) then
                  i1a_3 = 0
                  go to 201
               end if
c     ------ eliminate bad data fragments ------
               where (i1a_3 .eq. -1) i1a_3 = 0 
               hdm_ndm = -1
            endif            
c     ------ place data correctly in 9 x 16 x 16 data array ------
c     
            i1a_3_all(2:8,5:11,1:16) = i1a_3(:,:,1:16) 
c     ------ coordinate transformation ------
c     if (sc .eq. 2) then
            if (.False.) then
               allocate (i1a_dummy(9,16,32))
               do i = 1, 9
                  do j = 1, 16
                     i1a_dummy(i,j,:) = i1a_3_all(10-i,17-j,:)
                  end do 
               end do
               i1a_3_all = i1a_dummy(:,:,:)
               deallocate (i1a_dummy)  
            endif
c     
c     ====== end HDM ======      
         end select
c==============================================================c
c     Either the NDM(5,5,9) or HDM(7,7,32) has been placed     c
c     correctly within the full (16,9,32) data array           c
c     Thus there is no longer any differences between NDM      c
c     and HDM mode                                             c     
c==============================================================c
c
c     ------ instrument i13int ------
 201     fd1d1a_3(1:16) = i1a_3_int(1:16)*ki3int/volt3(:,sc)**2
c     
c     ------ instrument i1b ------
         fd1d1b = i1b*k1b/volt1b(:,shift+1,sc)**2
c     
c     ------ 3d distribution ------
         do i = 1, 9            ! elevation
            do j = 1, 16        ! azimuth
               do k = 1, 16     ! frequency  
                  fd3d(i,j,k) = (i1a_3_all(i,j,k)*ki3(sc,i))/
     $                 (volt3(k,sc)**4*cos(elevation3(i,sc)))
               end do
            end do
         end do
c     ------ cgs units ------
         fd3d = fd3d*1e-20  
c============================================================c
c     Start write out (I3)                                   c
c============================================================c      
         write(16,*) mode,shift,perihelion_shift,hdm_ndm,i1a_3_on 
c     write(16,*) bitrate,mode,dismod,fformat,shift,perihelion_shift,
c     $        i1a_3_on,i2a_b
         write(16,*) heldissun,helcarlng, helcarlat, helcarrot
         write(16,*) eardissun, earcarlng, earcarlat, hseangle,
     $        earcarrot
         write(16,*) helvrad, helvnorm
         write(16,*) mi1a
         write(16,*) malph
         write(16,*) mi1b
         write(16,*) b, bsig
         write(16,*) '1-D i1a integrated'
         write(16,'(32(1x,e15.8))') fd1d1a_3*1e-20  !km to cm
         write(16,*) '1-D i1a velocities (km / s)'
         write(16,'(16(1x,e15.8),16(1x,e15.8))') volt3(:,sc), 
     $        fd1d1a_3(17:32)
         write(16,*) '1-D i1b'
         write(16,'(32(1x,e15.8))') fd1d1b*1e-20    !km to cm
         write(16,*) '1-D i1b velocities (km / s)'
         write(16,'(32(1x,e15.8))') volt1b(:,shift+1,sc)
         write(16,*)
     $   '3-D protonen+alpha distribution function (cm^-6 s^-3) and x,y,
     $   z velocity components (km /s)'
         write(16,33)'Maximum of distribution',maxloc(FD3D),maxval(fd3d)
         do k = 1, 16         ! frequency
            do j = 1, 9       ! elevation
               do i = 1, 16   ! azimuth
                  if (fd3d(j,i,k) .ne. 0.) then
                     write(16,1000) i,j,k,fd3d(j,i,k),i1a_3_all(j,i,k),
     $                    vx(j,i,k,shift+1,perihelion_shift+1,sc),
     $                    vy(j,i,k,shift+1,perihelion_shift+1,sc),
     $                    vz(j,i,k,shift+1,perihelion_shift+1,sc),
c  dog    
     $	                volt1a(k,shift+1,sc),
     $                    azimuth(i,shift+1,perihelion_shift+1,sc),
     $                    elevation(j,sc)
c  dog

                  endif
               end do
            end do 
         end do
c     
c     ====== end i3 =======
      end select
c============================================================c
c     Handle electron data correctly to calculate the        c
c     distribution function (did not change this)            c
c============================================================c
 7    if (minval(i2b) .lt. 0) eminus = -1 
      indata = reshape(i2b,(/128/)) 
      compressed = 255         
      do i = 1, 128             
         do j = 1, dim          
            if (norm(j) .eq. indata(i)) then  
               compressed(i) = j
               exit
            end if
         end do
      end do
c
c     Spacecraft dependend decompression
      select case (sc)
      case (1)
         electron =  h1corr(compressed)
      case (2)
         electron = h2corr(compressed)
      end select
c    
c     Correction for electrostatic spacecraft potential, velocity
c dog      phi_sc = 3.85 * (2.13 - dlog(mi1a(1) * heldissun**2)) 
c dog      v_sw2  = v_sc2 - si_me2 * phi_sc 
c dog      v      = sqrt(v_sw2)	
c dog      v4     = v_sw2**2          
c
c dog
      if (mi1a(1).gt.1.e-10)then 
	  phi_sc = 3.85 * (2.13 - dlog(mi1a(1) * heldissun**2)) !correct for
	  write(18,*)directory//'/'//
     2          file(1:2)//'y'//file(3:8)//out//'_'//endung//'.'//shiftc
	 else
	  phi_sc = 10
	endif
	
      if (mi1a(1).lt.1.e-10)then
c	write(6,*)'mila(1).lt.0., density< 1.e-10'
c       read(5,*)
	write(17,*)'mila(1).lt.0., density< 1.e-10   ' ,directory//'/'//
     2          file(1:2)//'y'//file(3:8)//out//'_'//endung//'.'//shiftc
	
 
	go to 3
	 end if

	  

	 v_sw2  = v_sc2 - si_me2 * phi_sc ! spacecraft potential see programs
				                      ! by K. Ivory
	 sq=.true.
	 do n=1,6
	   if (v_sw2(n).le.0)then
	   sq=.false.
	   endif
	 enddo
	if (sq.eqv.(.false.)) then
	  write(17,*)'v<0 ' ,directory//'/'//
     2          file(1:2)//'y'//file(3:8)//out//'_'//endung//'.'//shiftc
	  goto 3
	endif
	 v      = sqrt(v_sw2)	          ! return argument: velocity (km/s)
	 v4     = v_sw2**2	                ! 4th power 
c---end correction for electrostatic spacecraft potential, velocity
c dog
      electron_c = reshape(electron(1:128),shape) 
c     
c     2D distribution function	 
      do i = 1, 16		
         electron_distribution(:,i) = electron_c(:,i) / v4(i) * cr(sc)
         vx_e(:,i) = v(i)*cos(azimuthe(sc,shift+1,perihelion_shift+1,:))
         vy_e(:,i) = v(i)*sin(azimuthe(sc,shift+1,perihelion_shift+1,:))
      end do
c     in cgs!
      electron_distribution = electron_distribution / 1.e12 
c============================================================c
c     Continue write out                                     c
c============================================================c      
      write(16,*)
     $  '2-D electron distribution function (cm^-6 s^-3) and & radial ve
     $  locity (km/s)' 
      write(16,*)'Maximum of distribution',
     $     maxloc(electron_distribution), maxval(electron_distribution)
      select case (eminus)
      case (-1)                
         write(16,*) 'no electron data available'
      case (0)   
	 do k = 1, 8 
            do i = 1, 16
               if (electron_distribution(k,i) .gt. 0.d0) then
                  write(16,1001) k,i,electron_distribution(k,i),
     $                 electron_c(k,i), vx_e(k,i)/1.e3, vy_e(k,i)/1.e3 !1.e3 m to km
               endif
            end do 
         end do
      end select
c
      go to 3 
c      
      stop
      end
c      
      function out_0(i)
c     
c     writes the correct 0 into the output file. If minutes or seconds
c     have only a single digit, i.e, 1, then the routine produces a
c     character variable '01', otherwise1 i.e., 11, the character
c     variable '11' remains. Numbers larger than 99 as well as
c     negative numbers are not allowed.
c     
      implicit none
      integer i
      character out_0 * 2
      if (i .lt. 0) go to 99
      if (i .gt. 99) go to 99
      if (i .lt. 10) then
         out_0(1:1) = '0'
         write(out_0(2:2),'(i1)') i
      else
         write(out_0(1:2),'(i2)') i
      end if
      return
 99   write(6,*) 'Wrong number in out_0'
      stop
      end
      

