C Electron_Moment_Avg.for C H. K. Hills, Hughes STX, April 1995. C Entered as B44511 in NSSDC's TRF (Technical Reference File). C This program was written and checked on a Micro-VAX. c Oct. 6, 1994, program started, ignoring angles but averaging other parameters. c March, 1995 additions of average angle calculations. c Mar. 24. Corrected the formats, made last angle conform to 0 - 360 range c instead of +- 180, to match input range. c April 13, 1995 Converted sigma from radians to degrees for computed angles. c April 21, 1995 Converted some calculations to double precision (for sigma). c April 26, 1995 Added error return to data write statement. C This program was written to operate on ISEE 3 data set 78-079A-01O, which is C a data set of solar wind electron moments from the Los Alamos experiment. C The original data set has nominal 168-s resolution; this program generates C hourly averages of the data. The original data set included eight data C parameters: electron density; bulk flow speed; flow azimuth; Tperp; Tpar; C the azimuth of the max Temp; the heat flux; and the azimuth of the heat flux. C Because of the possible distortions induced by simply averaging the azimuth C angles, each angle is averaged as follows: the plane components of a unit C vector are computed for each given angle and these components are averaged C over the hour. Then an average angle is computed from the average components. C Sigma is computed and output for each parameter and for the calculated average C angles. The sum of the squares of the averages of the two components is also C given, for each of the three averaged angles. The deviation of this value C from 1.0 is also a measure of the scatter in the original input points. C If someone wants detailed or highly accurate angular information, they C should be looking at the higher time resolution data. C The program checks and requires that time increase monotonically. If not, it C writes an error message and stops. C Procedure: C1. Zero each parameter counter and sum. C2. Read 1st record, establish the year, month, day, and hour. C3. Increment each parameter counter and sum (iff each passes reasonable check). C4. For subsequent records, require all these time values to be the same. C5. When they are not, compute each average, store in output file, zero each C parameter counter and sum, remember the new time, go to 3. Program Electron_Moment_Avg real*4 phiavg(3),sigma(14) real*4 sigma_phi(3) character*1 yesno,iecho Integer*4 itime,ltime,parmcount(14) Real*4 utime,time Integer*2 hour,uhour,express(8) character*60 dataname,infilename integer*4 iyymmdd real*8 xcomp,ycomp,sumsquare(14),sq_avg(14),plasma(8) real*8 parmsum(14),parmavg(14),rsq(3),difference,fill c equivalence(utime,uhour) data express/1,1,0,1,1,0,1,0/ ! used to skip/select angle cases fill = 9.9999 ! *** dummy value for testing *** angle=180./3.14159265 ! mult by this to convert to degrees icnt=0 ! input records echoed icnt7=0 ! output records written icthis=0 ! count input records on this current file only assign 54 to iswitch assign 501 to iprint write(6,600) 600 format(/,' For use with ISEE 3 Solar Wind Electron Moments data', 1 ' set (78-079A-01O)',/,' from the Los Alamos Experiment.', 1 ' H. K. Hills, February, 1995.',/, 2 ' FORTRAN source for this program is entered as',/, 3 ' B44511 in NSSDC''s Technical Reference File (TRF).',//, 4 ' This makes hourly averages from 168-sec data set.',/, 5 ' Sigma values, and number of samples in average are included.', 6 /, 6 ' Angles are computed by averaging their cartesian components,', 6 /,' then finding the angle for the resulting average.', 7 //, ' User enters output file name, then input file name.',/, 8 ' Only one output file is made, but can have many input files.', 8 /,' At EOF on input, get chance to enter another input file.',/, 9 ' If there are no more files to process, enter ^Z to stop.',/, 9 ' 365 full days yield 8760 records @85 Bytes = 744.6 KBytes.',//) write(6,*) 'Enter name of output file (up to 60 characters).' read(5,10,end=30) dataname open(unit=7,file=dataname,access='sequential', 1 status='new',form='formatted',iostat=ios,blocksize=234, 2 recordsize=234) if(ios.ne.0) go to 3007 write(7,2099) 2099 format(' Date Hr Density Speed Az_Sp', 1 ' Tperp Tpar Az_T Heatflux Az_HF', 2 ' Sig_Den Sig_Spd Sig_Az', 3 ' Sig_Tperp Sig_Tpar Sig_Az Sig_Heat Sig_Az', 4 ' N1 N2 N3 N4 N5 N6 N7 N8 Mag1 Mag2 Mag3') write(6,*) 'Enter number of records to process. Use -1 for all.' read(5,*) iprocstop if(iprocstop.eq.-1) assign 56 to iswitch write(6,*) 'Enter Y if you want to echo input on screen.' read(5,499) iecho 499 format(A1) if(iecho.eq.'Y'.or.iecho.eq.'y') assign 500 to iprint write(6,*) 'Enter name of first input file (up to 60 characters).' 111 read(5,10,end=30) infilename 10 format(A60) open(unit=10,file=infilename,status='old',access='sequential', 1 iostat=ios,readonly) if(ios.ne.0) go to 3010 write(6,*) 'Opened File ',infilename iend=0 icthis=0 go to 301 300 continue write(6,*)icnt,' input records; ',icnt7,' output records written.' C 365 full days yield 8760 records. 3001 write(6,*) 'Enter name of next input file (up to 60 characters).' go to 111 C1. Zero each parameter counter and sum. 301 continue do 21 k=1,14 parmcount(k)=0 parmsum(k)=0. sumsquare(k)=0. parmavg(k)=-9.999 sq_avg(k)=-9.999 21 continue do 22 k=1,3 phiavg(k)=-9.999 22 continue C2. Read 1st record, establish the year, month, day, and hour. read(10,1000,err=1330,end=330) iyymmdd,sec,utime,plasma itime=iyymmdd time=utime uhour=utime/10000. hour=uhour go to 304 ! bypass time checking (time was just set = utime) 2 read(10,1000,err=1331,end=331) iyymmdd,sec,utime,plasma 1000 format(1x,I6,1x,F8.1,1x,F8.1,5(1x,1PE11.4),/, 1 3(1PE11.4,1x)) uhour=utime/10000. C4. For subsequent records, require all these time values to be the same. 4 continue if(iyymmdd.lt.itime) go to 333 ! ERR FLAG if day backs up if(iyymmdd.eq.itime.and.utime.lt.time) go to 333 !ERR FLAG; time backs up if(iyymmdd.eq.itime.and.uhour.eq.hour) then 304 continue go to iprint ! either 500 or 501 c 500 write(6,2000) iyymmdd,utime,plasma 2000 format(1x,I6,F8.0,8(1PE8.1)) c write(8,2001) iyymmdd,sec,utime,plasma 2001 format(1x,I6,1x,F8.1,1x,F8.0,8(1PE11.4)) 501 continue c icnt=icnt+1 icthis=icthis+1 go to iswitch ! skip check of record number if no limit requested 54 if(icthis.ge.iprocstop) go to 331 56 continue C plasma(x) is an azimuth angle for x = 3, 6, 8 do 110 k=1,8 if(plasma(k).eq.fill.or.plasma(k).eq.-fill)go to 110 !skip if fill value parmcount(k)=parmcount(k)+1 parmsum(k)=parmsum(k)+plasma(k) sumsquare(k)=sumsquare(k)+(plasma(k)**2) 110 continue C Now here handle the angle components kazi = -1 do 210 k=3,8 if(express(k).eq.1) go to 210 ! skip except 3,6,8 kazi=kazi + 2 ! (9,10), (11,12), (13,14) for (x,y)1, (x,y)2, (x,y)3 if(plasma(k).eq.fill.or.plasma(k).eq.-fill)go to 210 !skip if fill value xcomp=Dcosd(plasma(k)) ycomp=Dsind(plasma(k)) parmsum(8+kazi)=parmsum(8+kazi)+xcomp parmsum(9+kazi)=parmsum(9+kazi)+ycomp parmcount(8+kazi)=parmcount(8+kazi)+1 parmcount(9+kazi)=parmcount(9+kazi)+1 sumsquare(8+kazi)=sumsquare(8+kazi)+xcomp*xcomp sumsquare(9+kazi)=sumsquare(9+kazi)+ycomp*ycomp 210 continue time=utime go to 2 endif C Now time has changed (either iyymmdd or uhour or both) if(iyymmdd.eq.itime.and.utime.lt.time) go to 333 ! ERR; time has backed up C If got to here, then either the day or the hour (or both) have increased. C5. When day or hour advances, compute each average, store in output file, zero C each parameter counter and sum, remember the new time, then go to location C just after read and continue. 350 continue C Compute each average, including angles and (x,y) components of angles do 410 jk=1,14 if(parmcount(jk).eq.0) go to 410 parmavg(jk)=parmsum(jk)/parmcount(jk) ! = c also compute here the average of components_squared sq_avg(jk)=sumsquare(jk)/parmcount(jk) ! = 410 continue Compute phi average from component averages c Indices (9,10), (11,12), (13,14) for (x,y)1, (x,y)2, (x,y)3 kav=0 do 510 jk=9,13,2 if(parmcount(jk).eq.0) go to 510 y=parmavg(jk+1) x=parmavg(jk) kav=kav+1 phiavg(kav)=atan2d(y,x) ! arctan2d(y,x); kav = 1,2,3 rsq(kav)=parmavg(jk)*parmavg(jk) + parmavg(jk+1)*parmavg(jk+1) C This is DP of x*x + y*y radius-squared. Would be 1 if no variation in values. 510 continue if(phiavg(3).lt.0.) phiavg(3) = phiavg(3) + 360. ! makes range like input Do 521 kk=1,14 if(express(k).eq.0) go to 521 ! skip indices 3,6,8 c write(6,*) kk, sq_avg(kk), parmavg(kk) difference=sq_avg(kk)-(parmavg(kk)*parmavg(kk)) if(difference.lt.0.) then write(6,*)'SQRT(negative value) ',iyymmdd,hour,kk, sq_avg(kk), 1 parmavg(kk), parmcount(kk), parmsum(kk) stop endif sigma(kk)=Dsqrt(difference) 521 continue sigma_phi(1)=(abs(parmavg(9))*sigma(10) 1 + abs(parmavg(10))*sigma(9))/rsq(1) sigma_phi(2)=(abs(parmavg(11))*sigma(12) 1 + abs(parmavg(12))*sigma(11))/rsq(2) sigma_phi(3)=(abs(parmavg(13))*sigma(14) 1 + abs(parmavg(14))*sigma(13))/rsq(3) sigma(3)=sigma_phi(1)*angle if(sigma(3).gt.9999.90) then sigma(6) = 9999.9 write(6,*)'Sigma(3) overflow before ',iyymmdd,sec, 1 icnt,' rcrds in.' endif sigma(6)=sigma_phi(2)*angle if(sigma(6).gt.9999.90) then sigma(6) = 9999.9 write(6,*)'Sigma(3) overflow before ',iyymmdd,sec, 1 icnt,' rcrds in.' endif sigma(8)=sigma_phi(3)*angle if(sigma(8).gt.9999.90) then sigma(8) = 9999.9 write(6,*)'Sigma(3) overflow before ',iyymmdd,sec, 1 icnt,' rcrds in.' endif C Now write output; averages and number of samples in each average. write(7,2112,err=3456)(itime,hour,(parmavg(j),j=1,2),phiavg(1), 1 (parmavg(j),j=4,5),phiavg(2),parmavg(7),phiavg(3), 2 (sigma(n),n=1,8),(parmcount(n),n=1,8),rsq) 2112 format(1x,I6,1x,I2,2(1PE10.3),0PF6.1,2(1PE10.3),0PF6.1, 1 1PE10.3,0PF6.1, 2 2(1PE10.3),0PF6.1,2(1PE10.3),0PF6.1, 3 1PE10.3,0PF6.1,8I3,3F6.3) icnt7=icnt7+1 C Remember the new time itime=iyymmdd time=utime uhour=utime/10000. hour=uhour C Zero the counters do 221 k=1,14 parmcount(k)=0 parmsum(k)=0. sumsquare(k)=0. sq_avg(k)=-9.999 parmavg(k)=-9.999 221 continue if(iend.eq.1) go to 300 go to 4 3007 write(6,*)' Problem opening output file. IOSTAT=',ios stop 3010 write(6,*)' Problem opening input file. IOSTAT=',ios go to 3001 330 continue write(6,*) ' EOF found while seeking first record of new file.' stop 1330 write(6,*) ' Error found while seeking first record of new file.' stop 331 iend=1 !EOF; process current data, ask for next infile and start again. go to 350 1331 write(6,*)' Input Err after ',icnt,' total input records.' stop ! MAY NOT WANT TO STOP HERE****** c go to 2 ! now set to read again 333 continue write(6,*) 'Time changed backwards, at ',itime,'; record ',icthis stop 3456 continue write(6,*) 'Write error. ',iyymmdd,sec,icnt,' Input records read.' stop 30 continue stop end