*********************************************************************** * helioglib62 - heliog library v6.2 * * 17/34 GHz Dual Frequency Observation * *********************************************************************** *********************************************************************** * rdfrm17 - read 1-frame raw data of 17 GHz * * and convert to 84x85 array * * 1-sec data, 600-frame event data, event data * * iunit : input,integer fortran file unit number * * for raw data file * * filehd : input,character*80 raw data filename header * * in the case of event-data, raw data filename * * nframe : input,integer frame number * * frmhd : output,character*16 frame header * * ndjst : output,integer observing date in jst (yyyymmdd) * * msjst : output,observing time in jst (hhmmsssss) * * natt : output,integer attenuator in db * * nfalt : output,integer frequency alternation * * 0 : fix, 1 : alternate * * nfreq : output,integer observing frequency in GHz * * nstat : output,integer frame status * * 0 : not recorded * * 1 : calibration data * * 2 : solar observation data * * 3 : event data * * 15 : re-transferred data * * icorrc : output,integer*2(84,85) rcp/cos uv-map * * icorrs : output,integer*2(84,85) rcp/sin uv-map * * icorlc : output,integer*2(84,85) lcp/cos uv-map * * icorls : output,integer*2(84,85) lcp/sin uv-map * * * * icurf : the file number which is currently opend * * icurr : the record number to be read by the next read * * statement * * September 7, 1998 Y. Hanaoka * *********************************************************************** c subroutine rdfrm17(iunit,filehd,nframe, - frmhd,ndjst,msjst,natt,nfalt,nfreq,nstat, - icorrc,icorrs,icorlc,icorls) c character*80 infile,filehd character*3 fnum character*28 infohd character*600 dflag character*16 frmhd character*8 tchar character*2 ext1 character*6 ext2 character*1 bcdc character*1 hex(0:15)/'0','1','2','3','4','5','6','7', - '8','9','a','b','c','d','e','f'/ integer nt(16) integer*2 iraw(14112) integer*2 icorrc(84,85),icorrs(84,85), - icorlc(84,85),icorls(84,85) integer icurf/0/ c ***** function bcd -> integer ***** ibcd(bcdc)=ichar(bcdc)/16*10+min(9,mod(ichar(bcdc),16)) c ***** file-type test ***** if (icurf.eq.0) then open(iunit,file=filehd, - access='direct',recl=28248,status='old',iostat=nfstat) if (nfstat.ne.0) then write(6,'('' data-file is 600-frame type'')') else write(6,'('' data-file is event-data type'')') c write(6,'('' open file :'',a40)') filehd write(6,'('' open file :'',a60)') filehd icurf=1 endif endif c ***** initialize ***** if (nfstat.ne.0) then ifile=(nframe-1)/600+1 iframe=nframe-(ifile-1)*600 if (icurf.ne.ifile) then if (icurf.ne.0) close(iunit) icurf=ifile write(fnum,'(i3.3)') icurf ilfl=index(filehd,' ')-1 infile=filehd(1:ilfl)//fnum open(iunit,file=infile, - access='direct',recl=28248,status='old') read(iunit,rec=1) infohd,dflag c write(6,'('' open file :'',a40)') infile write(6,'('' open file :'',a60)') infile end if c nstat=ichar(dflag(iframe:iframe))-48 if (nstat.eq.0) then write(6,'('' nframe ='',i6,'' nstat ='',i6)') nframe,nstat return endif nrec=iframe+1 else nstat=3 nrec=nframe end if c ***** read frame header and data ***** c read(iunit,rec=nrec) nfrm,frmhd,iraw if (nfrm .ne. nframe) then write(6,'('' frame number error : frame ='',i6,'' -> '',i6)') - nframe,nfrm endif tchar=frmhd(1:8) ext1=frmhd(9:10) ext2=frmhd(11:16) nyr=ibcd(tchar(1:1)) nmo=ibcd(tchar(2:2)) ndt=ibcd(tchar(3:3)) nhr=ibcd(tchar(4:4)) nmin=ibcd(tchar(5:5)) nmsec=ibcd(tchar(6:6))*1000 - +ibcd(tchar(7:7))*10+min(9,ichar(tchar(8:8))/16) c write(6,'('' frame :'',i6,'' dd/mm/yy/hh:mm:ss.sss : '', c - i2.2,''/'',i2.2,''/'',i2.2,'' '', c - i2.2,'':'',i2.2,'':'',f6.3)') c - nfrm,nyr,nmo,ndt,nhr,nmin,real(nmsec)/1000. c do 300 i=1,8 nt(i*2-1)=ichar(tchar(i:i))/16 nt(i*2) =mod(ichar(tchar(i:i)),16) 300 continue nbcder=0 do 310 i=1,15 310 if (nt(i) .gt. 9) nbcder=1 if (nbcder.eq.1) then write(6,'('' time stamp error : frame ='',i6,'' '',15a1, - '' -> '',i2.2,i2.2,i2.2,i2.2,i2.2,f6.3)') - nframe,(hex(nt(i)),i=1,15), - nyr,nmo,ndt,nhr,nmin,real(nmsec)/1000. end if c if (nyr .gt. 90) then nyr=nyr+1900 else nyr=nyr+2000 endif ndjst=nyr*10000+nmo*100+ndt msjst=nhr*10000000+nmin*100000+nmsec c if (mod(ichar(ext2(6:6)),2).eq.0) then natt=0 else natt=10 end if nfalt=mod(nt(16)/2,2) if (mod(nt(16),2).eq.1) then nfreq=34 else nfreq=17 end if c ***** create 84x85 arrays ***** call arrng(iraw,icorrc,icorrs) call arrng(iraw(7057),icorlc,icorls) c return end c c *********************************************************************** * arrng - arrange raw 1-frame data to 84x85 arrays * * iraw : input,integer(14112) raw 1-frame data * * icorc : output,integer*2(84,85) correlation (cos) and * * intensity * * icors : output,integer*2(84,85) correlation (sin) * * correlation(n,m) : n,m = 1-84 or n01-n28/e28-e01/w01-w28 * * icorc(n,n)=0. icors(n,n)=0. * * (icorc(n,n)=32760. icors(n,n)=0.) * * intensity icorc(n,85) : n = 1-84 or n01-n28/e28-e01/w01-w28 * * October 10,1995 Y. Hanaoka * *********************************************************************** c subroutine arrng(iraw,icorc,icors) c integer*2 iraw(7056) integer*2 icorc(84,85),icors(84,85) c ***** flux ***** do 900 i=1,84 900 icorc(i,85)=iraw(i) c ***** N - N ***** n=0 do 100 j=1,27 do 100 i=j+1,28 n=n+1 icorc(i,j)=iraw(3136+n*2-1+84) 100 continue c n=0 do 120 j=1,27 do 120 i=j+1,28 n=n+1 icors(i,j)=iraw(3136+n*2+84) 120 continue c ***** N - EW ***** n=0 do 200 j=1,28 do 200 i=29,84 n=n+1 icorc(i,j)=iraw(n*2-1+84) 200 continue c n=0 do 220 j=1,28 do 220 i=29,84 n=n+1 icors(i,j)=iraw(n*2+84) 220 continue c ***** EW - EW ***** n=0 do 300 j=29,83 do 300 i=j+1,84 n=n+1 icorc(i,j)=iraw(3892+n*2-1+84) 300 continue c n=0 do 320 j=29,83 do 320 i=j+1,84 n=n+1 icors(i,j)=iraw(3892+n*2+84) 320 continue c ***** remainder ***** do 400 i=1,84 400 icorc(i,i)=0 do 410 i=1,84 410 icors(i,i)=0 c do 420 j=1,83 do 420 i=j+1,84 420 icorc(j,i)=icorc(i,j) do 430 j=1,83 do 430 i=j+1,84 430 icors(j,i)=icors(i,j) c return end *********************************************************************** * rddat17 - read raw data and create 84x85 array * * for 17 GHz data * * iunit : input,integer fortran file unit number * * for raw data file * * filehd : input,character*80 raw data filename header * * nfrst : input,integer first frame number for integration * * nfrend : input,integer last frame number for integration * * nantst : input,integer(84) antenna status * * ndjst : output,integer observing date in jst (yyyymmdd) * * msjst : output,integer center of observing time * * in jst (hhmmsssss) * * msjsts : output,integer start of observing time * * in jst (hhmmsssss) * * msjste : output,integer end of observing time * * in jst (hhmmsssss) * * natt : output,integer attenuator in db * * nfalt : output,integer frequency alternation * * 0 : fix, 1 : alternate * * nfreq : output,integer observing frequency in GHz * * nstat : output,integer frame status * * 0 : not recorded * * 1 : calibration data * * 2 : solar observation data * * 3 : event data * * 15 : re-transferred data * * corrc : output,real(84,85) rcp/cos uv-map * * corrs : output,real(84,85) rcp/sin uv-map * * corlc : output,real(84,85) lcp/cos uv-map * * corls : output,real(84,85) lcp/sin uv-map * * nframe : output,integer number of valid frames * * 0 : all frames are invalid; in this case, other outputs * * are dummy * * June 5, 1998 Y. Hanaoka * *********************************************************************** c subroutine rddat17(iunit,filehd,nfrst,nfrend,nantst, - ndjst,msjst,msjsts,msjste,natt,nfalt,nfreq,nstat, - corrc,corrs,corlc,corls,nframe) c character*80 filehd character*16 frmhd integer nantst(1) integer*2 icorrc(84,85),icorrs(84,85), - icorlc(84,85),icorls(84,85) real corrc(84,85),corrs(84,85),corlc(84,85),corls(84,85) c do 100 j=1,85 do 100 i=1,84 100 corrc(i,j)=0. do 110 j=1,85 do 110 i=1,84 110 corrs(i,j)=0. do 120 j=1,85 do 120 i=1,84 120 corlc(i,j)=0. do 130 j=1,85 do 130 i=1,84 130 corls(i,j)=0. c ***** read data ***** c nframe=0 do 200 iframe=nfrst,nfrend call rdfrm17(iunit,filehd,iframe, - frmhd,ndjst,msjst,natt,nfalt,nfreq,nstat, - icorrc,icorrs,icorlc,icorls) if (nstat.eq.0) goto 200 if (nframe.eq.0) then * start time * msjsts=msjst mdjsts=msjst/10000000*3600000+mod(msjst,10000000)/100000*60000 - +mod(msjst,100000) end if * end time * nhr=msjst/10000000 nmin=mod(msjst,10000000)/100000 * 1 sec data : +1000 msec, event data : +50 msec * if (nstat.eq.3) then nmsec=mod(msjst,100000)+50 else nmsec=mod(msjst,100000)+1000 end if if (nmsec.ge.60000) then nmsec=nmsec-60000 nmin=nmin+1 if (nmin.eq.60) then nmin=0 nhr=nhr+1 end if end if msjste=nhr*10000000+nmin*100000+nmsec mdjste=nhr*3600000+nmin*60000+nmsec c nframe=nframe+1 c ***** integration ***** c do 220 j=1,84 do 220 i=1,84 220 corrc(i,j)=corrc(i,j)+real(icorrc(i,j)) do 230 j=1,85 do 230 i=1,84 230 corrs(i,j)=corrs(i,j)+real(icorrs(i,j)) do 240 j=1,84 do 240 i=1,84 240 corlc(i,j)=corlc(i,j)+real(icorlc(i,j)) do 250 j=1,85 do 250 i=1,84 250 corls(i,j)=corls(i,j)+real(icorls(i,j)) c * the values of flux are expressed in unsigned 2-byte integer * do 260 i=1,84 if (icorrc(i,85).ge.0) then corrc(i,85)=corrc(i,85)+real(icorrc(i,85)) else corrc(i,85)=corrc(i,85)+real(icorrc(i,85))+65536. end if 260 continue do 270 i=1,84 if (icorlc(i,85).ge.0) then corlc(i,85)=corlc(i,85)+real(icorlc(i,85)) else corlc(i,85)=corlc(i,85)+real(icorlc(i,85))+65536. end if 270 continue c 200 continue c ***** normalize ***** c if (nframe.eq.0) then write(6,'('' there is no valid frame'')') return endif framen=real(nframe) do 300 j=1,85 do 300 i=1,84 300 corrc(i,j)=corrc(i,j)/framen do 310 j=1,85 do 310 i=1,84 310 corrs(i,j)=corrs(i,j)/framen do 320 j=1,85 do 320 i=1,84 320 corlc(i,j)=corlc(i,j)/framen do 330 j=1,85 do 330 i=1,84 330 corls(i,j)=corls(i,j)/framen c ***** 0-set for correlations of unuse antennas ***** c do 510 j=1,84 if (nantst(j).eq.1) then do 520 i=1,84 corrc(i,j)=0. corrs(i,j)=0. corlc(i,j)=0. corls(i,j)=0. corrc(j,i)=0. corrs(j,i)=0. corlc(j,i)=0. corls(j,i)=0. 520 continue end if 510 continue c ***** observation time ***** c msjst=(mdjsts+mdjste)/2 msjst=msjst/3600000*10000000+mod(msjst,3600000)/60000*100000 - +mod(msjst,60000) c c write(6,'('' '',8i10)') ndjst,msjst,msjsts,msjste, c - natt,nfalt,nfreq,nstat return end *********************************************************************** * snapuv17 - create 2-dimensional uv-map from raw data * * and calibrate * * iunit : input,integer fortran file unit number * * for raw data file * * filehd : input,character*80 raw data filename header * * nfrst : input,integer first frame number for integration * * nfrend : input,integer last frame number for integration * * nantst : input,integer(84) antenna status * * ical : input,integer calibration menu * * ical=999 : no calibration * * ical=0 : use icald as calibration data * * ical<>0 : calculate calibration data integrating * * abs(ical) frames * * if ical<0, calculated antenna phases are * * shifted, using icald as reference data * * icald : input,integer*2(336) initial calibration data * * : output obtained calibration data * * if integrate frames, cal data for the last frame * * rlph : input,real(84) instrumental phase differences of * * antennas * * if rlph(1)=-1, no phase correction is done * * if rlph(1)=1., phase differences are initialized * * taper : input,real taper at 160d * * ndjst : output,integer observing date in jst (yyyymmdd) * * msjst : output,integer center of observing time * * in jst (hhmmsssss) * * msjsts : output,integer start of observing time * * in jst (hhmmsssss) * * msjste : output,integer end of observing time * * in jst (hhmmsssss) * * natt : output,integer attenuator in db * * nfalt : output,integer frequency alternation * * 0 : fix, 1 : alternate * * nfreq : output,integer observing frequency in GHz * * nstat : output,integer frame status * * 0 : not recorded * * 1 : calibration data * * 2 : solar observation data * * 3 : event data * * 15 : re-transferred data * * flxr : output,real(84) rcp total flux * * flxl : output,real(84) lcp total flux * * sumra : output,real(512,512) rcp/cos uv-map * * sumrb : output,real(512,512) rcp/sin uv-map * * sumla : output,real(512,512) lcp/cos uv-map * * sumlb : output,real(512,512) lcp/sin uv-map * * nfrmst : output,integer status * * 0 : all frames have been read * * 1 : there is an invalid frame; when nfrmst=1, other output * * is not valid * * November 27, 1998 Y. Hanaoka * *********************************************************************** c subroutine snapuv17(iunit,filehd, - nfrst,nfrend,nantst,ical,icald,rlph,taper, - ndjst,msjst,msjsts,msjste,natt,nfalt,nfreq,nstat, - flxr,flxl,sumra,sumrb,sumla,sumlb,nfrmst) c character*80 filehd integer nantst(1) integer*2 icald(336),icald0(336) real corrc(84,85),corrs(84,85),corlc(84,85),corls(84,85), - cora(512,512),corb(512,512), - flxr(84),flxl(84), - sumra(512,512),sumrb(512,512),sumla(512,512),sumlb(512,512) real calib(168),rlph(1) c c do 300 j=1,512 do 300 i=1,512 sumra(i,j)=0. sumrb(i,j)=0. sumla(i,j)=0. sumlb(i,j)=0. 300 continue do 310 i=1,84 flxr(i)=0. flxl(i)=0. 310 continue do 320 i=1,336 320 icald0(i)=icald(i) c do 350 iframe=nfrst,nfrend call rddat17(iunit,filehd,iframe,iframe,nantst, - ndjst,msjst,jsts,jste,natt,nfalt,nfreq,nstat, - corrc,corrs,corlc,corls,nframe) c ***** if there is a no data frame, exit ***** if (nframe.eq.0) then nfrmst=1 go to 900 end if c ***** start and end time of integration ***** if (iframe.eq.nfrst) then msjsts=jsts mdjsts=jsts/10000000*3600000+mod(jsts,10000000)/100000*60000 - +mod(jsts,100000) end if c if (iframe.eq.nfrend) then msjste=jste mdjste=jste/10000000*3600000+mod(jste,10000000)/100000*60000 - +mod(jste,100000) end if c ***** get calibration data and calculate r-l phase shift ***** dcns=0. dcew=0. rlns=0. rlew=0. if (ical.ne.999) then if (ical.ne.0) then call rddcal17(iunit,filehd,iframe, - iabs(ical),nantst,icald,nframe) if (ical.lt.0) then call dcshft(icald0,icald,dcns,dcew) ccc write(6,'('' disk center shift ='',2f10.2)') dcns,dcew end if c do 580 i=1,336 c 580 icald(i)=icaldc(i) end if if (rlph(1).ne.-1.) then if (rlph(1).eq.1.) call gtrlph(ndjst,rlph) c write(6,'('' '',7f10.2)') rlph call phshft(rlph,icald,rlns,rlew) ccc write(6,'('' phase shift ='',2f10.2)') rlns,rlew end if end if c c ***** create and calibrate 2d uv-map (rcp) ***** call fi2d(corrc,corrs,cora,corb) if (ical.ne.999) then do 600 i=1,84 600 calib(i)=real(icald(i))/100. do 610 i=85,168 610 calib(i)=real(icald(i))/10000. call addapa(cora,corb,calib,dcns,dcew) end if c do 360 i=1,84 360 flxr(i)=flxr(i)+corrc(i,85) do 370 j=1,512 do 370 i=1,512 sumra(i,j)=sumra(i,j)+cora(i,j) sumrb(i,j)=sumrb(i,j)+corb(i,j) 370 continue c ***** create and calibrate 2d uv-map (lcp) ***** call fi2d(corlc,corls,cora,corb) if (ical.ne.999) then do 620 i=1,84 620 calib(i)=real(icald(168+i))/100. do 630 i=85,168 630 calib(i)=real(icald(168+i))/10000. call addapa(cora,corb,calib,dcns+rlns,dcew+rlew) end if c ***** integration ***** do 380 i=1,84 380 flxl(i)=flxl(i)+corlc(i,85) do 390 j=1,512 do 390 i=1,512 sumla(i,j)=sumla(i,j)+cora(i,j) sumlb(i,j)=sumlb(i,j)+corb(i,j) 390 continue c 350 continue c c ***** center of the integration time ***** msjst=(mdjsts+mdjste)/2 msjst=msjst/3600000*10000000+mod(msjst,3600000)/60000*100000 - +mod(msjst,60000) c ***** normalize ***** framen=real(nfrend-nfrst+1) do 700 i=1,84 flxr(i)=flxr(i)/framen flxl(i)=flxl(i)/framen 700 continue do 710 j=1,512 do 710 i=1,512 sumra(i,j)=sumra(i,j)/framen sumrb(i,j)=sumrb(i,j)/framen sumla(i,j)=sumla(i,j)/framen sumlb(i,j)=sumlb(i,j)/framen 710 continue c call filter1(taper,sumra) call filter1(taper,sumrb) call filter1(taper,sumla) call filter1(taper,sumlb) c nfrmst=0 900 return end c c *********************************************************************** * rddcal17 - read raw heliograph data and calibrate * * iunit : input,integer fortran file unit number * * for raw data file * * filehd : input,character*80 raw data filename header * * iframe : input,integer frame number for calibration * * itgr : input,integer number of frames for integration * * nantst : input,integer(84) antenna status * * icald : output,integer*2(336) antenna phases and gains * * nframe : output,integer number of frames which is actually * * used * * June 5, 1998 Y. Hanaoka * *********************************************************************** c subroutine rddcal17(iunit,filehd,iframe,itgr,nantst, - icald,nframe) c character*80 filehd integer nantst(84) real corrc(84,85),corrs(84,85),corlc(84,85),corls(84,85) c integer iposns(28),iposew(57),icv(9) real phansr(200),phansl(200),phaewr(200),phaewl(200), - amansr(200),amansl(200),amaewr(200),amaewl(200) integer*2 icald(336) c data iposns/0,1,2,3,4,5,6,7,8,9,10,11, - 12,14,16,20,24,28, - 32,40,48,56, - 64,80,96,112,128,144/ data iposew/-160,-144,-128,-112,-96,-80,-64, - -56,-48,-40,-32, - -28,-24,-20,-16,-14,-12, - -11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1, - 0,1,2,3,4,5,6,7,8,9,10,11, - 12,14,16,20,24,28, - 32,40,48,56, - 64,80,96,112,128,144,160/ data icv/1,2,4,6,8,12,16,24,32/ c nantns=28 nantew=57 ncv=9 c call rddat17(iunit,filehd,iframe-(itgr-1)/2, - iframe+itgr/2,nantst, - ndjst,msjst,jsts,jste,natt,nfalt,nfreq,nstat, - corrc,corrs,corlc,corls,nframe) c write(6,'('' nframe ='',i6)') nframe c c inic=iframe-nfrst inic=0 call phcal(nantns,iposns,ncv,icv,1,corrc,corrs,inic,phansr, - sums,nmax) ccc write(6,'('' rcp/ph/ns-cal : error ='',f12.4,'' n ='',i8)') ccc - sums,nmax call phcal(nantew,iposew,ncv,icv,2,corrc,corrs,inic,phaewr, - sums,nmax) ccc write(6,'('' rcp/ph/ew-cal : error ='',f12.4,'' n ='',i8)') ccc - sums,nmax call amcal(nantns,iposns,ncv,icv,1,corrc,corrs,inic,amansr, - sums,nmax) ccc write(6,'('' rcp/am/ns-cal : error ='',f12.4,'' n ='',i8)') ccc - sums,nmax call amcal(nantew,iposew,ncv,icv,2,corrc,corrs,inic,amaewr, - sums,nmax) ccc write(6,'('' rcp/am/ew-cal : error ='',f12.4,'' n ='',i8)') ccc - sums,nmax c call phcal(nantns,iposns,ncv,icv,1,corlc,corls,inic,phansl, - sums,nmax) ccc write(6,'('' lcp/ph/ns-cal : error ='',f12.4,'' n ='',i8)') ccc - sums,nmax call phcal(nantew,iposew,ncv,icv,2,corlc,corls,inic,phaewl, - sums,nmax) ccc write(6,'('' lcp/ph/ew-cal : error ='',f12.4,'' n ='',i8)') ccc - sums,nmax call amcal(nantns,iposns,ncv,icv,1,corlc,corls,inic,amansl, - sums,nmax) ccc write(6,'('' lcp/am/ns-cal : error ='',f12.4,'' n ='',i8)') ccc - sums,nmax call amcal(nantew,iposew,ncv,icv,2,corlc,corls,inic,amaewl, - sums,nmax) ccc write(6,'('' lcp/am/ew-cal : error ='',f12.4,'' n ='',i8)') ccc - sums,nmax c icald(1)=0 do 600 i=2,28 600 icald(i)=int(phansr(i-1)*100.) do 610 i=1,56 610 icald(i+28)=int(phaewr(i)*100.) icald(85)=10000 c do 620 i=2,28 icald(i+84)=int(amansr(i-1)*10000.) if (nantst(i).eq.1) icald(i+84)=0 620 continue do 630 i=1,28 icald(i+112)=int(amaewr(i)*10000.) if (nantst(i+56).eq.1) icald(i+112)=0 630 continue do 640 i=1,28 icald(i+140)=int(amaewr(i+28)*10000.) if (nantst(57-i).eq.1) icald(i+140)=0 640 continue c icald(169)=0 do 700 i=2,28 700 icald(i+168)=int(phansl(i-1)*100.) do 710 i=1,56 710 icald(i+196)=int(phaewl(i)*100.) c icald(253)=10000 do 720 i=2,28 icald(i+252)=int(amansl(i-1)*10000.) if (nantst(i).eq.1) icald(i+252)=0 720 continue do 730 i=1,28 icald(i+280)=int(amaewl(i)*10000.) if (nantst(i+56).eq.1) icald(i+280)=0 730 continue do 740 i=1,28 icald(i+308)=int(amaewl(i+28)*10000.) if (nantst(57-i).eq.1) icald(i+308)=0 740 continue c return end *********************************************************************** * fi2d : create 2-dimensional uv-map(512x512) from * * whole correlation array(84x85) * * corc : input,real(84,84) cos of all correlations * * cors : input,real(84,85) sin of all correlations * * cora : output,real(512,512) cos of 2-d correlations * * corb : output,real(512,512) sin of 2-d correlations * * June 22,1992 Y. Hanaoka * *********************************************************************** c subroutine fi2d(corc,cors,cora,corb) c integer ipos(-28:28) real corc(84,85),cors(84,85),cora(512,512),corb(512,512) data ipos/-160,-144,-128,-112,-96,-80,-64, - -56,-48,-40,-32, - -28,-24,-20,-16,-14,-12, - -11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1, - 0,1,2,3,4,5,6,7,8,9,10,11, - 12,14,16,20,24,28, - 32,40,48,56, - 64,80,96,112,128,144,160/ c do 100 l=1,512 do 100 i=1,512 100 cora(i,l)=0. do 110 l=1,512 do 110 i=1,512 110 corb(i,l)=0. c ***** create lower half of uv map (cos) ***** do 200 l=-27,-1 do 200 i=-28,28 iy=iabs(l)+1 if (i.lt.0) then cora(-ipos(i)+257,ipos(l)+257)=corc(-i+56,iy) else if (i.eq.0) then cora(-ipos(i)+257,ipos(l)+257)=corc(1,iy) else cora(-ipos(i)+257,ipos(l)+257)=corc(-i+57,iy) end if 200 continue c ***** create center of uv map (cos) ***** cora(257,257)=corc(1,1) do 210 i=-28,-1 c cora(-ipos(i)+257,257)=corc(i+57,1) cora(-ipos(i)+257,257)=corc(-i+56,1) 210 continue do 220 i=1,28 cora(-ipos(i)+257,257)=corc(i+56,1) 220 continue c ***** create upper half of uv map (cos) ***** do 230 l=1,27 do 230 i=-28,28 iy=iabs(l)+1 if (i.lt.0) then ix=i+57 cora(-ipos(i)+257,ipos(l)+257)=corc(i+57,iy) else if (i.eq.0) then cora(-ipos(i)+257,ipos(l)+257)=corc(1,iy) else cora(-ipos(i)+257,ipos(l)+257)=corc(i+56,iy) end if 230 continue c ***** create lower half of uv map (sin) ***** do 250 l=-27,-1 do 250 i=-28,28 iy=iabs(l)+1 if (i.lt.0) then corb(-ipos(i)+257,ipos(l)+257)=-cors(-i+56,iy) else if (i.eq.0) then corb(-ipos(i)+257,ipos(l)+257)=cors(1,iy) else corb(-ipos(i)+257,ipos(l)+257)=-cors(-i+57,iy) end if 250 continue c ***** create center of uv map (sin) ***** corb(257,257)=cors(1,1) do 260 i=-28,-1 c corb(-ipos(i)+257,257)=cors(i+57,1) corb(-ipos(i)+257,257)=-cors(-i+56,1) 260 continue do 270 i=1,28 corb(-ipos(i)+257,257)=cors(i+56,1) 270 continue c ***** create upper half of uv map (sin) ***** do 280 l=1,27 do 280 i=-28,28 iy=iabs(l)+1 if (i.lt.0) then corb(-ipos(i)+257,ipos(l)+257)=cors(i+57,iy) else if (i.eq.0) then corb(-ipos(i)+257,ipos(l)+257)=-cors(1,iy) else corb(-ipos(i)+257,ipos(l)+257)=cors(i+56,iy) end if 280 continue c return end c *********************************************************************** * addapa - add antenna phases and gains to 2-d uv map * * cora : input,real(512,512) raw uv map (cos) * * output,real(512,512) calibrated uv map (cos) * * corb : input,real(512,512) raw uv map (sin) * * output,real(512,512) calibrated uv map (sin) * * calib : input,real(168) calibration data * * calib( 1- 28) antenna phase n01-n28 (degree) * * calib( 29- 56) antenna phase w01-w28 (degree) * * calib( 57- 84) antenna phase e01-e28 (degree) * * calib( 85-112) antenna gain n01-n28 * * calib(113-140) antenna gain w01-w28 * * calib(141-168) antenna gain e01-e28 * * rlns : input,real r-l phase difference(n-s) * * rlew : input,real r-l phase difference(e-w) * * November 3,1992 Y. Hanaoka * *********************************************************************** c subroutine addapa(cora,corb,calib,rlns,rlew) c integer ipos(29) real cora(512,512),corb(512,512),calib(1), - pee(28),pew(28),pen(28),aee(28),aew(28),aen(28) data ipos/0,1,2,3,4,5,6,7,8,9,10,11,12, - 14,16,20,24,28,32,40,48,56,64, - 80,96,112,128,144,160/ c pi=3.14159 rd=57.2958 c ***** separate calibration data ***** do 100 i=1,28 100 pen(i)=calib(i) do 110 i=1,28 110 pew(i)=calib(i+28) do 120 i=1,28 120 pee(i)=calib(i+56) c do 150 i=1,28 c 150 aen(i)=1. aen(i)=calib(i+84) if (aen(i).lt.1.e-4) aen(i)=1. 150 continue do 160 i=1,28 c 160 aew(i)=1. aew(i)=calib(i+28+84) if (aew(i).lt.1.e-4) aew(i)=1. 160 continue do 170 i=1,28 c 170 aee(i)=1. aee(i)=calib(i+56+84) if (aee(i).lt.1.e-4) aee(i)=1. 170 continue c c***** gain offset ***** c offset=0. c do 200 i=2,28 c 200 offset=offset+alog(aen(i)) c do 210 i=1,28 c 210 offset=offset+alog(aee(i)) c do 220 i=1,28 c 220 offset=offset+alog(aew(i)) c offset=exp(offset/84.) c write(6,'('' amplitude offset ='',f10.4)') offset offset=1. c ***** calibration ***** c do 300 il=1,28 c ***** calibrate N-N correlation ***** l=257+ipos(il) am=sqrt(cora(257,l)*cora(257,l)+corb(257,l)*corb(257,l)) if (am.lt.1.e-4) then ph=0. else ph=acos(cora(257,l)/am) if (corb(257,l).lt.0.) ph=2*pi-ph end if c ph=ph-pen(il)/rd ph=ph-pen(il)/rd-rlns*real(ipos(il))/rd am=am/aen(il)*offset*offset cora(257,l)=am*cos(ph) corb(257,l)=am*sin(ph) cora(257,514-l)=cora(257,l) corb(257,514-l)=-corb(257,l) c ***** calibrate N-W correlation ***** do 320 ix=1,28 i=257-ipos(ix+1) am=sqrt(cora(i,l)*cora(i,l)+corb(i,l)*corb(i,l)) if (am.lt.1.e-4) then ph=0. else ph=acos(cora(i,l)/am) if (corb(i,l).lt.0.) ph=2*pi-ph end if c ph=ph+(pew(ix)-pen(il))/rd ph=ph+(pew(ix)-pen(il))/rd-rlew*real(-ipos(ix+1))/rd - -rlns*real(ipos(il))/rd am=am/aew(ix)/aen(il)*offset*offset cora(i,l)=am*cos(ph) corb(i,l)=am*sin(ph) c if (il.ne.1) then cora(514-i,514-l)=cora(i,l) corb(514-i,514-l)=-corb(i,l) c end if 320 continue c ***** calibrate N-E correlation ***** do 340 ix=1,28 if (il.ne.1) then i=257+ipos(ix+1) am=sqrt(cora(i,l)*cora(i,l)+corb(i,l)*corb(i,l)) if (am.lt.1.e-4) then ph=0. else ph=acos(cora(i,l)/am) if (corb(i,l).lt.0.) ph=2*pi-ph end if c ph=ph+(pee(ix)-pen(il))/rd ph=ph+(pee(ix)-pen(il))/rd-rlew*real(ipos(ix+1))/rd - -rlns*real(ipos(il))/rd am=am/aee(ix)/aen(il)*offset*offset cora(i,l)=am*cos(ph) corb(i,l)=am*sin(ph) c if (il.ne.1) then cora(514-i,514-l)=cora(i,l) corb(514-i,514-l)=-corb(i,l) end if 340 continue c c write(6,'('' done '',i2)') il 300 continue c return end c c *********************************************************************** * phshft - calculate rcp-lcp phase shift due to polarization * * rlph : input,real(84) instrumental phase differences of * * antennas * * icald : input,integer(336) calibration data * * rlns : output,real r-l phase difference(n-s) (degree) * * rlew : output,real r-l phase difference(e-w) (degree) * * November 1,1992 Y. Hanaoka * *********************************************************************** c subroutine phshft(rlph,icald,rlns,rlew) c integer*2 icald(1) integer ipos(29) real rlph(1),rlsh(84) data ipos/0,1,2,3,4,5,6,7,8,9,10,11,12, - 14,16,20,24,28,32,40,48,56,64, - 80,96,112,128,144,160/ c ***** N-kei shift ***** sxx=0. sxy=0. rlsh0=0. do 200 i=1,12 if (rlph(i+1).lt.900.) then rlsh(i+1)=(real(icald(i+1))-real(icald(i+1+168)))/100. - -rlph(i+1) do 220 while(abs(rlsh(i+1)-rlsh0).gt.180.) 220 rlsh(i+1)=rlsh(i+1)-sign(360.,rlsh(i+1)-rlsh0) rlsh0=rlsh(i+1) c sxx=sxx+real(ipos(i+1))*real(ipos(i+1)) c sxy=sxy+real(ipos(i+1))*rlsh(i+1) sxx=sxx+real(ipos(i+1)) sxy=sxy+rlsh(i+1) else rlsh(i+1)=999. end if 200 continue rlns=sxy/sxx c sxx=0. sxy=0. do 230 i=1,27 if (rlph(i+1).lt.900.) then rlsh(i+1)=(real(icald(i+1))-real(icald(i+1+168)))/100. - -rlph(i+1) do 240 while(abs(rlsh(i+1)-rlns*real(ipos(i+1))).gt.180.) 240 rlsh(i+1)=rlsh(i+1)-sign(360.,rlsh(i+1)-rlns*real(ipos(i+1))) sxx=sxx+real(ipos(i+1)) sxy=sxy+rlsh(i+1) else rlsh(i+1)=999. end if 230 continue rlns=sxy/sxx c ***** EW-kei shift ***** sxx=0. sxy=0. rlsh0=0. do 250 i=1,12 if (rlph(i+28).lt.900.) then rlsh(i+28)=(real(icald(i+28))-real(icald(i+28+168)))/100. - -rlph(i+28) do 270 while(abs(rlsh(i+28)-rlsh0).gt.180.) 270 rlsh(i+28)=rlsh(i+28)-sign(360.,rlsh(i+28)-rlsh0) rlsh0=rlsh(i+28) c sxx=sxx+real(ipos(i+1))*real(ipos(i+1)) c sxy=sxy+real(ipos(i+1))*rlsh(i+28) sxx=sxx+real(ipos(i+1)) sxy=sxy+rlsh(i+28) else rlsh(i+28)=999. end if 250 continue c rlsh0=0. do 260 i=1,12 if (rlph(i+56).lt.900.) then rlsh(i+56)=(real(icald(i+56))-real(icald(i+56+168)))/100. - -rlph(i+56) do 280 while(abs(rlsh(i+56)-rlsh0).gt.180.) 280 rlsh(i+56)=rlsh(i+56)-sign(360.,rlsh(i+56)-rlsh0) rlsh0=rlsh(i+56) c sxx=sxx+real(-ipos(i+1))*real(-ipos(i+1)) c sxy=sxy+real(-ipos(i+1))*rlsh(i+56) sxx=sxx-real(-ipos(i+1)) sxy=sxy-rlsh(i+56) else rlsh(i+56)=999. end if 260 continue rlew=sxy/sxx c sxx=0. sxy=0. do 350 i=1,28 if (rlph(i+28).lt.900.) then rlsh(i+28)=(real(icald(i+28))-real(icald(i+28+168)))/100. - -rlph(i+28) do 370 while(abs(rlsh(i+28)-rlew*real(ipos(i+1))).gt.180.) 370 rlsh(i+28)=rlsh(i+28)-sign(360., - rlsh(i+28)-rlew*real(ipos(i+1))) sxx=sxx+real(ipos(i+1)) sxy=sxy+rlsh(i+28) else rlsh(i+28)=999. end if 350 continue c do 360 i=1,28 if (rlph(i+56).lt.900.) then rlsh(i+56)=(real(icald(i+56))-real(icald(i+56+168)))/100. - -rlph(i+56) do 380 while(abs(rlsh(i+56)-rlew*real(-ipos(i+1))).gt.180.) 380 rlsh(i+56)=rlsh(i+56)-sign(360., - rlsh(i+56)-rlew*real(-ipos(i+1))) sxx=sxx-real(-ipos(i+1)) sxy=sxy-rlsh(i+56) else rlsh(i+56)=999. end if 360 continue rlew=sxy/sxx c c write(6,'('' '',7f10.2)') rlsh return end c c *********************************************************************** * dcshft - calculate disk center shift * * icald : input,integer(336) initial calibration data * * icaldc : input,integer(336) obtained calibration data * * dcns : output,real phase difference(n-s) (degree) * * dcew : output,real phase difference(e-w) (degree) * * November 5,1992 Y. Hanaoka * *********************************************************************** c subroutine dcshft(icald,icaldc,dcns,dcew) c integer*2 icald(1),icaldc(1) integer ipos(29) real dcsh(84) data ipos/0,1,2,3,4,5,6,7,8,9,10,11,12, - 14,16,20,24,28,32,40,48,56,64, - 80,96,112,128,144,160/ c ***** N-kei shift ***** sxx=0. sxy=0. dcsh0=0. do 200 i=1,12 if (icald(i+1+84).ne.0) then dcsh(i+1)=(real(icald(i+1))-real(icaldc(i+1)))/100. do 220 while(abs(dcsh(i+1)-dcsh0).gt.180.) 220 dcsh(i+1)=dcsh(i+1)-sign(360.,dcsh(i+1)-dcsh0) dcsh0=dcsh(i+1) sxx=sxx+real(ipos(i+1)) sxy=sxy+dcsh(i+1) else dcsh(i+1)=999. end if 200 continue dcns=sxy/sxx c sxx=0. sxy=0. do 230 i=1,27 if (icald(i+1+84).ne.0) then dcsh(i+1)=(real(icald(i+1))-real(icaldc(i+1)))/100. do 240 while(abs(dcsh(i+1)-dcns*real(ipos(i+1))).gt.180.) 240 dcsh(i+1)=dcsh(i+1)-sign(360.,dcsh(i+1)-dcns*real(ipos(i+1))) sxx=sxx+real(ipos(i+1)) sxy=sxy+dcsh(i+1) else dcsh(i+1)=999. end if 230 continue dcns=sxy/sxx c ***** EW-kei shift ***** sxx=0. sxy=0. dcsh0=0. do 250 i=1,12 if (icald(i+28+84).ne.0) then dcsh(i+28)=(real(icald(i+28))-real(icaldc(i+28)))/100. do 270 while(abs(dcsh(i+28)-dcsh0).gt.180.) 270 dcsh(i+28)=dcsh(i+28)-sign(360.,dcsh(i+28)-dcsh0) dcsh0=dcsh(i+28) sxx=sxx+real(ipos(i+1)) sxy=sxy+dcsh(i+28) else dcsh(i+28)=999. end if 250 continue c dcsh0=0. do 260 i=1,12 if (icald(i+56+84).ne.0) then dcsh(i+56)=(real(icald(i+56))-real(icaldc(i+56)))/100. do 280 while(abs(dcsh(i+56)-dcsh0).gt.180.) 280 dcsh(i+56)=dcsh(i+56)-sign(360.,dcsh(i+56)-dcsh0) dcsh0=dcsh(i+56) sxx=sxx-real(-ipos(i+1)) sxy=sxy-dcsh(i+56) else dcsh(i+56)=999. end if 260 continue dcew=sxy/sxx c sxx=0. sxy=0. do 350 i=1,28 if (icald(i+28+84).ne.0) then dcsh(i+28)=(real(icald(i+28))-real(icaldc(i+28)))/100. do 370 while(abs(dcsh(i+28)-dcew*real(ipos(i+1))).gt.180.) 370 dcsh(i+28)=dcsh(i+28)-sign(360., - dcsh(i+28)-dcew*real(ipos(i+1))) sxx=sxx+real(ipos(i+1)) sxy=sxy+dcsh(i+28) else dcsh(i+28)=999. end if 350 continue c do 360 i=1,28 if (icald(i+56+84).ne.0) then dcsh(i+56)=(real(icald(i+56))-real(icaldc(i+56)))/100. do 380 while(abs(dcsh(i+56)-dcew*real(-ipos(i+1))).gt.180.) 380 dcsh(i+56)=dcsh(i+56)-sign(360., - dcsh(i+56)-dcew*real(-ipos(i+1))) sxx=sxx-real(-ipos(i+1)) sxy=sxy-dcsh(i+56) else dcsh(i+56)=999. end if 360 continue dcew=sxy/sxx c c write(6,'('' '',7f10.2)') dcsh return end *********************************************************************** * rdfrm34 - read 1-frame raw data of 34 GHz * * and convert to 84x85 array * * 1-sec data, 600-frame event data, event data * * iunit : input,integer fortran file unit number * * for raw data file * * filehd : input,character*80 raw data filename header * * in the case of event-data, raw data filename * * nframe : input,integer frame number for * * frmhd : output,character*16 frame header * * ndjst : output,integer observing date in jst (yyyymmdd) * * msjst : output,observing time in jst (hhmmsssss) * * natt : output,integer attenuator in db (always 0 for 34 GHz)* * nfalt : output,integer frequency alternation * * 0 : fix, 1 : alternate * * nfreq : output,integer observing frequency in GHz * * nstat : output,integer frame status * * 0 : not recorded * * 1 : calibration data * * 2 : solar observation data * * 3 : event data * * 15 : re-transferred data * * icorc : output,integer*2(84,85) cos uv-map * * icors : output,integer*2(84,85) sin uv-map * * * * icurf : the file number which is currently opend * * icurr : the record number to be read by the next read * * statement * * November 24, 1998 Y. Hanaoka * *********************************************************************** c subroutine rdfrm34(iunit,filehd,nframe, - frmhd,ndjst,msjst,natt,nfalt,nfreq,nstat,icorc,icors) c character*80 infile,filehd character*3 fnum character*28 infohd character*600 dflag character*16 frmhd character*8 tchar character*2 ext1 character*6 ext2 character*1 bcdc character*1 hex(0:15)/'0','1','2','3','4','5','6','7', - '8','9','a','b','c','d','e','f'/ integer nt(16) integer*2 iraw(7056) integer*2 icorc(84,85),icors(84,85) integer icurf/0/ c ***** function bcd -> integer ***** ibcd(bcdc)=ichar(bcdc)/16*10+min(9,mod(ichar(bcdc),16)) c ***** file-type test ***** if (icurf.eq.0) then open(iunit,file=filehd, - access='direct',recl=14136,status='old',iostat=nfstat) if (nfstat.ne.0) then write(6,'('' data-file is 600-frame type'')') else write(6,'('' data-file is event-data type'')') write(6,'('' open file :'',a60)') filehd icurf=1 endif endif c ***** initialize ***** if (nfstat.ne.0) then ifile=(nframe-1)/600+1 iframe=nframe-(ifile-1)*600 if (icurf.ne.ifile) then if (icurf.ne.0) close(iunit) icurf=ifile write(fnum,'(i3.3)') icurf ilfl=index(filehd,' ')-1 infile=filehd(1:ilfl)//fnum open(iunit,file=infile, - access='direct',recl=14136,status='old') read(iunit,rec=1) infohd,dflag write(6,'('' open file :'',a60)') infile end if c nstat=ichar(dflag(iframe:iframe))-48 if (nstat.eq.0) then ccc write(6,'('' nframe ='',i6,'' nstat ='',i6)') nframe,nstat return endif nrec=iframe+1 else nstat=3 nrec=nframe end if c ***** read frame header and data ***** c read(iunit,rec=nrec) nfrm,frmhd,iraw if (nfrm .ne. nframe) then write(6,'('' frame number error : frame ='',i6,'' -> '',i6)') - nframe,nfrm endif tchar=frmhd(1:8) ext1=frmhd(9:10) ext2=frmhd(11:16) nyr=ibcd(tchar(1:1)) nmo=ibcd(tchar(2:2)) ndt=ibcd(tchar(3:3)) nhr=ibcd(tchar(4:4)) nmin=ibcd(tchar(5:5)) nmsec=ibcd(tchar(6:6))*1000 - +ibcd(tchar(7:7))*10+min(9,ichar(tchar(8:8))/16) c write(6,'('' frame :'',i6,'' dd/mm/yy/hh:mm:ss.sss : '', c - i2.2,''/'',i2.2,''/'',i2.2,'' '', c - i2.2,'':'',i2.2,'':'',f6.3)') c - nfrm,nyr,nmo,ndt,nhr,nmin,real(nmsec)/1000. c do 300 i=1,8 nt(i*2-1)=ichar(tchar(i:i))/16 nt(i*2) =mod(ichar(tchar(i:i)),16) 300 continue nbcder=0 do 310 i=1,15 310 if (nt(i) .gt. 9) nbcder=1 if (nbcder.eq.1) then write(6,'('' time stamp error : frame ='',i6,'' '',15a1, - '' -> '',i2.2,i2.2,i2.2,i2.2,i2.2,f6.3)') - nframe,(hex(nt(i)),i=1,15), - nyr,nmo,ndt,nhr,nmin,real(nmsec)/1000. end if c if (nyr .gt. 90) then nyr=nyr+1900 else nyr=nyr+2000 endif ndjst=nyr*10000+nmo*100+ndt msjst=nhr*10000000+nmin*100000+nmsec c * attenuator is always 0 dB at 34 GHz regardless of the status natt=0 nfalt=mod(nt(16)/2,2) if (mod(nt(16),2).eq.1) then nfreq=34 else nfreq=17 end if c ***** create 84x85 arrays ***** call arrng(iraw,icorc,icors) c return end *********************************************************************** * rddat34 - read raw data and create 84x85 array * * for 34 GHz data * * iunit : input,integer fortran file unit number * * for raw data file * * filehd : input,character*80 raw data filename header * * nfrst : input,integer first frame number for integration * * nfrend : input,integer last frame number for integration * * nantst : input,integer(84) antenna status * * ndjst : output,integer observing date in jst (yyyymmdd) * * msjst : output,integer center of observing time * * in jst (hhmmsssss) * * msjsts : output,integer start of observing time * * in jst (hhmmsssss) * * msjste : output,integer end of observing time * * in jst (hhmmsssss) * * natt : output,integer attenuator in db * * nfalt : output,integer frequency alternation * * 0 : fix, 1 : alternate * * nfreq : output,integer observing frequency in GHz * * corc : output,real(84,85) cos uv-map * * cors : output,real(84,85) sin uv-map * * nframe : output,integer number of valid frames * * 0 : all frames are invalid; in this case, other outputs * * are dummy * * June 5, 1998 Y. Hanaoka * *********************************************************************** c subroutine rddat34(iunit,filehd,nfrst,nfrend,nantst, - ndjst,msjst,msjsts,msjste,natt,nfalt,nfreq,nstat, - corc,cors,nframe) c character*80 filehd character*16 frmhd integer nantst(1) integer*2 icorc(84,85),icors(84,85) real corc(84,85),cors(84,85) c do 100 j=1,85 do 100 i=1,84 100 corc(i,j)=0. do 110 j=1,85 do 110 i=1,84 110 cors(i,j)=0. c ***** read data ***** c nframe=0 do 200 iframe=nfrst,nfrend call rdfrm34(iunit,filehd,iframe, - frmhd,ndjst,msjst,natt,nfalt,nfreq,nstat,icorc,icors) if (nstat.eq.0) goto 200 if (nframe.eq.0) then * start time * msjsts=msjst mdjsts=msjst/10000000*3600000+mod(msjst,10000000)/100000*60000 - +mod(msjst,100000) end if * end time * nhr=msjst/10000000 nmin=mod(msjst,10000000)/100000 * 1 sec data : +1000 msec, event data : +50 msec * if (nstat.eq.3) then nmsec=mod(msjst,100000)+50 else nmsec=mod(msjst,100000)+1000 end if if (nmsec.ge.60000) then nmsec=nmsec-60000 nmin=nmin+1 if (nmin.eq.60) then nmin=0 nhr=nhr+1 end if end if msjste=nhr*10000000+nmin*100000+nmsec mdjste=nhr*3600000+nmin*60000+nmsec c nframe=nframe+1 c ***** integration ***** c do 220 j=1,84 do 220 i=1,84 220 corc(i,j)=corc(i,j)+real(icorc(i,j)) do 230 j=1,85 do 230 i=1,84 230 cors(i,j)=cors(i,j)+real(icors(i,j)) c * the values of flux are expressed in unsigned 2-byte integer * do 260 i=1,84 if (icorc(i,85).ge.0) then corc(i,85)=corc(i,85)+real(icorc(i,85)) else corc(i,85)=corc(i,85)+real(icorc(i,85))+65536. end if 260 continue c 200 continue c ***** normalize ***** c if (nframe.eq.0) then write(6,'('' there is no valid frame'')') return endif framen=real(nframe) do 300 j=1,85 do 300 i=1,84 300 corc(i,j)=corc(i,j)/framen do 310 j=1,85 do 310 i=1,84 310 cors(i,j)=cors(i,j)/framen c ***** 0-set for correlations of unuse antennas ***** c do 510 j=1,84 if (nantst(j).eq.1) then do 520 i=1,84 corc(i,j)=0. cors(i,j)=0. corc(j,i)=0. cors(j,i)=0. 520 continue end if 510 continue c ***** observation time ***** c msjst=(mdjsts+mdjste)/2 msjst=msjst/3600000*10000000+mod(msjst,3600000)/60000*100000 - +mod(msjst,60000) c c write(6,'('' '',8i10)') ndjst,msjst,msjsts,msjste, c - natt,nfalt,nfreq,nstat return end *********************************************************************** * snapuv34 - create 2-dimensional uv-map from raw data * * and calibrate * * for 34 GHz data * * iunit : input,integer fortran file unit number * * for raw data file * * filehd : input,character*80 raw data filename header * * nfrst : input,integer first frame number for integration * * nfrend : input,integer last frame number for integration * * nantst : input,integer(84) antenna status * * ical : input,integer calibration menu * * ical=999 : no calibration * * ical=0 : use icald as calibration data * * ical<>0 : calculate calibration data integrating * * abs(ical) frames * * if ical<0, calculated antenna phases are * * shifted, using icald as reference data * * icald : input,integer*2(168) initial calibration data * * : output obtained calibration data * * if integrate frames, cal data for the last frame * * taper : input,real taper at 160d * * ndjst : output,integer observing date in jst (yymmdd) * * msjst : output,integer center of observing time * * in jst (hhmmsssss) * * msjsts : output,integer start of observing time * * in jst (hhmmsssss) * * msjste : output,integer end of observing time * * in jst (hhmmsssss) * * natt : output,integer attenuator in db * * nfalt : output,integer frequency alternation * * 0 : fix, 1 : alternate * * nfreq : output,integer observing frequency in GHz * * nstat : output,integer frame status * * 0 : not recorded * * 1 : calibration data * * 2 : solar observation data * * 3 : event data * * 15 : re-transferred data * * flx : output,real(84) total flux * * suma : output,real(512,512) cos uv-map * * sumb : output,real(512,512) sin uv-map * * nfrmst : output,integer status * * 0 : all frames have been read * * 1 : there is an invalid frame; when nfrmst=1, other output * * is not valid * * November 27, 1998 Y. Hanaoka * *********************************************************************** c subroutine snapuv34(iunit,filehd,nfrst,nfrend, - nantst,ical,icald,taper, - ndjst,msjst,msjsts,msjste,natt,nfalt,nfreq,nstat, - flx,suma,sumb,nfrmst) c character*80 filehd integer nantst(1) integer*2 icald(168),icald0(168) real corc(84,85),cors(84,85),cora(512,512),corb(512,512), - flx(84),suma(512,512),sumb(512,512) real calib(168) c c do 300 j=1,512 do 300 i=1,512 suma(i,j)=0. sumb(i,j)=0. 300 continue do 310 i=1,84 flx(i)=0. 310 continue do 320 i=1,168 320 icald0(i)=icald(i) c do 350 iframe=nfrst,nfrend call rddat34(iunit,filehd,iframe,iframe,nantst, - ndjst,msjst,jsts,jste,natt,nfalt,nfreq,nstat, - corc,cors,nframe) c ***** if there is a no data frame, exit ***** if (nframe.eq.0) then nfrmst=1 go to 900 end if c ***** start and end time of integration ***** if (iframe.eq.nfrst) then msjsts=jsts mdjsts=jsts/10000000*3600000+mod(jsts,10000000)/100000*60000 - +mod(jsts,100000) end if c if (iframe.eq.nfrend) then msjste=jste mdjste=jste/10000000*3600000+mod(jste,10000000)/100000*60000 - +mod(jste,100000) end if c ***** get calibration data and calculate r-l phase shift ***** dcns=0. dcew=0. rlns=0. rlew=0. if (ical.ne.999) then if (ical.ne.0) then call rddcal34(iunit,filehd,iframe, - iabs(ical),nantst,icald,nframe) if (ical.lt.0) then call dcshft(icald0,icald,dcns,dcew) ccc write(6,'('' disk center shift ='',2f10.2)') dcns,dcew end if c do 580 i=1,168 c 580 icald(i)=icaldc(i) end if end if c c ***** create and calibrate 2d uv-map ***** call fi2d(corc,cors,cora,corb) if (ical.ne.999) then do 600 i=1,84 600 calib(i)=real(icald(i))/100. do 610 i=85,168 610 calib(i)=real(icald(i))/10000. call addapa(cora,corb,calib,dcns,dcew) end if c ***** integration ***** do 360 i=1,84 360 flx(i)=flx(i)+corc(i,85) do 370 j=1,512 do 370 i=1,512 suma(i,j)=suma(i,j)+cora(i,j) sumb(i,j)=sumb(i,j)+corb(i,j) 370 continue c 350 continue c c ***** center of the integration time ***** msjst=(mdjsts+mdjste)/2 msjst=msjst/3600000*10000000+mod(msjst,3600000)/60000*100000 - +mod(msjst,60000) c ***** normalize ***** framen=real(nfrend-nfrst+1) do 700 i=1,84 flx(i)=flx(i)/framen 700 continue do 710 j=1,512 do 710 i=1,512 suma(i,j)=suma(i,j)/framen sumb(i,j)=sumb(i,j)/framen 710 continue c call filter1(taper,sumra) call filter1(taper,sumrb) call filter1(taper,sumla) call filter1(taper,sumlb) c nfrmst=0 900 return end c c *********************************************************************** * rddcal34 - read raw heliograph data and calibrate * * for 34 GHz data * * iunit : input,integer fortran file unit number * * for raw data file * * filehd : input,character*80 raw data filename header * * iframe : input,integer frame number for calibration * * itgr : input,integer number of frames for integration * * nantst : input,integer(84) antenna status * * icald : output,integer*2(168) antenna phases and gains * * nframe : output,integer number of frames which is actually * * used * * June 5, 1998 Y. Hanaoka * *********************************************************************** c subroutine rddcal34(iunit,filehd,iframe,itgr,nantst, - icald,nframe) c character*80 filehd integer nantst(84) real corc(84,85),cors(84,85) c integer iposns(28),iposew(57),icv(9) real phans(200),phaew(200),amans(200),amaew(200) integer*2 icald(168) c data iposns/0,1,2,3,4,5,6,7,8,9,10,11, - 12,14,16,20,24,28, - 32,40,48,56, - 64,80,96,112,128,144/ data iposew/-160,-144,-128,-112,-96,-80,-64, - -56,-48,-40,-32, - -28,-24,-20,-16,-14,-12, - -11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1, - 0,1,2,3,4,5,6,7,8,9,10,11, - 12,14,16,20,24,28, - 32,40,48,56, - 64,80,96,112,128,144,160/ data icv/1,2,4,6,8,12,16,24,32/ c nantns=28 nantew=57 ncv=9 c call rddat34(iunit,filehd,iframe-(itgr-1)/2, - iframe+itgr/2,nantst, - ndjst,msjst,jsts,jste,natt,nfalt,nfreq,nstat, - corc,cors,nframe) c write(6,'('' nframe ='',i6)') nframe c c inic=iframe-nfrst inic=0 call phcal(nantns,iposns,ncv,icv,1,corc,cors,inic,phans, - sums,nmax) ccc write(6,'('' ph/ns-cal : error ='',f12.4,'' n ='',i8)') ccc - sums,nmax call phcal(nantew,iposew,ncv,icv,2,corc,cors,inic,phaew, - sums,nmax) ccc write(6,'('' ph/ew-cal : error ='',f12.4,'' n ='',i8)') ccc - sums,nmax call amcal(nantns,iposns,ncv,icv,1,corc,cors,inic,amans, - sums,nmax) ccc write(6,'('' am/ns-cal : error ='',f12.4,'' n ='',i8)') ccc - sums,nmax call amcal(nantew,iposew,ncv,icv,2,corc,cors,inic,amaew, - sums,nmax) ccc write(6,'('' am/ew-cal : error ='',f12.4,'' n ='',i8)') ccc - sums,nmax c icald(1)=0 do 600 i=2,28 600 icald(i)=int(phans(i-1)*100.) do 610 i=1,56 610 icald(i+28)=int(phaew(i)*100.) icald(85)=10000 c do 620 i=2,28 icald(i+84)=int(amans(i-1)*10000.) if (nantst(i).eq.1) icald(i+84)=0 620 continue do 630 i=1,28 icald(i+112)=int(amaew(i)*10000.) if (nantst(i+56).eq.1) icald(i+112)=0 630 continue do 640 i=1,28 icald(i+140)=int(amaew(i+28)*10000.) if (nantst(57-i).eq.1) icald(i+140)=0 640 continue c return end *********************************************************************** * rdpar - read component file * * iunit : input,integer fortran file unit number * * for component file * * infpar : input,character*80 parameter filename * * dec1 : output,real declination at -12h JST (degree) * * dec2 : output,real declination at 12h JST (degree) * * dec3 : output,real declination at 36h JST (degree) * * ha1 : output,real hour angle at -12h JST (second) * * ha2 : output,real hour angle at 12h JST (second) * * ha3 : output,real hour angle at 36h JST (second) * * solr1 : output,real solar radius at -12h JST (arcsecond) * * solr2 : output,real solar radius at 12h JST (arcsecond) * * solr3 : output,real solar radius at 36h JST (arcsecond) * * solp1 : output,real solar polar angle at -12h JST (degree) * * solp2 : output,real solar polar angle at 12h JST (degree) * * solp3 : output,real solar polar angle at 36h JST (degree) * * solb1 : output,real solar b0 at -12h JST (degree) * * solb2 : output,real solar b0 at 12h JST (degree) * * solb3 : output,real solar b0 at 36h JST (degree) * * nantst : output,integer(84) antenna status * * October 13,1995 Y. Hanaoka * *********************************************************************** c subroutine rdpar(iunit,infpar, - dec1,dec2,dec3,ha1,ha2,ha3, - solr1,solr2,solr3,solp1,solp2,solp3,solb1,solb2,solb3, - nantst) c character*80 infpar integer nantst(84) c open(iunit,file=infpar,access='sequential',status='old') read(iunit,'(3f10.5,6f10.2,6f10.4,84i1)') dec1,dec2,dec3, - ha1,ha2,ha3,solr1,solr2,solr3,solp1,solp2,solp3, - solb1,solb2,solb3,(nantst(i),i=1,84) close(iunit) c write(6,'('' antenna status =''/'' '',28i2/'' '',28i2 c - /'' '',28i2))') c - (nantst(i),i=1,84) c return end *********************************************************************** * puthdr - write standard header information to fits file * * iunit : input,integer fortran unit number for fits file * * nbit : input,integer bitpix of fits header * * naxis : input,integer naxis of fits header * * naxes : input,integer() naxis1,naxis2,... of fits header * * ndjst : input,integer observation date in jst (yyyymmdd) * * msjst : input,integer center of observation time * * in jst (hhmmsssss) * * msjsts : input,integer start of observation time * * in jst (hhmmsssss) * * msjste : input,integer end of observation time * * in jst (hhmmsssss) * * nfrst : input,integer first frame number * * nfrend : input,integer last frame number * * npol : input,integer polarization r:1,l:-1,r+l:0,r-l:2, * * l-r:-2,r|l:10,(r-l)/(r+l):12 * * natt : input,integer attenuator 0 or 10 (dB) * * nfalt : input,integer frequency alternation * * 0 : fix, 1 : alternate * * nfreq : input,integer observation frequency * * nstat : input,integer frame status * * 0 : not recorded * * 1 : calibration data * * 2 : solar observation data * * 3 : event data * * 15 : re-transferred data * * dattyp : input,character*18 data type * * * ndut : input,integer observation date in ut (yyyymmdd) * * msut : input,integer center of observation time * * in ut (hhmmsssss) * * msuts : input,integer start of observation time * * in ut (hhmmsssss) * * msute : input,integer end of observation time * * in ut (hhmmsssss) * * June 5, 1998 Y. Hanaoka * *********************************************************************** c subroutine puthdr(iunit,nbit,naxis,naxes, - ndjst,msjst,msjsts,msjste,nfrst,nfrend, - npol,natt,nfalt,nfreq,nstat,dattyp) c integer naxes(2) character*18 dtobs,tmobs,strobs,endobs,jstdt,jsttm,jsts,jste, - pol,att,alt,freq,frstat,dattyp character*44 cmnt c * msut : center of the observing time * msuts : start of the observing time * msute : end of the observing time * ndut : date at the start of the observing time ndut=ndjst msut=msjst call jst2ut(ndut,msut) ndut=ndjst msute=msjste call jst2ut(ndut,msute) ndut=ndjst msuts=msjsts call jst2ut(ndut,msuts) c ***** c cmnt=' ' c ***** write primary header information ***** call ftpprh(iunit,1,nbit,naxis,naxes,0,0,0,istat) c write(dtobs,'(i4.4,''-'',i2.2,''-'',i2.2)') ndut/10000, - mod(ndut,10000)/100,mod(ndut,100) c write(dtobs,'(i2.2,''/'',i2.2,''/'',i2.2)') mod(ndut,100), c - mod(ndut,10000)/100,ndut/10000 write(tmobs,'(i2.2,'':'',i2.2,'':'',i2.2,''.'',i3.3)') - msut/10000000,mod(msut,10000000)/100000, - mod(msut,100000)/1000,mod(msut,1000) write(strobs,'(i2.2,'':'',i2.2,'':'',i2.2,''.'',i3.3)') - msuts/10000000,mod(msuts,10000000)/100000, - mod(msuts,100000)/1000,mod(msuts,1000) write(endobs,'(i2.2,'':'',i2.2,'':'',i2.2,''.'',i3.3)') - msute/10000000,mod(msute,10000000)/100000, - mod(msute,100000)/1000,mod(msute,1000) write(jstdt,'(i4.4,''-'',i2.2,''-'',i2.2)') ndjst/10000, - mod(ndjst,10000)/100,mod(ndjst,100) c write(jstdt,'(i2.2,''/'',i2.2,''/'',i2.2)') mod(ndjst,100), c - mod(ndjst,10000)/100,ndjst/10000 write(jsttm,'(i2.2,'':'',i2.2,'':'',i2.2,''.'',i3.3)') - msjst/10000000,mod(msjst,10000000)/100000, - mod(msjst,100000)/1000,mod(msjst,1000) c write(jsts,'(i2.2,'':'',i2.2,'':'',f6.3)') msjsts/10000000, c - mod(msjsts,10000000)/100000,real(mod(msjsts,100000))/1000. write(jsts,'(i2.2,'':'',i2.2,'':'',i2.2,''.'',i3.3)') - msjsts/10000000,mod(msjsts,10000000)/100000, - mod(msjsts,100000)/1000,mod(msjsts,1000) c write(jste,'(i2.2,'':'',i2.2,'':'',f6.3)') msjste/10000000, c - mod(msjste,10000000)/100000,real(mod(msjste,100000))/1000. write(jste,'(i2.2,'':'',i2.2,'':'',i2.2,''.'',i3.3)') - msjste/10000000,mod(msjste,10000000)/100000, - mod(msjste,100000)/1000,mod(msjste,1000) c if (npol.eq.1) then write(pol,'(''rcp'')') else if (npol.eq.-1) then write(pol,'(''lcp'')') else if (npol.eq.0) then write(pol,'(''r+l'')') else if (npol.eq.2) then write(pol,'(''r-l'')') else if (npol.eq.-2) then write(pol,'(''l-r'')') else if (npol.eq.10) then write(pol,'(''r|l'')') else if (npol.eq.12) then write(pol,'(''(r-l)/(r+l)'')') else write(pol,'(''unknown'')') end if write(att,'(i2.2,''dB'')') natt if (nfalt.eq.0) then write(alt,'(''fix'')') else write(alt,'(''alt'')') endif write(freq,'(i2.2,''GHz'')') nfreq if (nstat.eq.1) then write(frstat,'(''calib'')') else if (nstat.eq.3) then write(frstat,'(''event'')') else write(frstat,'(''1-sec obs'')') endif endif c call ftpkys(iunit,'date-obs',dtobs, cmnt,istat) call ftpkys(iunit,'time-obs',tmobs, cmnt,istat) call ftpkys(iunit,'strt-obs',strobs,cmnt,istat) call ftpkys(iunit,'end-obs',endobs, cmnt,istat) call ftpkys(iunit,'jstdate', jstdt, cmnt,istat) call ftpkys(iunit,'jsttime', jsttm, cmnt,istat) call ftpkys(iunit,'jst-strt',jsts, cmnt,istat) call ftpkys(iunit,'jst-end', jste, cmnt,istat) call ftpkyj(iunit,'startfrm',nfrst, cmnt,istat) call ftpkyj(iunit,'endfrm', nfrend,cmnt,istat) c call ftpkys(iunit,'polariz', pol, cmnt,istat) call ftpkys(iunit,'att-10db',att, cmnt,istat) call ftpkys(iunit,'obs-mode',alt, cmnt,istat) call ftpkys(iunit,'obs-freq',freq, cmnt,istat) call ftpkys(iunit,'frm-stat',frstat,cmnt,istat) call ftpkys(iunit,'data-typ',dattyp,cmnt,istat) c call ftpkys(iunit,'object','sun',cmnt,istat) call ftpkys(iunit,'telescop','radioheliograph',cmnt,istat) call ftpkys(iunit,'origin','nobeyama radio obs',cmnt,istat) c call ftpkys(iunit,'hdrident','HeliogFITS 2.0',cmnt,istat) c if (istat.ne.0) write(6,'('' puthdr : istat ='',i6)') istat return end c c *********************************************************************** * gethdr - read standard header information to fits file * * iunit : input,integer fortran unit number for fits file * * nbit : output,integer bitpix of fits header * * naxis : output,integer naxis of fits header * * naxes : output,integer() naxis1,naxis2,... of fits header * * ndut : output,integer observation date in ut (yyyymmdd) * * msut : output,integer center of observation time * * msuts : output,integer start of observation time * * in ut (hhmmsssss) * * msute : output,integer end of observation time * * in ut (hhmmsssss) * * ndjst : output,integer observation date in jst (yyyymmdd) * * msjst : output,integer center of observation time * * in jst (hhmmsssss) * * msjsts : output,integer start of observation time * * in jst (hhmmsssss) * * msjste : output,integer end of observation time * * in jst (hhmmsssss) * * nfrst : output,integer first frame number * * nfrend : output,integer last frame number * * npol : output,integer polarization r:1,l:-1,r+l:0,r-l:2 * * natt : output,integer attenuator 0 or 10 (dB) * * nfalt : output,integer frequency alternation * * 0 : fix, 1 : alternate * * nfreq : output,integer observation frequency 17 (GHz) * * nstat : output,integer frame status * * 1 : calibration data * * 2 : solar observation data * * 3 : event data * * dattyp : output,character*18 data type * * June 17, 1998 Y. Hanaoka * *********************************************************************** c subroutine gethdr(iunit,nbit,naxis,naxes,ndut,msut,msuts,mutue, - ndjst,msjst,msjsts,msjste,nfrst,nfrend, - npol,natt,nfalt,nfreq,nstat,dattyp) c integer naxes(2) logical simple,extend character*18 dattyp,hdrident character*44 cmnt c ***** read primary header information ***** call ftgprh(iunit,simple,nbit,naxis,naxes,npc,ngc,extend,istat) c call ftgkys(iunit,'hdrident',hdrident, cmnt,istath) if (istath.eq.0) then if (hdrident(1:14) .eq. 'HeliogFITS 2.0') then call gethdr20(iunit,ndut,msut,msuts,mutue, - ndjst,msjst,msjsts,msjste,nfrst,nfrend, - npol,natt,nfalt,nfreq,nstat,dattyp) endif else call gethdr10(iunit,ndut,msut, - ndjst,msjst,msjsts,msjste,nfrst,nfrend, - npol,natt,nfalt,nfreq,dattyp) msuts=-1 msute=-1 nstat=-1 endif c if (istat.ne.0) write(6,'('' gethdr : istat ='',i6)') istat return end c c *********************************************************************** * gethdr20 - read ver2.0 FITS header information * * iunit : input,integer fortran unit number for fits file * * ndut : output,integer observation date in ut (yyyymmdd) * * msut : output,integer center of observation time * * msuts : output,integer start of observation time * * in ut (hhmmsssss) * * msute : output,integer end of observation time * * in ut (hhmmsssss) * * ndjst : output,integer observation date in jst (yyyymmdd) * * msjst : output,integer center of observation time * * in jst (hhmmsssss) * * msjsts : output,integer start of observation time * * in jst (hhmmsssss) * * msjste : output,integer end of observation time * * in jst (hhmmsssss) * * nfrst : output,integer first frame number * * nfrend : output,integer last frame number * * npol : output,integer polarization r:1,l:-1,r+l:0,r-l:2 * * natt : output,integer attenuator 0 or 10 (dB) * * nfalt : output,integer frequency alternation * * 0 : fix, 1 : alternate * * nfreq : output,integer observation frequency 17 (GHz) * * nstat : output,integer frame status * * 1 : calibration data * * 2 : solar observation data * * 3 : event data * * dattyp : output,character*18 data type * * June 5, 1998 Y. Hanaoka * *********************************************************************** c subroutine gethdr20(iunit,ndut,msut,msuts,msute, - ndjst,msjst,msjsts,msjste,nfrst,nfrend, - npol,natt,nfalt,nfreq,nstat,dattyp) c character*18 dtobs,tmobs,strobs,endobs,jstdt,jsttm,jsts,jste, - pol,att,alt,freq,frstat,dattyp character*44 cmnt c call ftgkys(iunit,'date-obs',dtobs, cmnt,istat) if (istat.eq.0) then read(dtobs,'(i4,''-'',i2,''-'',i2)') id1,id2,id3 ndut=id1*10000+id2*100+id3 else istat=0 write(6,'('' input ndut : '')') read(5,*) ndut end if call ftgkys(iunit,'time-obs',tmobs, cmnt,istat) if (istat.eq.0) then read(tmobs,'(i2,'':'',i2,'':'',f6.3)') id1,id2,rd3 msut=id1*10000000+id2*100000+int(rd3*1000.+0.5) else istat=0 write(6,'('' input msut : '')') read(5,*) msut end if call ftgkys(iunit,'strt-obs',strobs, cmnt,istat) if (istat.eq.0) then read(strobs,'(i2,'':'',i2,'':'',f6.3)') id1,id2,rd3 msuts=id1*10000000+id2*100000+int(rd3*1000.+0.5) else istat=0 write(6,'('' input msuts : '')') read(5,*) msuts end if call ftgkys(iunit,'end-obs',endobs, cmnt,istat) if (istat.eq.0) then read(endobs,'(i2,'':'',i2,'':'',f6.3)') id1,id2,rd3 msute=id1*10000000+id2*100000+int(rd3*1000.+0.5) else istat=0 write(6,'('' input msute : '')') read(5,*) msute end if call ftgkys(iunit,'jstdate', jstdt, cmnt,istat) if (istat.eq.0) then if (jstdt(3:3) .eq. '/') then read(jstdt,'(i2,''/'',i2,''/'',i2)') id1,id2,id3 ndjst=(id3+1900)*10000+id2*100+id1 else read(jstdt,'(i4,''-'',i2,''-'',i2)') id1,id2,id3 ndjst=id1*10000+id2*100+id3 endif else istat=0 write(6,'('' input ndjst : '')') read(5,*) ndjst end if call ftgkys(iunit,'jsttime', jsttm, cmnt,istat) if (istat.eq.0) then read(jsttm,'(i2,'':'',i2,'':'',f6.3)') id1,id2,rd3 msjst=id1*10000000+id2*100000+int(rd3*1000.+0.5) else istat=0 write(6,'('' input msjst : '')') read(5,*) msjst end if call ftgkys(iunit,'jst-strt',jsts, cmnt,istat) if (istat.eq.0) then read(jsts,'(i2,'':'',i2,'':'',f6.3)') id1,id2,rd3 msjsts=id1*10000000+id2*100000+int(rd3*1000.+0.5) else istat=0 write(6,'('' input msjsts : '')') read(5,*) msjsts end if call ftgkys(iunit,'jst-end', jste, cmnt,istat) if (istat.eq.0) then read(jste,'(i2,'':'',i2,'':'',f6.3)') id1,id2,rd3 msjste=id1*10000000+id2*100000+int(rd3*1000.+0.5) else istat=0 write(6,'('' input msjste : '')') read(5,*) msjste end if call ftgkyj(iunit,'startfrm',nfrst, cmnt,istat) if (istat.ne.0) then istat=0 write(6,'('' input nfrst : '')') read(5,*) nfrst end if call ftgkyj(iunit,'endfrm', nfrend,cmnt,istat) if (istat.ne.0) then istat=0 write(6,'('' input nfrend : '')') read(5,*) nfrend end if call ftgkys(iunit,'polariz',pol, cmnt,istat) if (istat.eq.0) then if (pol(1:3).eq.'rcp') then npol=1 else if (pol(1:3).eq.'lcp') then npol=-1 else if (pol(1:3).eq.'r+l') then npol=0 else if (pol(1:3).eq.'r-l') then npol=2 else if (pol(1:3).eq.'l-r') then npol=-2 else if (pol(1:3).eq.'r|l') then npol=10 else if (pol(1:11).eq.'(r-l)/(r+l)') then npol=12 else npol=99 end if else istat=0 write(6,'('' input npol : '')') read(5,*) npol end if call ftgkys(iunit,'att-10db',att, cmnt,istat) if (istat.eq.0) then read(att,'(i2)') natt else istat=0 write(6,'('' input natt : '')') read(5,*) natt end if call ftgkys(iunit,'obs-mode',alt, cmnt,istat) if (istat.eq.0) then if (alt(1:3).eq.'fix') then nfalt=0 else if (alt(1:3).eq.'alt') then nfalt=1 endif else istat=0 write(6,'('' input nfalt : '')') read(5,*) nfalt end if call ftgkys(iunit,'obs-freq',freq, cmnt,istat) if (istat.eq.0) then read(freq,'(i2)') nfreq else istat=0 write(6,'('' input nfreq : '')') read(5,*) nfreq end if call ftgkys(iunit,'frm-stat',frstat, cmnt,istat) if (istat.eq.0) then if (frstat(1:9).eq.'1-sec obs') then nstat=2 else if (frstat(1:5).eq.'event') then nfstat=3 else if (frstat(1:5).eq.'calib') then nfstat=1 endif else istat=0 write(6,'('' input nstat : '')') read(5,*) nstat end if call ftgkys(iunit,'data-typ',dattyp,cmnt,istat) if (istat.ne.0) then istat=0 write(6,'('' input dattyp : '')') read(5,'(a)') dattyp end if c if (istat.ne.0) write(6,'('' gethdr : istat ='',i6)') istat return end c c *********************************************************************** * gethdr10 - read ver1.0 FITS header information * * iunit : input,integer fortran unit number for fits file * * ndut : output,integer observation date in ut (yyyymmdd) * * msut : output,integer center of observation time * * ndjst : output,integer observation date in jst (yyyymmdd) * * msjst : output,integer center of observation time * * in jst (hhmmsssss) * * msjsts : output,integer start of observation time * * in jst (hhmmsssss) * * msjste : output,integer end of observation time * * in jst (hhmmsssss) * * nfrst : output,integer first frame number * * nfrend : output,integer last frame number * * npol : output,integer polarization r:1,l:-1,r+l:0,r-l:2 * * natt : output,integer attenuator 0 or 10 (dB) * * nfalt : output,integer frequency alternation * * 0 : fix, 1 : alternate * * nfreq : output,integer observation frequency 17 (GHz) * * dattyp : output,character*18 data type * * June 5, 1998 Y. Hanaoka * *********************************************************************** c subroutine gethdr10(iunit,ndut,msut, - ndjst,msjst,msjsts,msjste,nfrst,nfrend, - npol,natt,nfalt,nfreq,dattyp) c character*18 dtobs,tmobs,jstdt,jsttm,jsts,jste, - pol,att,alt,freq,dattyp character*44 cmnt c call ftgkys(iunit,'date-obs',dtobs, cmnt,istat) if (istat.eq.0) then read(dtobs,'(i2,''/'',i2,''/'',i2)') id1,id2,id3 ndut=(id3+1900)*10000+id2*100+id1 else istat=0 write(6,'('' input ndut : '')') read(5,*) ndut end if call ftgkys(iunit,'time-obs',tmobs, cmnt,istat) if (istat.eq.0) then read(tmobs,'(i2,'':'',i2,'':'',f6.3)') id1,id2,rd3 msut=id1*10000000+id2*100000+int(rd3*1000.+0.5) else istat=0 write(6,'('' input msut : '')') read(5,*) msut end if call ftgkys(iunit,'jstdate', jstdt, cmnt,istat) if (istat.eq.0) then if (jstdt(3:3) .eq. '/') then read(jstdt,'(i2,''/'',i2,''/'',i2)') id1,id2,id3 ndjst=(id3+1900)*10000+id2*100+id1 else read(jstdt,'(i4,''-'',i2,''-'',i2)') id1,id2,id3 ndjst=id1*10000+id2*100+id3 endif else istat=0 write(6,'('' input ndjst : '')') read(5,*) ndjst end if call ftgkys(iunit,'jsttime', jsttm, cmnt,istat) if (istat.eq.0) then read(jsttm,'(i2,'':'',i2,'':'',f6.3)') id1,id2,rd3 msjst=id1*10000000+id2*100000+int(rd3*1000.+0.5) else istat=0 write(6,'('' input msjst : '')') read(5,*) msjst end if call ftgkys(iunit,'jst-strt',jsts, cmnt,istat) if (istat.eq.0) then read(jsts,'(i2,'':'',i2,'':'',f6.3)') id1,id2,rd3 msjsts=id1*10000000+id2*100000+int(rd3*1000.+0.5) else istat=0 write(6,'('' input msjsts : '')') read(5,*) msjsts end if call ftgkys(iunit,'jst-end', jste, cmnt,istat) if (istat.eq.0) then read(jste,'(i2,'':'',i2,'':'',f6.3)') id1,id2,rd3 msjste=id1*10000000+id2*100000+int(rd3*1000.+0.5) else istat=0 write(6,'('' input msjste : '')') read(5,*) msjste end if call ftgkyj(iunit,'startfrm',nfrst, cmnt,istat) if (istat.ne.0) then istat=0 write(6,'('' input nfrst : '')') read(5,*) nfrst end if call ftgkyj(iunit,'endfrm', nfrend,cmnt,istat) if (istat.ne.0) then istat=0 write(6,'('' input nfrend : '')') read(5,*) nfrend end if call ftgkys(iunit,'polariz',pol, cmnt,istat) if (istat.eq.0) then if (pol(1:3).eq.'rcp') then npol=1 else if (pol(1:3).eq.'lcp') then npol=-1 else if (pol(1:3).eq.'r+l') then npol=0 else if (pol(1:3).eq.'r-l') then npol=2 else if (pol(1:3).eq.'l-r') then npol=-2 else if (pol(1:3).eq.'r|l') then npol=10 else if (pol(1:11).eq.'(r-l)/(r+l)') then npol=12 else npol=99 end if else istat=0 write(6,'('' input npol : '')') read(5,*) npol end if call ftgkys(iunit,'att-10db',att, cmnt,istat) if (istat.eq.0) then read(att,'(i2)') natt else istat=0 write(6,'('' input natt : '')') read(5,*) natt end if call ftgkys(iunit,'obs-mode',alt, cmnt,istat) if (istat.eq.0) then if (alt(1:3).eq.'fix') then nfalt=0 else if (alt(1:3).eq.'alt') then nfalt=1 endif else istat=0 write(6,'('' input nfalt : '')') read(5,*) nfalt end if call ftgkys(iunit,'obs-freq',freq, cmnt,istat) if (istat.eq.0) then read(freq,'(i2)') nfreq else istat=0 write(6,'('' input nfreq : '')') read(5,*) nfreq end if call ftgkys(iunit,'data-typ',dattyp,cmnt,istat) if (istat.ne.0) then istat=0 write(6,'('' input dattyp : '')') read(5,'(a)') dattyp end if c if (istat.ne.0) write(6,'('' gethdr10 : istat ='',i6)') istat return end c c *********************************************************************** * putcor - write coordinate parameters to fits file * * iunit : input,integer fortran unit number for fits file * * crval1 : input,real crval1 * * crval2 : input,real crval2 * * crpix1 : input,real crpix1 * * crpix2 : input,real crpix2 * * cdelt1 : input,real cdelt1 * * cdelt2 : input,real cdelt2 * * ictype : input integer * * ctype=1 : ctype1/2='solar-west'/'solar-north' * * June 8, 1998 Y. Hanaoka * *********************************************************************** c subroutine putcor(iunit,crval1,crval2,crpix1,crpix2, - cdelt1,cdelt2,ictype) c character*44 cmnt c cmnt=' ' c call ftpkyf(iunit,'crval1',crval1,2,'arcsec',istat) call ftpkyf(iunit,'crval2',crval2,2,'arcsec',istat) call ftpkyf(iunit,'crpix1',crpix1,2,cmnt,istat) call ftpkyf(iunit,'crpix2',crpix2,2,cmnt,istat) call ftpkyf(iunit,'cdelt1',cdelt1,5,'arcsec',istat) call ftpkyf(iunit,'cdelt2',cdelt2,5,'arcsec',istat) if (ictype .eq. 1) then call ftpkys(iunit,'ctype1','solar-west',cmnt,istat) call ftpkys(iunit,'ctype2','solar-north',cmnt,istat) endif c if (istat.ne.0) write(6,'('' putcor : istat ='',i6)') istat return end c c *********************************************************************** * putast - write antenna status to fits file * * iunit : input,integer fortran unit number for fits file * * nantst : input,integer(84) antenna status 0:use 1:unuse * * ANTST-N1 : n01-n14 * * ANTST-N2 : n15-n28 * * ANTST-E1 : e28-e15 * * ANTST-E2 : e14-e01 * * ANTST-W1 : w01-w14 * * ANTST-W2 : w15-w28 * * April 28, 1993 Y. Hanaoka * *********************************************************************** c subroutine putast(iunit,nantst) c integer nantst(1) character*14 antst c c cmnt=' ' c do 100 i=1,14 100 write(antst(i:i),'(i1)') nantst(i) call ftpkys(iunit,'ANTST-N1',antst, - 'antenna status n01-n14',istat) c do 110 i=1,14 110 write(antst(i:i),'(i1)') nantst(i+14) call ftpkys(iunit,'ANTST-N2',antst, - 'antenna status n15-n28',istat) c do 120 i=1,14 120 write(antst(i:i),'(i1)') nantst(i+28) call ftpkys(iunit,'ANTST-E1',antst, - 'antenna status e28-e15',istat) c do 130 i=1,14 130 write(antst(i:i),'(i1)') nantst(i+42) call ftpkys(iunit,'ANTST-E2',antst, - 'antenna status e14-e01',istat) c do 140 i=1,14 140 write(antst(i:i),'(i1)') nantst(i+56) call ftpkys(iunit,'ANTST-W1',antst, - 'antenna status w01-w14',istat) c do 150 i=1,14 150 write(antst(i:i),'(i1)') nantst(i+70) call ftpkys(iunit,'ANTST-W2',antst, - 'antenna status w15-w28',istat) c if (istat.ne.0) write(6,'('' puthdr : istat ='',i6)') istat return end c c *********************************************************************** * getast - read antenna status from fits file * * iunit : input,integer fortran unit number for fits file * * nantst : output,integer(84) antenna status 0:use 1:unuse * * ANTST-N1 : n01-n14 * * ANTST-N2 : n15-n28 * * ANTST-E1 : e28-e15 * * ANTST-E2 : e14-e01 * * ANTST-W1 : w01-w14 * * ANTST-W2 : w15-w28 * * April 28, 1993 Y. Hanaoka * *********************************************************************** c subroutine getast(iunit,nantst) c integer nantst(1) character*14 antst character*44 cmnt c call ftgkys(iunit,'ANTST-N1',antst,cmnt,istat) if (istat.ne.0) then istat=0 write(6,'('' input antenna status n01-n14 : '')') read(5,'(a)') antst end if do 100 i=1,14 100 read(antst(i:i),'(i1)') nantst(i) c call ftgkys(iunit,'ANTST-N2',antst,cmnt,istat) if (istat.ne.0) then istat=0 write(6,'('' input antenna status n15-n28 : '')') read(5,'(a)') antst end if do 110 i=1,14 110 read(antst(i:i),'(i1)') nantst(i+14) c call ftgkys(iunit,'ANTST-E1',antst,cmnt,istat) if (istat.ne.0) then istat=0 write(6,'('' input antenna status e28-e15 : '')') read(5,'(a)') antst end if do 120 i=1,14 120 read(antst(i:i),'(i1)') nantst(i+28) c call ftgkys(iunit,'ANTST-E2',antst,cmnt,istat) if (istat.ne.0) then istat=0 write(6,'('' input antenna status e14-e01 : '')') read(5,'(a)') antst end if do 130 i=1,14 130 read(antst(i:i),'(i1)') nantst(i+42) c call ftgkys(iunit,'ANTST-W1',antst,cmnt,istat) if (istat.ne.0) then istat=0 write(6,'('' input antenna status w01-w14 : '')') read(5,'(a)') antst end if do 140 i=1,14 140 read(antst(i:i),'(i1)') nantst(i+56) c call ftgkys(iunit,'ANTST-W2',antst,cmnt,istat) if (istat.ne.0) then istat=0 write(6,'('' input antenna status w15-w28 : '')') read(5,'(a)') antst end if do 150 i=1,14 150 read(antst(i:i),'(i1)') nantst(i+70) c if (istat.ne.0) write(6,'('' getast : istat ='',i6)') istat return end c c *********************************************************************** * puteph - write ephemeris parameters to fits file * * iounit : input,integer fortran unit number for fits file * * dec1 : input,real declination at -12h JST (degree) * * dec2 : input,real declination at 12h JST (degree) * * dec3 : input,real declination at 36h JST (degree) * * ha1 : input,real hour angle at -12h JST (second) * * ha2 : input,real hour angle at 12h JST (second) * * ha3 : input,real hour angle at 36h JST (second) * * solr1 : input,real solar radius at -12h JST (arcsecond) * * solr2 : input,real solar radius at 12h JST (arcsecond) * * solr3 : input,real solar radius at 36h JST (arcsecond) * * solp1 : input,real solar polar angle at -12h JST (degree) * * solp2 : input,real solar polar angle at 12h JST (degree) * * solp3 : input,real solar polar angle at 36h JST (degree) * * solb1 : input,real solar b0 at -12h JST (degree) * * solb2 : input,real solar b0 at 12h JST (degree) * * solb3 : input,real solar b0 at 36h JST (degree) * * April 28, 1993 Y. Hanaoka * *********************************************************************** c subroutine puteph(iounit,dec1,dec2,dec3,ha1,ha2,ha3, - solr1,solr2,solr3,solp1,solp2,solp3,solb1,solb2,solb3) c c character*44 cmnt c cmnt=' ' c call ftpkyf(iounit,'DECLI1',dec1,5, - 'declination at -12h JST (degree)',istat) call ftpkyf(iounit,'DECLI2',dec2,5, - 'declination at 12h JST (degree)',istat) call ftpkyf(iounit,'DECLI3',dec3,5, - 'declination at 36h JST (degree)',istat) call ftpkyf(iounit,'HOURA1',ha1,2, - 'hour angle at -12h JST (second)',istat) call ftpkyf(iounit,'HOURA2',ha2,2, - 'hour angle at 12h JST (second)',istat) call ftpkyf(iounit,'HOURA3',ha3,2, - 'hour angle at 36h JST (second)',istat) call ftpkyf(iounit,'SOLR1',solr1,2, - 'solar radius at -12h JST (arcsecond)',istat) call ftpkyf(iounit,'SOLR2',solr2,2, - 'solar radius at 12h JST (arcsecond)',istat) call ftpkyf(iounit,'SOLR3',solr3,2, - 'solar radius at 36h JST (arcsecond)',istat) call ftpkyf(iounit,'SOLP1',solp1,4, - 'solar polar angle at -12h JST (degree)',istat) call ftpkyf(iounit,'SOLP2',solp2,4, - 'solar polar angle at 12h JST (degree)',istat) call ftpkyf(iounit,'SOLP3',solp3,4, - 'solar polar angle at 36h JST (degree)',istat) call ftpkyf(iounit,'SOLB1',solb1,4, - 'solar b0 at -12h JST (degree)',istat) call ftpkyf(iounit,'SOLB2',solb2,4, - 'solar b0 at 12h JST (degree)',istat) call ftpkyf(iounit,'SOLB3',solb3,4, - 'solar b0 at 36h JST (degree)',istat) c if (istat.ne.0) write(6,'('' puteph : istat ='',i6)') istat return end c c *********************************************************************** * geteph - read ephemeris parameters to fits file * * iounit : input,integer fortran unit number for fits file * * dec1 : output,real declination at -24h JST (degree) * * dec2 : output,real declination at 0h JST (degree) * * dec3 : output,real declination at 24h JST (degree) * * ha1 : output,real hour angle at -24h JST (second) * * ha2 : output,real hour angle at 0h JST (second) * * ha3 : output,real hour angle at 24h JST (second) * * solr1 : output,real solar radius at -24h JST (arcsecond) * * solr2 : output,real solar radius at 0h JST (arcsecond) * * solr3 : output,real solar radius at 24h JST (arcsecond) * * solp1 : output,real solar polar angle at -24h JST (degree) * * solp2 : output,real solar polar angle at 0h JST (degree) * * solp3 : output,real solar polar angle at 24h JST (degree) * * solb1 : output,real solar b0 at -24h JST (degree) * * solb2 : output,real solar b0 at 0h JST (degree) * * solb3 : output,real solar b0 at 24h JST (degree) * * April 28, 1993 Y. Hanaoka * *********************************************************************** c subroutine geteph(iounit,dec1,dec2,dec3,ha1,ha2,ha3, - solr1,solr2,solr3,solp1,solp2,solp3,solb1,solb2,solb3) c character*44 cmnt c call ftgkye(iounit,'DECLI1',dec1,cmnt,istat) call ftgkye(iounit,'DECLI2',dec2,cmnt,istat) call ftgkye(iounit,'DECLI3',dec3,cmnt,istat) call ftgkye(iounit,'HOURA1',ha1,cmnt,istat) call ftgkye(iounit,'HOURA2',ha2,cmnt,istat) call ftgkye(iounit,'HOURA3',ha3,cmnt,istat) call ftgkye(iounit,'SOLR1',solr1,cmnt,istat) call ftgkye(iounit,'SOLR2',solr2,cmnt,istat) call ftgkye(iounit,'SOLR3',solr3,cmnt,istat) call ftgkye(iounit,'SOLP1',solp1,cmnt,istat) call ftgkye(iounit,'SOLP2',solp2,cmnt,istat) call ftgkye(iounit,'SOLP3',solp3,cmnt,istat) call ftgkye(iounit,'SOLB1',solb1,cmnt,istat) call ftgkye(iounit,'SOLB2',solb2,cmnt,istat) call ftgkye(iounit,'SOLB3',solb3,cmnt,istat) c if (istat.ne.0) write(6,'('' geteph : istat ='',i6)') istat return end c c *********************************************************************** * putoep - write ephemeris parameters to fits file * * iounit : input,integer fortran unit number for fits file * * solr : input,real solar radius 1.0125*optical(arcsecond) * * solp : input,real solar polar angle (degree) * * solb : input,real solar b0 (degree) * * dec : input,real declination (degree) * * ha : input,real hour angle (second) * * az : input,real azimuth (degree) * * alt : input,real altitude (degree) * * za : input,real zenith angle (degree) * * pmat : input,real(4) projection matrix * * April 3, 1997 Y. Hanaoka * *********************************************************************** c subroutine putoep(iounit,solr,solp,solb,dec,ha,az,alt,za,pmat) c real pmat(4) c call ftpkyf(iounit,'SOLR',solr,2, - 'optical solar radius (arcsecond)',istat) call ftpkyf(iounit,'SOLP',solp,4, - 'solar polar angle (degree)',istat) call ftpkyf(iounit,'SOLB',solb,4, - 'solar b0 (degree)',istat) call ftpkyf(iounit,'DEC',dec,4, - 'declination (degree)',istat) call ftpkyf(iounit,'HOURA',ha,2, - 'hour angle (second)',istat) call ftpkyf(iounit,'azimuth',az,4, - 'azimuth (degree)',istat) call ftpkyf(iounit,'altitude',alt,4, - 'altitude (degree)',istat) call ftpkyf(iounit,'zangle',za,4, - 'zenithangle (degree)',istat) call ftpkyf(iounit,'pmat1',pmat(1),5, - 'projection matrix',istat) call ftpkyf(iounit,'pmat2',pmat(2),5, - 'projection matrix',istat) call ftpkyf(iounit,'pmat3',pmat(3),5, - 'projection matrix',istat) call ftpkyf(iounit,'pmat4',pmat(4),5, - 'projection matrix',istat) c if (istat.ne.0) write(6,'('' putoep : istat ='',i6)') istat return end *********************************************************************** * gtrlph - get RCP-LCP antenna phase differences * * ndjst : input,integer observation date in JST (yyyymmdd) * * rlph : output,real(84) antenna phase differences * * April 19,1999 Y. Hanaoka * *********************************************************************** c subroutine gtrlph(ndjst,rlph) c real rlph(84),rlph1(84),rlph2(84),rlph3(84) c ***** sun0628rlphdif ***** data rlph1/ - 0.00, -4.03, -75.21, -21.17, -21.71, -15.37, 25.35, - -35.85, -30.48, -47.45, -35.35, -5.55, -41.21, -33.28, - -35.76, 8.39, -28.47, -65.54, -83.06, -21.76, -78.84, - -0.05, -48.55,-133.69, -5.46, -39.98, -74.39, -47.20, - -23.58, -26.46, -55.01, -66.11, -21.89, -39.60, 3.58, - -3.46, -40.04, -36.36, -20.18, -63.51, -52.00, -48.86, - -60.03, -21.27, -17.42, -66.35, -16.73, -17.76, -15.70, - -34.83, -33.05, -50.46, -34.46, -23.05, -32.89, -57.34, - 0.00, -25.77, -16.49, 27.45, -35.71, -38.96, 3.78, - -57.58, -32.48, -49.23, -29.52, -43.15, -35.21, -31.82, - -32.89, -9.94, -25.04, -49.20, -84.94, -17.64, -74.44, - 10.00, -65.67, -19.46, -76.04, -53.33, 2.40, -28.47/ c ***** sun0706_2rlphdif ***** data rlph2/ - 0.00, -1.61, -1.27, -3.22, 0.43, -0.88, -3.18, - -3.20, -5.30, -1.75, 12.32, -5.02, -2.62, -4.71, - -2.35, -2.04, -4.23, -5.66, -3.39, -3.79, 4.01, - 0.02, 1.48, 5.60, 0.82, 1.58, 999., -0.28, - -0.58, 2.38, 3.97, 6.05, 1.37, -0.20, -1.94, - -6.73, -2.18, -4.12, -4.87, -0.71, 1.58, 5.40, - 1.99, 8.10, 10.66, 9.29, 3.35, -0.98, -1.76, - 4.80, 5.35, 11.09, 6.86, 8.57, 12.33, 15.91, - 2.28, 3.70, 2.87, 1.53, 4.16, 7.74, 8.87, - 9.64, 5.73, 4.15, 3.87, 2.71, 8.43, 12.43, - 4.50, 1.89, 2.86, -1.53, 6.70, 7.60, 7.89, - 0.54, -1.25, -4.63, 7.65, -5.21, -11.35, -1.32/ c ***** s971101rlphdif ***** data rlph3/ - 0.00, -3.19, -37.25, -12.59, -10.78, -7.12, 10.83, - -19.86, -16.73, -23.88, -11.40, -5.42, -23.18, -18.19, - -21.88, 3.24, -19.30, -36.46, -46.66, -13.21, -40.29, - -2.09, -25.55, -65.66, -3.71, -24.06, -39.07, -27.55, - -14.49, -11.82, -27.66, -31.79, -10.36, -21.36, -1.43, - -6.44, -19.19, -21.85, -12.40, -31.97, -24.24, -22.12, - -30.07, -6.07, -3.68, -28.55, -5.19, -10.82, -6.32, - -14.62, -13.74, -13.38, -7.08, 1.44, -0.52, -0.97, - -1.55, -12.06, -7.54, 14.25, -16.78, -16.88, 5.14, - -26.22, -14.78, -23.97, -12.64, -19.86, -13.17, -10.65, - -13.87, -4.44, -12.67, -27.13, -38.88, -5.86, -34.63, - 4.75, -38.54, -12.16, -23.87, -1.56, -0.50, -2.94/ c if (ndjst .le. 19920629) then do 200 i=1,84 200 rlph(i)=rlph1(i) else if (ndjst .ge. 19971030 .and. ndjst .le. 19971106) then do 210 i=1,84 210 rlph(i)=rlph3(i) else if (ndjst .ge. 19980121 .and. ndjst .le. 19980126) then do 220 i=1,84 220 rlph(i)=rlph3(i) else do 230 i=1,84 230 rlph(i)=rlph2(i) endif c return end *********************************************************************** * phcal : phase calibration subroutine * * assumption : phase s1 = 0, a0 = 0 * * using a part of correlations * * weighted errors in least square method * * double precision * * * * nant : input,integer number of antennas * * ipos : input,integer antenna positions * * dimension of ipos : number of antennas (nant) * * ncv : input,integer number of correlation vector * * icv : input,fourier component vector * * dimension of icv : number of different vector (ncv) * * nsew : input,integer (1)ns or (2)ew * * corc : input,real cos of correlation * * cors : input,real sin of correlation * * inic : input,integer condition initial value * * pha : input,real initial values if inic ne 0 * * pha : output,real antenna phases * * rerr : output,real residual error * * nmax : output,integer number of iterations * * * * phm : phase matrix * * amo : observed correlation amplitudes * * pho : observed correlation phases * * phc : calculated errors of correlation phases * * phx : calculated antenna phases and fourier phases * * vw : work area * * November 16,1992 Y. Hanaoka * *********************************************************************** c c nrow : number of correlations c ncol=nant-1+ncv-1 c subroutine phcal(nant,ipos,ncv,icv,nsew,corc,cors,inic, - pha,rerr,nmax) c implicit real*8 (a-h,o-z) integer ipos(1),icv(1),erst(100),nvco(2000),nac1(2000), - nac2(2000) real*8 phm(2000,200),pho(2000),amo(2000),phc(2000), - phx(200),vw(2200,200),ang(50) real corc(84,85),cors(84,85),pha(1),rerr c external pherr c common /phmatr/nrow,ncol,phm,amo,pho,nac1,nac2 c ***** matrix for calibration ***** c ncol=nant-1+ncv-1 call phmat(nant,ipos,ncv,icv,nrow,nrow1,nrow2,nvco,nac1,nac2,phm) c ***** observed correlation phases ***** c call phobs(corc,cors,nsew,nrow,nac1,nac2,amo,pho) c ***** initial estimation ***** c do 400 i=1,nant 400 if (ipos(i).eq.0) ncent=i c if (inic.eq.0) then c do 420 i=1,nant 420 erst(i)=0 do 440 ic=1,ncol 440 phx(ic)=0. c * initial estimation for east antennas c if (nrow2.gt.nrow1) then iv=0 do 610 ir=nrow1+1,nrow2 do 620 i=nant-1,nant-ncent+1,-1 620 if (phm(ir,i).eq.-1.) go to 630 630 iant2=i do 640 i=iant2,nant-ncent+1,-1 640 if (phm(ir,i).eq.1.) go to 650 i=0 650 iant1=i if (nvco(ir)-1.ne.iv) nsamp=0 iv=nvco(ir)-1 c if (erst(iant2).eq.0) then if (iv.eq.0) then if (iant1.eq.0) then phx(iant2)=-pho(ir) else phx(iant2)=phx(iant1)-pho(ir) end if else if (iant1.eq.0) then phx(iant2)=phx(nant-1+iv)-pho(ir) else phx(iant2)=phx(iant1)+phx(nant-1+iv)-pho(ir) end if end if phx(iant2)=dmod(dmod(phx(iant2)+180.,360.d0)+360.,360.d0)-180. erst(iant2)=1 else nsamp=nsamp+1 if (iant1.eq.0) then ang(nsamp)=pho(ir)+phx(iant2) else ang(nsamp)=pho(ir)-phx(iant1)+phx(iant2) end if c write(6,'('' '',i5,5f10.2)') nsamp,(ang(isamp),isamp=1,nsamp) call dangmn(nsamp,ang,phx(nant-1+iv)) end if 610 continue do 670 i=nant,ncol 670 phx(i)=-phx(i) c c write(6,'('' initial estimation of antenna phases : '' c - '' positions ='',i5,'' -'',i5)') ipos(ncent-1),ipos(1) c write(6,'('' '',10f8.2)') (phx(i),i=nant-ncent+1,nant-1) c end if c * initial estimation for north or west antennas c if (nrow1.ne.0) then c iv=0 do 500 ic=nant,ncol 500 phx(ic)=0. do 510 ir=1,nrow1 do 520 i=nant-ncent,1,-1 520 if (phm(ir,i).eq.-1.) go to 530 530 iant2=i do 540 i=iant2,1,-1 540 if (phm(ir,i).eq.1.) go to 550 i=0 550 iant1=i if (nvco(ir)-1.ne.iv) nsamp=0 iv=nvco(ir)-1 c if (erst(iant2).eq.0) then if (iv.eq.0) then if (iant1.eq.0) then phx(iant2)=-pho(ir) else phx(iant2)=phx(iant1)-pho(ir) end if else if (iant1.eq.0) then phx(iant2)=phx(nant-1+iv)-pho(ir) else phx(iant2)=phx(iant1)+phx(nant-1+iv)-pho(ir) end if end if phx(iant2)=dmod(dmod(phx(iant2)+180.,360.d0)+360.,360.d0)-180. erst(iant2)=1 else nsamp=nsamp+1 if (iant1.eq.0) then ang(nsamp)=pho(ir)+phx(iant2) else ang(nsamp)=pho(ir)-phx(iant1)+phx(iant2) end if call dangmn(nsamp,ang,phx(nant-1+iv)) end if 510 continue c c write(6,'('' initial estimation of antenna phases : '' c - '' positions ='',i5,'' -'',i5)') ipos(ncent+1),ipos(nant) c write(6,'('' '',10f8.2)') (phx(i),i=1,nant-ncent) c end if c c write(6,'('' initial estimation of correlation phases : '' c - '' vectors ='',i5,'' -'',i5)') icv(2),icv(ncv) c write(6,'('' '',10f8.2)') (phx(i),i=nant,ncol) c else do 580 i=1,ncol 580 phx(i)=pha(i) end if c ***** least square calculation ***** c nmax=nrow*ncol*100 epsr=1.d-3 call dnolf1(phx,ncol,pherr,nrow,epsr,nmax, - phc,sums,vw,2200,icon) do 700 i=1,ncol 700 pha(i)=dmod(dmod(phx(i)+180.,360.d0)+360.,360.d0)-180. rerr=sums c if (icon.ne.0) write(6,'('' phcal : icon of dnolf1 ='',i6)') - icon c write(6,'('' icon ='',i6,'' sums ='',f12.4,'' n ='',i8)') c - icon,sums,nmax c write(6,'('' phc =''/'' '',10f8.2)') (phc(i),i=1,nrow) c if (nrow2.gt.nrow1) then c write(6,'('' antenna phases : '' c - '' positions ='',i5,'' -'',i5)') ipos(ncent-1),ipos(1) c write(6,'('' '',10f8.2)') (pha(i),i=nant-ncent+1,nant-1) c end if c if (nrow1.ne.0) then c write(6,'('' antenna phases :'', c - '' positions ='',i5,'' -'',i5)') ipos(ncent+1),ipos(nant) c write(6,'('' '',10f8.2)') (pha(i),i=1,nant-ncent) c end if c write(6,'('' correlation phases : '' c - '' vectors ='',i5,'' -'',i5)') icv(2),icv(ncv) c write(6,'('' '',10f8.2)') (pha(i),i=nant,ncol) c c write(6, c - '(/'' ant1 ant2 amplitude phase error'')') c do 900 ir=1,nrow c 900 write(6,'('' '',2i5,3f13.2,)') c - nac1(ir),nac2(ir),amo(ir),pho(ir), c - phc(ir)/amo(ir)*1000. c return end c *********************************************************************** * pherr - error between observed and calculated phases * * phx : input, real estimated antenna errors and fourier * * phases * * phe : output, real errors between observed and calculated * * phases * * June 22,1992 Y. Hanaoka * *********************************************************************** c subroutine pherr(phx,phe) c implicit real*8 (a-h,o-z) real*8 phm(2000,200),amo(2000),pho(2000),phx(1),phe(1) integer nac1(2000),nac2(2000) common /phmatr/nrow,ncol,phm,amo,pho,nac1,nac2 c do 100 il=1,nrow phe(il)=0. do 110 ir=1,ncol 110 phe(il)=phe(il)+phm(il,ir)*phx(ir) phe(il)=(dmod(dmod(phe(il)-pho(il)+180.,360.d0)+360., - 360.d0)-180.) - *amo(il)/1000. c phe(il)=dlog(dabs((dmod(dmod(phe(il)-pho(il)+180.,360.d0) c - +360.,360.d0)-180.))+1.)*20. c - *amo(il)/1000. c phe(il)=dsin((phe(il)-pho(il))/2.d0/57.2958)*100. c - *amo(il)/1000. 100 continue c return end c *********************************************************************** * phmat - phase matrix * * nant : input, integer number of antennas * * ipos : input, integer antenna positions * * ncv : input, integer number of fourier component * * vectors * * icv : input, integer fourier component vector * * nrow : output, integer number of correlations * * nrow1 : output, integer number of north of east correlations * * nrow2 : output, integer number of west correlations * * nvco : output, integer fourier vector for each correlation * * nac1 : output, integer antenna number 1 for each correlation * * nac2 : output, integer antenna number 2 for each correlation * * phm : output, real phase matrix * * June 22,1992 Y. Hanaoka * *********************************************************************** c subroutine phmat(nant,ipos,ncv,icv, - nrow,nrow1,nrow2,nvco,nac1,nac2,phm) c implicit real*8 (a-h,o-z) integer ipos(nant),icv(ncv),nvco(nrow),nac1(ncv),nac2(ncv) real*8 phm(2000,200) c ncol=nant-1+ncv-1 c do 100 ic=1,ncol do 100 ir=1,2000 100 phm(ir,ic)=0. c do 200 i=1,nant 200 if (ipos(i).eq.0) ncent=i c ir=0 ***** N-N or W-W elements ***** do 300 iv=1,ncv do 310 iant=ncent,nant-1 do 320 jant=iant+1,nant if (ipos(jant)-ipos(iant).ne.icv(iv)) go to 320 ir=ir+1 nac1(ir)=iant-ncent nac2(ir)=jant-ncent if (iant.ne.ncent) phm(ir,iant-ncent)=1. phm(ir,jant-ncent)=-1. nvco(ir)=iv if (iv.ne.1) phm(ir,nant-1+iv-1)=1. go to 310 320 continue 310 continue 300 continue nrow1=ir c ***** E-E elements ***** do 400 iv=1,ncv do 410 iant=ncent,2,-1 do 420 jant=iant-1,1,-1 if (ipos(iant)-ipos(jant).ne.icv(iv)) go to 420 ir=ir+1 nac1(ir)=iant-ncent nac2(ir)=jant-ncent if (iant.ne.ncent) phm(ir,nant-iant)=1. phm(ir,nant-jant)=-1. nvco(ir)=iv if (iv.ne.1) phm(ir,nant-1+iv-1)=-1. go to 410 420 continue 410 continue 400 continue nrow2=ir c ***** E-W elements ***** do 500 iv=2,ncv do 510 iant=ncent-1,1,-1 do 520 jant=ncent+1,nant if (ipos(jant)-ipos(iant).ne.icv(iv)) go to 520 ir=ir+1 nac1(ir)=iant-ncent nac2(ir)=jant-ncent phm(ir,nant-iant)=1. phm(ir,jant-ncent)=-1. nvco(ir)=iv phm(ir,nant-1+iv-1)=1. go to 510 520 continue 510 continue 500 continue c nrow=ir c c write(6,'('' nrow ='',i4,'' nrow1 ='',i4,'' nrow2 ='',i4)') c - nrow,nrow1,nrow2 c do 900 ir=1,nrow c 900 write(6,'('' '',20f4.0)') (phm(ir,ic),ic=1,nant-1+ncv-1) c return end c *********************************************************************** * phobs - observed fourier components of the sun * * inf : input, character filename of correlation data * * nsew : input, integer (1)n-kei (2)e-kei (3)w-kei * * nrow : input, integer number of correlation * * nac1 : input, integer antenna number 1 * * nac2 : input, integer antenna number 1 * * amo : output, real*8 correlation amplitude * * pho : output, real*8 correlation phases * * June 22,1992 Y. Hanaoka * *********************************************************************** c subroutine phobs(corc,cors,nsew,nrow,nac1,nac2,amo,pho) c implicit real*8 (a-h,o-z) real*8 amo(nrow),pho(nrow) real corc(84,85),cors(84,85) integer nac1(nrow),nac2(nrow) c do 100 ir=1,nrow if (nsew.eq.1) then c=corc(nac1(ir)+1,nac2(ir)+1) s=cors(nac1(ir)+1,nac2(ir)+1) else if (nac1(ir).eq.0) then if (nac2(ir).lt.0) then c=corc(1,nac2(ir)+57) s=cors(1,nac2(ir)+57) else c=corc(1,nac2(ir)+56) s=cors(1,nac2(ir)+56) end if else if (nac1(ir).lt.0) then if (nac2(ir).gt.0) then c=corc(nac1(ir)+57,nac2(ir)+56) s=cors(nac1(ir)+57,nac2(ir)+56) else c=corc(nac1(ir)+57,nac2(ir)+57) s=cors(nac1(ir)+57,nac2(ir)+57) end if else if (nac1(ir).gt.0) then c=corc(nac1(ir)+56,nac2(ir)+56) s=cors(nac1(ir)+56,nac2(ir)+56) end if end if amo(ir)=sqrt(c*c+s*s) if (amo(ir).lt.1.e-4) then pho(ir)=0. else pho(ir)=acos(c/amo(ir))*57.2958 end if if (s.lt.0.) pho(ir)=360.-pho(ir) c write(6,'('' '',2i4,4f13.3,)') c - nac1(ir),nac2(ir),c,s,amo(ir),pho(ir) 100 continue c return end c *********************************************************************** * dangmn - mean of angles (real*8) * * nc : input,integer number of sample angles * * ang : input,real*8(nc) angles * * amean : output,real*8 mean of angles * * June 22,1992 Y. Hanaoka * *********************************************************************** c subroutine dangmn(nc,ang,amean) c implicit real*8 (a-h,o-z) real*8 ang(50) c anwid=360. do 100 j=1,nc anmax=-180. anmin=180. do 110 i=1,nc an=dmod(dmod(ang(i)-ang(j)+180.,360.d0)+360.,360.d0)-180. if (an.gt.anmax) anmax=an if (an.lt.anmin) anmin=an 110 continue if (anmax-anmin.lt.anwid) then anwid=anmax-anmin nan=j end if 100 continue c san=0. do 150 i=1,nc 150 san=san+dmod(dmod(ang(i)-ang(nan)+180.,360.d0)+360., - 360.d0)-180. amean=dmod(dmod(san/real(nc)+ang(nan)+180.,360.d0)+360., - 360.d0)-180. c return end *********************************************************************** * amcal : amplitude calibration * * assumption : log(amplitude) a0 = 0 * * using a part of correlations * * weighted errors in least square method * * double precision * * * * nant : input,integer number of antennas * * ipos : input,integer antenna positions * * dimension of ipos : number of antennas (nant) * * ncv : input,integer number of correlation vector * * icv : input,fourier component vector * * dimension of icv : number of different vector (ncv) * * nsew : input,integer (1)ns or (2)ew * * corc : input,real cos of correlation * * cors : input,real sin of correlation * * inic : input,integer condition initial value * * ama : input,real initial values if inic ne 0 * * ama : output,real antenna amplitudes * * rerr : output,real residual error * * nmax : output,integer number of iterations * * * * amm : amplitude matrix * * amo : observed correlation amplitudes * * amlo : log(observed correlation amplitudes) * * amc : calculated errors of correlation amplitudes * * amx : observed amplitude and calculated amplitudes * * vw : work space * * November 16,1992 Y. Hanaoka * *********************************************************************** c c nrow : number of correlations c ncol=nant-1+ncv c subroutine amcal(nant,ipos,ncv,icv,nsew,corc,cors,inic, - ama,rerr,nmax) c implicit real*8 (a-h,o-z) integer ipos(1),icv(1),nvco(2000),nac1(2000), - nac2(2000),nsamp(100) real*8 amm(2000,200),amlo(2000),amo(2000),amc(2000), - amx(2000),vw(2200,200) real corc(84,85),cors(84,85),ama(1),rerr c external amerr c common /ammatr/nrow,ncol,amm,amo,amlo,nac1,nac2 c ***** matrix of calibration ***** c ncol=nant-1+ncv call ammat(nant,ipos,ncv,icv,nrow,nvco,nac1,nac2,amm) c ***** observed correlation amplitudes ***** c call amobs(corc,cors,nsew,nrow,nac1,nac2,amo,amlo) c ***** initial estimation ***** c if (inic.eq.0) then c do 440 ic=1,nant-1 440 amx(ic)=1. do 500 ic=nant,ncol amx(ic)=0. nsamp(ic-nant+1)=0 500 continue c do 520 ir=1,nrow amx(nant-1+nvco(ir))=amx(nant-1+nvco(ir))+amlo(ir) nsamp(nvco(ir))=nsamp(nvco(ir))+1 520 continue do 540 i=1,ncv 540 amx(nant-1+i)=amx(nant-1+i)/nsamp(i) c c write(6,'('' initial estimation of correlation amplitudes : '' c - '' vectors ='',i5,'' -'',i5)') icv(1),icv(ncv) c write(6,'(5('' '',i5,f8.2))') c - (nsamp(i-nant+1),exp(amx(i)),i=nant,ncol) c else do 580 i=1,ncol 580 amx(i)=alog(ama(i)) end if c ***** least square calculation ***** c nmax=nrow*ncol*100 call dnolf1(amx,ncol,amerr,nrow,1.d-6,nmax, - amc,sums,vw,2200,icon) c if (icon.ne.0) write(6,'('' amcal : icon of dnolf1 ='',i6)') - icon c write(6,'('' icon ='',i6,'' sums ='',f12.4,'' n ='',i8)') c - icon,sums,nmax c write(6,'('' amc =''/'' '',10f8.2)') (amc(i),i=1,nrow) c do 600 i=1,nant 600 if (ipos(i).eq.0) ncent=i c do 610 i=1,ncol 610 ama(i)=exp(amx(i)) rerr=sums c c if (ncent.ne.1) then c write(6,'('' antenna amplitudes : '' c - '' positions ='',i5,'' -'',i5)') ipos(ncent-1),ipos(1) c write(6,'('' '',10f8.4)') (ama(i),i=nant-ncent+1,nant-1) c end if c if (ncent.ne.nant) then c write(6,'('' antenna amplitudes :'', c - '' positions ='',i5,'' -'',i5)') ipos(ncent+1),ipos(nant) c write(6,'('' '',10f8.4)') (ama(i),i=1,nant-ncent) c end if c write(6,'('' correlation amplitudes : '' c - '' vectors ='',i5,'' -'',i5)') icv(1),icv(ncv) c write(6,'('' '',10f8.2)') (ama(i),i=nant,ncol) c c write(6, c - '(/'' ant1 ant2 amplitude log(amplitude) error'')') c do 900 ir=1,nrow c 900 write(6,'('' '',2i5,f13.2,3f13.4,)') nac1(ir),nac2(ir), c - amo(ir),amlo(ir),amc(ir)/amo(ir)*1000. c return end c *********************************************************************** * amerr - error between observed and calculated amplitudes * * amx : input, real estimated antenna amplitudes and * * correlation amplitudes * * phe : output, real errors between observed and calculated * * amplitudes * * June 22,1992 Y. Hanaoka * *********************************************************************** c subroutine amerr(amx,ame) c implicit real*8 (a-h,o-z) real*8 amm(2000,200),amo(2000),amlo(2000),amx(1),ame(1) integer nac1(2000),nac2(2000) common /ammatr/nrow,ncol,amm,amo,amlo,nac1,nac2 c do 100 il=1,nrow ame(il)=0. do 110 ir=1,ncol 110 ame(il)=ame(il)+amm(il,ir)*amx(ir) ame(il)=(ame(il)-amlo(il))*amo(il)/1000. 100 continue c return end c *********************************************************************** * ammat - amplitude matrix * * nant : input, integer number of antennas * * ipos : input, integer antenna positions * * ncv : input, integer number of fourier component * * vectors * * icv : input, integer fourier component vector * * nrow : output, integer number of correlations * * nvco : output, integer fourier vector for each correlation * * nac1 : output, integer antenna number 1 for each correlation * * nac2 : output, integer antenna number 2 for each correlation * * amm : output, real amplitude matrix * * June 22,1992 Y. Hanaoka * *********************************************************************** c subroutine ammat(nant,ipos,ncv,icv,nrow,nvco,nac1,nac2,amm) c implicit real*8 (a-h,o-z) integer ipos(nant),icv(ncv),nvco(nrow),nac1(ncv),nac2(ncv) real*8 amm(2000,200) c ncol=nant-1+ncv c do 100 ic=1,ncol do 100 ir=1,2000 100 amm(ir,ic)=0. c do 200 i=1,nant 200 if (ipos(i).eq.0) ncent=i c ***** N-N or W-W elements ***** ir=0 do 300 iv=1,ncv do 310 iant=ncent,nant-1 do 320 jant=iant+1,nant if (ipos(jant)-ipos(iant).ne.icv(iv)) go to 320 ir=ir+1 nac1(ir)=iant-ncent nac2(ir)=jant-ncent if (iant.ne.ncent) amm(ir,iant-ncent)=1. amm(ir,jant-ncent)=1. nvco(ir)=iv amm(ir,nant-1+iv)=1. go to 310 320 continue 310 continue 300 continue c ***** E-E elements ***** do 400 iv=1,ncv do 410 iant=ncent,2,-1 do 420 jant=iant-1,1,-1 if (ipos(iant)-ipos(jant).ne.icv(iv)) go to 420 ir=ir+1 nac1(ir)=iant-ncent nac2(ir)=jant-ncent if (iant.ne.ncent) amm(ir,nant-iant)=1. amm(ir,nant-jant)=1. nvco(ir)=iv amm(ir,nant-1+iv)=1. go to 410 420 continue 410 continue 400 continue c ***** E-W elements ***** do 500 iv=2,ncv do 510 iant=ncent-1,1,-1 do 520 jant=ncent+1,nant if (ipos(jant)-ipos(iant).ne.icv(iv)) go to 520 ir=ir+1 nac1(ir)=iant-ncent nac2(ir)=jant-ncent amm(ir,nant-iant)=1. amm(ir,jant-ncent)=1. nvco(ir)=iv amm(ir,nant-1+iv)=1. go to 510 520 continue 510 continue 500 continue c nrow=ir c c do 900 ir=1,nrow c 900 write(6,'('' '',20f4.0)') (amm(ir,ic),ic=1,nant-1+ncv) c return end c *********************************************************************** * amobs - observed fourier components of the sun * * inf : input, character filename of correlation data * * nsew : input, integer (1)n-kei (2)e-kei (3)w-kei * * nrow : input, integer number of correlation * * nac1 : input, integer antenna number 1 * * nac2 : input, integer antenna number 1 * * amo : output, real*8 correlation amplitude * * amlo : output, real*8 log(correlation phases) * * June 22,1992 Y. Hanaoka * *********************************************************************** c subroutine amobs(corc,cors,nsew,nrow,nac1,nac2,amo,amlo) c implicit real*8 (a-h,o-z) real*8 amo(nrow),amlo(nrow) real corc(84,85),cors(84,85) integer nac1(nrow),nac2(nrow) c do 100 ir=1,nrow if (nsew.eq.1) then c=corc(nac1(ir)+1,nac2(ir)+1) s=cors(nac1(ir)+1,nac2(ir)+1) else if (nac1(ir).eq.0) then if (nac2(ir).lt.0) then c=corc(1,nac2(ir)+57) s=cors(1,nac2(ir)+57) else c=corc(1,nac2(ir)+56) s=cors(1,nac2(ir)+56) end if else if (nac1(ir).lt.0) then if (nac2(ir).gt.0) then c=corc(nac1(ir)+57,nac2(ir)+56) s=cors(nac1(ir)+57,nac2(ir)+56) else c=corc(nac1(ir)+57,nac2(ir)+57) s=cors(nac1(ir)+57,nac2(ir)+57) end if else if (nac1(ir).gt.0) then c=corc(nac1(ir)+56,nac2(ir)+56) s=cors(nac1(ir)+56,nac2(ir)+56) end if end if amo(ir)=sqrt(c*c+s*s) if (amo(ir).ne.0.) amlo(ir)=dlog(amo(ir)) c write(6,'('' '',2i4,4f13.3,)') c - nac1(ir),nac2(ir),c,s,amo(ir),amlo(ir) 100 continue c return end *********************************************************************** * uvobs - extract observable uv-points * * coordinate of the DC component : (257,257) * * nantst : input,integer(84) antenna status * * taper : input,real taper at 160d * * a : input/output,real(512,512) cos of Fourier map * * b : input/output,real(512,512) sin of Fourier map * * June 22,1992 Y. Hanaoka * * Mar 12,1996 koshix * * November 27,1998 Y. Hanaoka * *********************************************************************** c subroutine uvobs(nantst,taper,a,b) c integer ipos(-28:28),nantst(1) real a(512,512),b(512,512),c(512,512) data ipos/-160,-144,-128,-112,-96,-80,-64, - -56,-48,-40,-32, - -28,-24,-20,-16,-14,-12, - -11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1, - 0,1,2,3,4,5,6,7,8,9,10,11, - 12,14,16,20,24,28, - 32,40,48,56, - 64,80,96,112,128,144,160/ c ***** extract observable uv points ***** do 100 l=1,512 do 100 i=1,512 100 c(i,l)=0. c do 200 l=-27,27 do 200 i=-28,28 200 c(ipos(i)+257,ipos(l)+257)=a(ipos(i)+257,ipos(l)+257) c do 210 l=1,512 do 210 i=1,512 210 a(i,l)=c(i,l) c do 220 l=-27,27 do 220 i=-28,28 220 c(ipos(i)+257,ipos(l)+257)=b(ipos(i)+257,ipos(l)+257) c do 230 l=1,512 do 230 i=1,512 230 b(i,l)=c(i,l) c ***** 0-set for correlations of N-kei unuse antennas ***** do 300 iex=1,28 if (nantst(iex).ne.0) then do 400 i=1,512 a(i,ipos(iex-1)+257)=0. b(i,ipos(iex-1)+257)=0. 400 continue do 410 i=1,512 a(i,-ipos(iex-1)+257)=0. b(i,-ipos(iex-1)+257)=0. 410 continue end if 300 continue c ***** 0-set for correlations of E-kei unuse antennas ***** do 310 iex=29,56 if (nantst(iex).ne.0) then do 420 l=258,512 a(ipos(57-iex)+257,l)=0. b(ipos(57-iex)+257,l)=0. 420 continue do 430 l=1,256 a(-ipos(57-iex)+257,l)=0. b(-ipos(57-iex)+257,l)=0. 430 continue end if 310 continue c ***** 0-set for correlations of W-kei unuse antennas ***** do 320 iex=57,84 if (nantst(iex).ne.0) then do 440 l=257,512 a(-ipos(iex-56)+257,l)=0. b(-ipos(iex-56)+257,l)=0. 440 continue do 450 l=1,257 a(ipos(iex-56)+257,l)=0. b(ipos(iex-56)+257,l)=0. 450 continue end if 320 continue c call filter1(taper,a) call filter1(taper,b) c return end *********************************************************************** * uvew1d - uv values of 1-dimension ew * * nsp : input,integer spacing * * msp : input,integer maximum spacing (<256*nsp) * * corc : input,real(84x85) cos of correlations * * cors : input,real(84x85) sin of correlations * * calib : input,real(168) calibration data * * rlew : input,real r-l phase difference(e-w) * * arft : output,real(512) 1-dimensional correlation * * July 16,1993 Y. Hanaoka * *********************************************************************** c subroutine uvew1d(nsp,msp,corc,cors,calib,rlew,arft,nsamp) c integer ipos(57),nsamp(256) real pha(57),ama(57),pho(256),amo(256), - arft(512),ph(50),calib(1) real corc(84,85),cors(84,85) data ipos/-160,-144,-128,-112,-96,-80,-64, - -56,-48,-40,-32, - -28,-24,-20,-16,-14,-12, - -11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1, - 0,1,2,3,4,5,6,7,8,9,10,11, - 12,14,16,20,24,28, - 32,40,48,56, - 64,80,96,112,128,144,160/ c rd=57.2958 c do 100 i=1,28 c pha(i)=calib(85-i) c pha(i+29)=calib(i+28) pha(i)=calib(85-i)+rlew*real(ipos(i)) pha(i+29)=calib(i+28)+rlew*real(ipos(i+29)) ama(i)=calib(85-i+84) ama(i+29)=calib(i+28+84) 100 continue pha(29)=0. ama(29)=1. c write(6,'('' antenna phases :''/(6f10.2))') (pha(i),i=1,57) c write(6,'('' antenna amplitudes :''/(6f10.2))') (ama(i),i=1,57) c c offset=0. c do 110 i=1,57 c 110 offset=offset+alog(ama(i)) c offset=exp(offset/57.) c write(6,'('' amplitude offset ='',f10.4)') offset c ***** data read and calibration ***** c do 200 i=1,512 200 arft(i)=0. do 210 i=1,256 amo(i)=0. pho(i)=0. nsamp(i)=0 210 continue miv=msp/nsp do 300 iv=1,miv nsamp(iv)=0 sam=0. do 310 iant=1,56 do 320 jant=iant+1,57 if (ipos(jant)-ipos(iant).ne.iv*nsp) go to 320 c if (iant.lt.29) then if (jant.lt.29) then c= corc(iant+28,jant+28) s=-cors(iant+28,jant+28) else if (jant.eq.29) then c=corc(iant+28,1) s=-cors(iant+28,1) else c=corc(iant+28,jant+27) s=cors(iant+28,jant+27) end if else if (iant.eq.29) then c=corc(1,jant+27) s=cors(1,jant+27) else if (iant.gt.29) then c=corc(iant+27,jant+27) s=cors(iant+27,jant+27) end if c am=sqrt(c*c+s*s) if (am.lt.1.e-4) then c ph(nsamp(iv))=0. go to 320 else nsamp(iv)=nsamp(iv)+1 ph(nsamp(iv))=acos(c/am)*rd end if if (s.lt.0.) ph(nsamp(iv))=360.-ph(nsamp(iv)) ph(nsamp(iv))=ph(nsamp(iv))-pha(iant)+pha(jant) c sam=sam+am/ama(iant)/ama(jant)*offset*offset sam=sam+am/ama(iant)/ama(jant) c write(6,'(3i4,4f13.3,)') iv,iant-29,jant-29,c,s,am, c - amod(amod(ph(nsamp(iv)),360.)+360.,360) c go to 310 320 continue 310 continue if (nsamp(iv).gt.0) then call angmn(nsamp(iv),ph,pho(iv)) amo(iv)=sam/real(nsamp(iv)) c write(6,'(i5,'' am ='',f10.2,'' ph ='',f10.2)') c - iv,amo(iv),pho(iv) arft(iv*2+1)=amo(iv)*cos(pho(iv)/rd) arft(iv*2+2)=amo(iv)*sin(pho(iv)/rd) end if c c write(6,'('' spacing ='',i6,'' sample ='',i6)') c - iv*nsp,nsamp(iv) 300 continue c arft(1)=corc(1,1) c return end c *********************************************************************** * angmn - mean of angles (real) * * nc : input,integer number of sample angles * * ang : input,real(nc) angles * * amean : output,real mean of angles * * June 30,1992 Y. Hanaoka * *********************************************************************** c subroutine angmn(nc,ang,amean) c real ang(50) c anwid=360. do 100 j=1,nc anmax=-180. anmin=180. do 110 i=1,nc an=amod(amod(ang(i)-ang(j)+180.,360.)+360.,360.)-180. if (an.gt.anmax) anmax=an if (an.lt.anmin) anmin=an 110 continue if (anmax-anmin.lt.anwid) then anwid=anmax-anmin nan=j end if 100 continue c san=0. do 150 i=1,nc 150 san=san+amod(amod(ang(i)-ang(nan)+180.,360.)+360., - 360.)-180. amean=amod(amod(san/real(nc)+ang(nan)+180.,360.)+360., - 360.)-180. c return end *********************************************************************** * caleph - calculate ephemeris parameters at the observing time * * dec1 : input,real declination at -24h JST (degree) * * dec2 : input,real declination at 0h JST (degree) * * dec3 : input,real declination at 24h JST (degree) * * ha1 : input,real hour angle at -24h JST (second) * * ha2 : input,real hour angle at 0h JST (second) * * ha3 : input,real hour angle at 24h JST (second) * * solr1 : input,real solar radius at -24h JST (arcsecond) * * solr2 : input,real solar radius at 0h JST (arcsecond) * * solr3 : input,real solar radius at 24h JST (arcsecond) * * solp1 : input,real solar polar angle at -24h JST (degree) * * solp2 : input,real solar polar angle at 0h JST (degree) * * solp3 : input,real solar polar angle at 24h JST (degree) * * solb1 : input,real solar b0 at -24h JST (degree) * * solb2 : input,real solar b0 at 0h JST (degree) * * solb3 : input,real solar b0 at 24h JST (degree) * * msjst : input,integer observation time in JST (hhmmsssss) * * dec : output,real declination (degree) * * ha : output,real hour angle (second) * * solr : output,real solar radius (arcsec) * * solp : output,real solar polar angle (degree) * * solb : output,real solar b0 (degree) * * az : output,real azimuth (degree) * * alt : output,real altitude (degree) * * za : output,real zenith angle (degree) * * pmat : output,real(4) projection matrix * * April 3,1997 Y. Hanaoka * *********************************************************************** c subroutine caleph(dec1,dec2,dec3,ha1,ha2,ha3, - solr1,solr2,solr3,solp1,solp2,solp3,solb1,solb2,solb3, - msjst,dec,ha,solr,solp,solb,az,alt,za,pmat) c real lat,pmat(4) c ***** second order fitting function ***** fit2nd(y1,y2,y3,x)=(y1+y3-2*y2)/2*x*x+(y3-y1)/2.*x+y2 c rd=57.29577951 pi=3.141592654 c long=138.47698 lat=35.07967 c the inclination of the ground is reduced from geographic latitude c dfrac=(real(msjst/10000000*3600) - +real(mod(msjst,10000000)/100000*60) - +real(mod(msjst,100000))/1000.-43200.)/86400. c dec=fit2nd(dec1,dec2,dec3,dfrac) ha=fit2nd(ha1-86400.,ha2,ha3+86400.,dfrac) solr=fit2nd(solr1,solr2,solr3,dfrac) solp=fit2nd(solp1,solp2,solp3,dfrac) solb=fit2nd(solb1,solb2,solb3,dfrac) c ***** calculate azimuth and altitude ***** call altaz(ha,dec,lat,az,alt,za) c ***** calculate projection matrix ***** ic=1 call promat(lat,az,alt,ic,pmat) c return end c *********************************************************************** * promat - projection matrix between earth and sun * * lat : input,real latitude(degree) * * ic : input,integer ic=1 : sun -> earth projection * * ic=-1 : earth -> sun projection * * pmat : output real projection matrix * * June 22,1992 Y. Hanaoka * *********************************************************************** c subroutine promat(lat,az,alt,ic,pmat) c real lat,pmat(4) c rd=57.29577951 c az1=180.-az c ***** vecter P ***** px=0. py=cos(lat/rd) pz=-sin(lat/rd) c ***** vector Z ***** zx=-sin(az1/rd)*cos(alt/rd) zy=-cos(az1/rd)*cos(alt/rd) zz=sin(alt/rd) c ***** vector X ***** xx=py*zz-pz*zy xy=pz*zx-px*zz xz=px*zy-py*zx xnorm=sqrt(xx*xx+xy*xy+xz*xz) xx=xx/xnorm xy=xy/xnorm xz=xz/xnorm c ****** vector Y ***** yx=zy*xz-zz*xy yy=zz*xx-zx*xz yz=zx*xy-zy*xx c pmat(1)=xx pmat(4)=yy if (ic.eq.1) then pmat(2)=xy pmat(3)=yx else pmat(2)=yx pmat(3)=xy end if c return end c c *********************************************************************** * solaaz - calculation of altitude and azimuth * * ha : input,real hour angle (second) * * dec : input,real declination (degree) * * lat : input,real latitude of observing site (degree) * * az : output,real azimuth (degree) * * alt : output,real altitude (degree) * * za : output,real zenith angle (degree) * * April 3,1997 Y. Hanaoka * *********************************************************************** c subroutine altaz ( ha, dec, lat, az, alt, za ) c real lat c rd=57.29577951 pi=3.141592654 c alt=sin(lat/rd)*sin(dec/rd) - +cos(lat/rd)*cos(dec/rd)*cos(ha/240./rd) if (abs(alt).ge.1.) then alt=sign(90.,alt) else alt=rd*asin(alt) end if c az=(-cos(lat/rd)*sin(dec/rd) - +sin(lat/rd)*cos(dec/rd)*cos(ha/240./rd))/cos(alt/rd) if (abs(az).ge.1.) then az=90.-sign(90.,az) else az=rd*acos(az) end if az=mod(sign(az,sin(ha/240./rd))+360.,360.) c za=(cos(dec/rd)*sin(lat/rd) - -sin(dec/rd)*cos(lat/rd)*cos(ha/240./rd))/cos(alt/rd) if (abs(za).ge.1.) then za=90.-sign(90.,za) else za=rd*acos(za) end if za=mod(sign(za,sin(ha/240./rd))+360.,360.) c return end c c *********************************************************************** * jst2ut - convert jst date/time to ut * * nd : input/output day in yyyymmdd * * mstime : input/output time in miliseconds * * June 5, 1998 Y. Hanaoka * *********************************************************************** c subroutine jst2ut(nd,mstime) c nyr=nd/10000 nmo=mod(nd,10000)/100 ndt=mod(nd,100) c if (nmo.lt.3) then nmo=nmo+12 nyr=nyr-1 end if nday=(1461*nyr)/4-nyr/100+nyr/400 - +(153*nmo-447)/5+ndt mstime=mstime-9*10000000 if (mstime.lt.0) then mstime=mstime+24*10000000 nday=nday-1 end if nday=nday-3 nday=nday+(nday*4+3)/146097-nday/146097 nyr=(nday*4+3)/1461 ndt=mod(nday+nyr-nyr/4,366) nmo=(ndt*5+2)/153 ndt=ndt-((nmo*153+2)/5-1) nmo=nmo+3 if (nmo .gt. 12) then nmo=nmo-12 nyr=nyr+1 end if c nd=nyr*10000+nmo*100+ndt c return end ***************************************************************** * filter1 - low pass filtering * * * * taper : input,real taper at 160d * * uv : input/output,real(512,512) * * Jun 13,1995 koshix * ***************************************************************** subroutine filter1(taper,uv) real taper,sigma,r real uv(512,512) if(taper .ne. 1.) then sigma=160./sqrt(alog(1./taper)) do 200 i=1,512 do 100 j=1,512 r=sqrt(real(i-257)*real(i-257)+real(j-257)*real(j-257)) uv(i,j)=uv(i,j)*exp(-1.*(r*r)/(sigma*sigma)) 100 continue 200 continue endif return end