program extract c c program used to extract part(s) of the large files of the Neckel c full disk and disk center itensity atlas c c The data files are supposed to have the followwing names: c for the full disk intensity : neckelfd1 to neckelfd10 c for the disk center intensity : neckeli1 to neckeli10 c c neckeli1 -> neckeli_3290 c neckeli2 -> neckeli_4000 c neckeli3 -> neckeli_5000 c neckeli4 -> neckeli_6000 c neckeli5 -> neckeli_7000 c neckeli6 -> neckeli_8000 c 4862040 Jul 16 12:16 neckeli_3290 c 5934206 Jul 16 12:33 neckeli_4000 c 4573151 Jul 16 12:34 neckeli_5000 c 4064968 Jul 16 13:01 neckeli_6000 c 3468745 Jul 16 13:44 neckeli_7000 has 111895 lines c 3073371 Jul 16 14:42 neckeli_8000 has 99141 lines c c Not all data files need be present. If a file is needed but not present c an error is reported. c c In extracting you have the choice between full disk and disk center c It is possible to extract a part of the spectrum that is covered in more c than one file. It runs over the boundaries of the files. c c remark about units: c The units in the inputfiles are W cm^-2 sr^-1 A^-1 (I_\lambda) c The more regular units are erg s^-1 cm^-2 sr^-1 Hz^-1 (I_\nu) c Therefore give the more regular values in the converted files. c The conversion formula takes care of all differences of Angstrom, c cm, Watt and erg. c parameter(mp=400000) character ch*1, filein*30, fileout*30 real wave(mp),i(mp),cnt(mp) real i2(mp),cnt2(mp) real wbfd(10),wbi(10),wefd(10),wei(10) real wb(10),we(10) c data wbfd/3290.0002,4000.0044,5000.0029,6000.0061,7000.0002, * 8000.0057,9000.0083,10000.0069,11000.0024,12000.0075/ data wbi/3290.0013,4000.0041,5000.0044,6000.0026,7000.0037, * 8000.0006,9000.0041,10000.0080,11000.0013,12000.0063/ data wefd/3999.9997,4999.9969,5999.9996,6999.9913,7999.9966, * 8999.9968,9999.9927,10999.9884,11999.9909,12509.9836/ data wei/3999.9995,4999.9984,5999.9960,6999.9948,7999.9916, * 8999.9948,9999.9964,10999.9873,11999.9897,12509.9824/ data c/2.99792458e10/ 3 continue print *,' Do you want to do full file conversions or just' print *,' extract part(s) of file(s) (f/p) p ' read(5,'(a)') ch if (ch.eq.'') ch='p' if (ch.ne.'f' .and. ch.ne.'F') then ch='p' else ch='f' endif if (ch.eq.'f') then c c Full file conversions c goto 98 97 print *,' Input file not present, try again' 98 print *,' Enter input file name ' read(5,'(a)') filein print *,' Enter output file name ' read(5,'(a)') fileout c open(unit=7,file=filein,form='formatted',status='old',err=97) open(unit=8,file=fileout,form='formatted',status='new') c do 10 k=1,mp read(7,*,end=11) wave(k),i(k),cnt(k) i2(k)=wave(k)**2.*0.1*i(k)/c cnt2(k)=wave(k)**2.*0.1*cnt(k)/c write(8,100) wave(k),i2(k),cnt2(k) 10 continue 100 format(0pf12.4,1pe12.4,e12.4) 11 close(7) close(8) else print *,' Enter output file name for extracted data ' read(5,'(a)') fileout open(unit=8,file=fileout,form='formatted',status='new') 4 print *,' Disk center or full disk intensities ? (c/f) c ' read(5,'(a)') ch if (ch.eq.'') ch='c' if (ch.ne.'c' .and. ch.ne.'f') then print *,'wrong option, try again' goto 4 endif goto 7 99 print *,' Needed input file not present, try again' 7 print *,' Enter start wavelength in Angstrom' read(5,*) start print *,' Enter end wavelength in Angstrom' read(5,*) eind nr=0 if (eind.lt.start) then nr=1 c c can be used in output routine to reverse the order c tus=start start=eind eind=tus endif c c Find out which files are needed c if (ch.eq.'c') then do 50 k=1,10 wb(k)=wbi(k) we(k)=wei(k) 50 continue else do 60 k=1,10 wb(k)=wbfd(k) we(k)=wefd(k) 60 continue endif c c Check all cases for begin and end wavelength location c if (eind.lt.wb(1) .or. start.gt.we(10)) stop 'waves out of range' if (start.lt.wb(1)) start=wb(1) if (eind.gt.we(10)) eind=we(10) do 70 k=1,10 if (start.gt.wb(k) .and. start.lt.we(k)) nst=k if (eind.gt.wb(k) .and. eind.lt.we(k)) ne=k 70 continue nsw=0 if ((ne-nst).gt.1) stop 'wave range too large' c c make input file name, applicable for both cases c if (ch.eq.'c') then if (nst.lt.10) then write(filein,110) nst else write(filein,111) nst endif else if (nst.lt.10) then write(filein,112) nst else write(filein,113) nst endif endif print *,filein 110 format('neckeli',i1) 111 format('neckeli',i2) 112 format('neckelfd',i1) 113 format('neckelfd',i2) open(unit=7,file=filein,form='formatted',status='old',err=99) c c Read complete file c do 80 k=1,mp read(7,*,end=81) wave(k),i(k),cnt(k) 80 continue 81 nb=k close(7) if (ne.eq.nst) then np=nb-1 else if (ch.eq.'c') then if (ne.lt.10) then write(filein,110) ne else write(filein,111) ne endif else if (ne.lt.10) then write(filein,112) ne else write(filein,113) ne endif endif print *,filein c c read on next file c open(unit=7,file=filein,form='formatted',status='old') do 90 k=np,mp read(7,*,end=91) wave(k),i(k),cnt(k) 90 continue 91 np=k-1 endif print *,' Input file(s) read' c c Select the needed points c do 200 k=1,np if (start.gt.wave(k)) nstrt=k if (eind.gt.wave(k)) neind=k+1 200 continue print *,' nstrt = ',nstrt,' neind = ',neind do 300 k=nstrt,neind i2(k)=wave(k)**2.*0.1*i(k)/c cnt2(k)=wave(k)**2.*0.1*cnt(k)/c write(8,100) wave(k),i2(k),cnt2(k) 300 continue close(7) print *,' Do you want to extract some more onto the same ' print *,' output file (y/n) n' read(5,'(a)') ch if (ch.eq.'') ch='n' if (ch.eq.'y') goto 4 close(8) endif print *,' Do you want to convert or extract some more? (y/n) y ' read(5,'(a)') ch if (ch.eq.'') ch='y' if (ch.eq.'y' .or. ch.eq.'Y') goto 3 stop end