c**************************************************************** C email from lindau c corrected for work with power station on pc 2002,3,30 c*************************************************************** program helios_3d_table2 USE PORTLIB include 'type_helios_3d_table.inc' ! type declaration real seconds c c c azimuth(data,shift,perishift,spacecraft) c elevation(data,spacecraft) c volt1a(data,shift,spacecraft) c volt1b(data,shift,spacecraft) c k1a(spacecraft,:) c vx(azimuth,elevation,energy,shift,perishift,spacecraft) c c CEM calibration data see blue book page 181 (BMFT-FB-W 81-015) c c azimuthe(spacecraft,shift,perishift,:) c include 'kalibrations_daten.inc' ! data file with azimuth, elevation, ! and energy data go to 1 ! c output file parameter 2 out(1:1) = 'h' ; out(4:4) = 'm' ; out(7:7) = 's' !define letters for !output file s = sec, m= min, h = hour c default values for input file and directory file = 'h d .txt' !reads raw data files, start in directory where ! these files are located. The raw data files ! are provided by the MPAe compressed on 3 CD's directory = 'Y D ' c==========input parameter: spacecraft id, year day of year ============ c c more or less selfexplaining c write(6,*) !get spacecraft id 4 write(6,*) !get spacecraft id 1 'spacecraft ? (= 1,2 other number will cancel the programn)' read(5,*) sc c=================== get system time (not necessary feature) ========== c call second(seconds) ; write(6,*) 'Runtime ', seconds/60., seconds= secnds(seconds) ; write(6,*) 'Runtime ', seconds/60., 1 ' min' c==================== if wrong spacecraft id stop ===================== c c may also be used to halt the program after one day was processed, c otherwise enter new input parameter (if spacecraft id changes, c this must be done in the previous statement.) c if ((sc .lt. 1) .or. (sc .gt. 2)) stop 11 write(6,*) 'year: yy?' ; read(5,'(a2)') year ! get year an day of year read(year,'(i2)') idummy ! check if correct year 1974 to 1986 if ((idummy .lt. 74) .or. (idummy .gt. 86)) go to 11 ! get correct value 12 write(6,*) 'day of the year: ddd (with 0: e.g., 002 for day 2)' read(5,'(a3)') doy read(doy,'(i3)') idummy ! check day of year if ((idummy .lt. 1) .or. (idummy .gt. 366)) go to 12! get correct value write(file(2:2),'(i1)') sc ! determine fielname with spacecraft id, file(3:4) = year ; file(6:8) = doy ! year and day of year c c====================make directory==================================== c c At this stage the directories are named with a capital letter Y c for the year followed by a 2 digit value (e.g., 75 for 1975), c next comes a capital D and a 3 digit value for the day of the c year. The directory is made from the input parameters. c directory(2:3) = year ; directory(5:7) = doy ! make new directory c call system('mkdir '//directory,status) status = SYSTEM('mkdir '//directory ) write(*,*)status if (status .ne. 0) then ! use old directory ? write(6,*) 1 '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 write(6,*)'***',file read(5,*) open (15, file=file) c c====================end make directory=============================== c c@@@@@@@@@@@@@@@@@@@@@@@@@@@ get spectra for one day @@@@@@@@@@@@@@@@@ c c ==================== begin read raw data ======================= c 3 read(15,*,end=4) yr_doy, msec !yy.ddd. millisec c c end = 4 ---> if end of file get new filename c read(15,*) bitrate, mode, fformat, dismod, shift, ! mode = hdm, ndm 1 perihelion_shift,i1a_3_on, i2a_b !fformat, dismod not used !the other parameter see file output_file_documentation.txt read(15,*) c see outpout_file_documentation.txt for the different parameter read(15,*) hellngasc, heldissun, helvrad, helvnorm, helcarlng, 1 helcarlat, helcarrot ! read(15,10) eardissun, earcarlng,earcarlat, hseangle,earcarrot ! 10 format(e15.8,e15.8,x,e15.8,e15.8,x,i6) read(15,*) read(15,*) read(15,*) read(15,*) gmaxel, gmaxaz, gmaxen, gmass ! gives the maximum for ! elvation, azimuth, and enegery and mass channel (I3 only) read(15,*) mi1a ! one fluid proton parameter I1a: density, velocity, ! temperature, azimuth-, and elevation of flow vector read(15,*) malph!one fluid alpha parameter I1a: density, velocity, ! temperature read(15,*) mi1b!one fluid proton parameter I1b: density, velocity, ! temperatur read(15,*) b, bsig !magnetic field vector bx,by,bz and deviation bsig read(15,*) i1b ! elcetron data read(15,*) i1aint ! integrated i1a read(15,*) i2b ! integrated i2b read(15,*) i1a ! count rate i1a i1a_all = 0 ; minus = 0 ; eminus = 0 ! clear data c ====================== end read procedure ============================== c c====================== start filename procedure ========================= c write hour, minutes, seconds in the filename and open th = int(msec * millisec2hour) ! get hour, tm = int((msec * millisec2hour - th)*60) ! minute, ts = int(((msec * millisec2hour - th)*60 - tm)*60) !second out(2:3) = out_0(th) ; out(5:6) = out_0(tm) ; out(8:9) = 1 out_0(ts) ! define output file parameter write(shiftc,'(i1)') shift ! include shift in file name: no shift = 0, ! shift = 1 select case (mode) ! include data mode in filename: ndm or hdm !mode = 0 = ndm; mode = 1 = hdm case (0) endung = 'ndm' case (1) endung = 'hdm' end select open(16, file = directory//'/'//file(1:2)//'y'//file(3:8)/ !open out- 1 /out//'_'//endung//'.'//shiftc) ! put file c as example: h1y75d002h23m59s49_hdm.0 for Helios 1 (h1), year c 1975, day 002, 23 hours, 59 minutes, 49 seconds, hdm mode no c shift; see outpout_file_documentation.txt write(6,*) directory//'/'//file(1:2)//'y'//file(3:8)//out/ 2 /'_'//endung//'.'//shiftc c====================================================================== c c Correction of the azimuth angle according to Pizzo et al 1983 and for c the effective sun diameter (azimuth and modulus of velocity) c c====================================================================== c correction for the effectiv sun diameter 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' c 0.25646 is the diameter of the sun from 1 AU in degree. A better c value is: (6.96 * 10^8) / (1.4959787 * 10^11) * 180 / pi = 0 c .26665670. The numbers are taken from the Explanatory c supplement to the Astronomical Almanac, 1992 azimuth_cor = atan2(vxcor,vycor) / pi180 - 0.25646 / 1 heldissun! azimuth correction depending on spacecraft distance c========================write correction to file ============= 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,*) 1 'correction for effective solar diameter: |v|, and azimuth' c==========================end correction============================== c c===========================start 3d analysis ========================= c c%%%%%%%%%%%%%%%%%%%%%% select mode %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% select case(mode) ! ndm or hdm mode c!!!!!!!!!!!!!!!!!!!!!!!!!!!! ndm mode !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! case (0) ! ndm data mode c--------------------------bad data----------------------------------- if (minval(i1a) .eq. -1) then i1a = 0 ; go to 7! only bad data, go to electrons end if c------------------------end bad data -------------------------------- c c-----------------------position maximum----------------------------- c index = maxloc(i1a) ! find maximum in 5,5,9 matrix c postion maximum in 16,9,32 matrix; avoid array boundary c violations, i.e. i1min = max(1,index(1)-2) --> if index(1)-2 c less than 1 ==> i1min = 1, else i1min = index(1)-2, etc... i1min = max(1,index(1)-2) ; i2min = max(1,index(2)-2) i3min = max(1,index(3)-2) ; i11min = max(1,gmaxaz(1)-2) i21min = max(1,gmaxel(1)-2) ; i31min = max(1,gmaxen(1)-2) i1max = min(7,index(1)+2) ; i2max = min(7,index(2)+2) i3max = min(32,index(3)+6) ; i11max = min(16,gmaxaz(1)+2) i21max = min(9,gmaxel(1)+2) ; i31max = min(32,gmaxen(1)+6) c c postion maximum in 16,9,32 matrix correctly c c-----------------------end position maximu-------------------------- c c----------------------set data in big matrix 16,9,32--------------- i1a_all(i11min:i11max,i21min:i21max,i31min:i31max) = 1 i1a(i1min:i1max,i2min:i2max,i3min:i3max) c------------------------end set data------------------------------ c c------------------------coordinate transformation----------------- allocate (i1a_dummy(16,9,32)) select case (sc) ! which spacecraft case (1) ! Helios 1 no coordinate transformation continue case (2) ! Helios 2 revers azimuth and elevation do i = 1, 16 ; do j = 1, 9 i1a_dummy(i,j,:) = i1a_all(17-i,10-j,:) end do ; end do i1a_all = i1a_dummy(:,:,:) end select deallocate (i1a_dummy) c-----------------------end coordinate transformation--------------- write(6,*)'***-2' read(5,*) c c!!!!!!!!!!!!!!!!!!!!end ndm mode!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c c!!!!!!!!!!!!!!!!!!!!!!!!!!!hdm mode!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! case(1) ! hdm data mode c-------------------------------------------------------------------- c if data frames are not correct (i.e. i1a(i,j,k)<0) ; check if c ndm mode applies, i.e. if around maximum a 5,5,9 matrix exists c without -1 c-------------------------------------------------------------------- if (minval(i1a) .eq. -1) minus = -1 c c it is possible that in the hdm mode a frame of data was not c transmitted. Then the data were denoted by -1. Check if it is c possible to handle the hdm mode in as a ndm format, and if so, c continue, otherwise get new data c c---------------------bad data--------------------------------------- select case (minus) ! bad data case (0) ! all data frames are correct continue case (-1) ! a data frame was not transmitted c c-----------------------position maximum----------------------------- c index = maxloc(i1a) ! find maximum, see ndm mode i1min = max(1,index(1)-2) ; i2min = max(1,index(2)-2) i3min = max(1,index(3)-2);i11min = max(1,gmaxaz(1)-2) i21min = max(1,gmaxel(1)-2) ; i31min = max(1 1 ,gmaxen(1)-2) i1max = min(7,index(1)+2) ; i2max = min(7,index(2)+2) i3max = min(32,index(3)+6) ; i11max = min(16 1 ,gmaxaz(1)+2) i21max = min(9,gmaxel(1)+2) ; i31max = min(32 6 ,gmaxen(1)+6) c--------------------------------end position maximum----------------- c c----------------------------bad bad data------------------------------ c even data in ndm format are bad, if (minval(i1a(i1min:i1max,i2min:i2max,i3min:i3max)) 1 .eq. -1) then i1a = 0 ; go to 7 !kill all data, process the electron data end if c----------------------------end bad bad data------------------------ c c if data in ndm format are okay, set bad data to 0 an continue in c hdm mode c c--------------------------good bad data----------------------------- where (i1a .eq. -1) i1a = 0 ! set -1 to 0 and continue end select ! minus (bad data) allocate (i1a_dummy(7,7,32)) c------------------------coordinate transformation----------------- select case (sc) ! which spacecraft case (1) ! Helios 1 write to i1a_all i1a_all(6:12,2:8,:) = i1a case (2) ! Helios 2 revers azimuth and elevation do i = 1, 7 ; do j = 1, 7 i1a_dummy(i,j,:) = i1a(8-i,8-j,:) end do ; end do i1a_all(6:12,2:8,:) = i1a_dummy end select ! sc c------------------------end coordinate transformation-------------- deallocate (i1a_dummy) end select ! mode c!!!!!!!!!!!!!!!!!!!!!!!!!!!end hdm mode!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c%%%%%%%%%%%%%%%%%%%%%% end select mode %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c c No longer any differences between ndm and hdm mode c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! fd3d = 0.d0 c______________________________________________________________________ c Attention: Kolesnik uses 4th power of velocity in BMFT-Fb-W-82 c -002 the second power is used. Here correct to second power. c_____________________________________________________________________ c---------------------------1d and 3d distribution functions---------- c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ fd1d1a = i1aint * k1aint / volt1a(:,shift+1,sc)**2 !calibrate c fd1d1b = i1b * k1b / volt1b(:,shift+1,sc)**2 ! integrated instruments c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% do i = 1, 16 ; do j = 1, 9 ; do k = 1, 32 !estimate 3-D distribution fd3d(i,j,k) = (i1a_all(i,j,k) * k1a(sc,j)) / (volt1a(k,shift!from 1 +1,sc)**4 * cos(elevation(j,sc)))!count rates end do ; end do ; end do fd3d = fd3d * 1e-20 ! in cgs c------------------------end 1d and 3d distribution functions---------- c======================end 3d analysis============================ c===========================electrons ============================= 7 if (minval(i2b) .lt. 0) eminus = -1 ! check if good data: i2b>=0 indata = reshape(i2b,(/128/)) ! shuffle 8,16 matrix into 128 array compressed = 255 ! define zeros, see conv_e1.pro and distfunc.pro do i = 1, 128 ! for details. Program was written by K. Ivory do j = 1, dim ! encode ascii data to compressed format if (norm(j) .eq. indata(i)) then compressed(i) = j exit end if end do end do c-----------------------------spacecraft dependend decompression------------ select case (sc) case (1) electron = h1corr(compressed) case (2) electron = h2corr(compressed) end select c--------------------------end spacecraft dependend decompression-------- c c--------correction for electrostatic spacecraft potential, velocity----- phi_sc = 3.85 * (2.13 - alog(mi1a(1) * heldissun**2)) !correct for if (mi1a(1).lt.1.e-10)then write(6,*)'mila(1).lt.0., density< 1.e-10' read(5,*) go to 3 end if v_sw2 = v_sc2 - si_me2 * phi_sc ! spacecraft potential see programs ! by K. Ivory v = sqrt(v_sw2) ! return argument: velocity (km/s) v4 = v_sw2**2 ! 4th power c---end correction for electrostatic spacecraft potential, velocity----- electron_c = reshape(electron(1:128),shape) !reshape uncompressed data ! from 128 array to 8,16 matrix c c!!!!!!!!!!!!!!!!!!!!!!!! 2d distribution function !!!!!!!!!!!!!!!!!!!!!!!!! c do i = 1, 16 ! calculate 2 D distribution function and l ! eclipticavelocity components if(v4(i).eq.0.)then write(6,*)'v4(i)=',v4(i) read(5,*) go to 3 endif 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 electron_distribution = electron_distribution / 1.e12 !in cgs c!!!!!!!!!!!!!!!!!!!! end 2d distribution function !!!!!!!!!!!!!!!!!!!!!!!!! c============================end electron=================================== c===============================start write out============================= c---------------one fluid parameter, magnetic field, spacecraft parameter-- write(16,*) mode,shift,perihelion_shift, minus, i1a_3_on write(16,*) heldissun,helcarlng, helcarlat, helcarrot write(16,*) eardissun, earcarlng, earcarlat, hseangle, 1 earcarrot write(16,*) helvrad, helvnorm write(16,*) mi1a write(16,*) malph write(16,*) mi1b write(16,*) b, bsig c------------end one fluid parameter, magnetic field, spacecraft parameter-- c c-----------------------------integrated-------------------------------- c c integrated i1a and i1b + velocities c write(16,*) '1-D i1a integrated' write(16,'(32(1x,e15.8))') fd1d1a * 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) c--------------------------end integrated------------------------------- c c-----------------------d data---------------------------------------- c c only values > 0 are printed, c write(16,*) 1'3-D protonen+alpha distribution function (cm^-6 s^-3) and x,y, 2z velocity components (km /s)' write(16,*) 'Maximum of distribution', maxloc(FD3D), maxval(fd3d) do k = 1, 32 ; do j = 1, 9 ; do i = 1, 16 if (fd3d(i,j,k) .ne. 0.) write(16 1 ,'(1x,i2,1x,i1,1x,i2,1x,1pe15.8,1x,i5,3(1x,1pe15.8))') i 1 ,j,k,fd3d(i,j,k),i1a_all(i,j,k), vx(i,j,k,shift+1 2 ,perihelion_shift+1,sc),vy(i,j,k,shift+1 2 ,perihelion_shift+1,sc),vz(i,j,k,shift+1 3 ,perihelion_shift+1,sc) end do ; end do ; end do c--------------------------end 3d data---------------------------------- c------------------------2d electron data------------------------------- write(16,*) 1 '2-D electron distribution function (cm^-6 s^-3) and & radial ve 2locity (km/s)' write(16,*) 'Maximum of distribution' 3 ,maxloc(electron_distribution), maxval(electron_distribution) select case (eminus) case (-1) ! bad data write(16,*) 'no electron data available' case (0) do k = 1, 8 ; do i = 1, 16 if (electron_distribution(k,i) .gt. 0.d0) write(16 1 ,'(1x,i2,1x,i2,1x,1pe15.8,1x,i5,2(1x,1pe15.8))') k,i 1 ,electron_distribution(k,i),electron_c(k,i) ,vx_e(k 2 ,i)/1.e3,vy_e(k,i)/1.e3 end do ; end do c /1.e3 m to km end select c get next data set at same day go to 3 c c@@@@@@@@@@@@@@@@@@@@@@@@@@@ end get spectra for one day @@@@@@@@@@@@@@@@@ c c========================initialize program ============================= 1 azimuth = azimuth * pi180 ! degree to radian elevation = elevation * pi180 ! degree to radian azimuthe = azimuthe * pi180 ! degree to radian c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c c calculate all possible velocitiy components c c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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(o,n,m,l,k,i) = volt1a(m,l,i) 1 * cos(azimuth(o,l,k,i)) * cos(elevation(n,i)) vy(o,n,m,l,k,i) = volt1a(m,l,i) 1 * sin(azimuth(o,l,k,i)) * cos(elevation(n,i)) vz(o,n,m,l,k,i) = volt1a(m,l,i) 1 * sin(elevation(n,i)) end do ; end do ; end do ; end do ; end do ; end do c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c get data c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! go to 2 stop end 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