From 8342c957dc5662ff6fadc084057c17ed21fdfd89 Mon Sep 17 00:00:00 2001 From: Mark Wieringa Date: Thu, 6 Nov 2025 17:16:05 +1100 Subject: [PATCH 1/2] add weather data - first pass --- CMakeLists.txt | 6 +++ inc/prog/fits.h | 24 ++++++----- src/prog/fits.for | 98 ++++++++++++++++++++++++++++++++++++++++---- src/prog/varlist.for | 8 ++-- src/subs/bsrch.for | 46 +++++++++++++++++++++ 5 files changed, 160 insertions(+), 22 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index b075a68f..b05faf50 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,6 +16,9 @@ else() set(MIRIAD_VARIANT "") endif() +# Debug option +option(DEBUG "Enable debug compilation" OFF) + message(STATUS "Building Miriad version ${PROJECT_VERSION}${MIRIAD_VARIANT}") @@ -64,6 +67,9 @@ add_compile_options($<$:-Wno-sometimes-uninitialized>) add_compile_options($<$:-fno-automatic>) add_compile_options($<$:-fallow-argument-mismatch>) add_compile_options($<$:-w>) +if (DEBUG) + add_compile_options($<$:-g>) +endif() # Dependencies find_package(X11 REQUIRED) diff --git a/inc/prog/fits.h b/inc/prog/fits.h index 78d1438d..e3cee211 100644 --- a/inc/prog/fits.h +++ b/inc/prog/fits.h @@ -15,18 +15,21 @@ c----------------------------------------------------------------------- integer uvSTOKES, uvFREQ, uvRA, uvDEC parameter (uvSTOKES=1, uvFREQ=2, uvRA=3, uvDEC=4) - integer MAXCONFG, MAXFREQ, MAXIF, MAXSRC + integer MAXCONFG, MAXFREQ, MAXIF, MAXSRC, MAXTIME parameter (MAXCONFG=40, * MAXSRC=10000, * MAXFREQ=MAXWIN, - * MAXIF=MAXFREQ) + * MAXIF=MAXFREQ, + * MAXTIME=8640) logical emok, inited, jok, llok, mosaic, systok, velcomp integer config, findx(MAXFREQ), freqid, freqids(MAXFREQ), * freqidx, mount(MAXCONFG), nants(MAXCONFG), nchan, * nconfig, nfreq, nif, nsrc, sindx(MAXSRC), srcid, - * srcids(MAXSRC), srcidx, velsys - real dnu, evec, jyperk, systemp, velref + * srcids(MAXSRC), srcidx, velsys, nwth + real dnu, evec, jyperk, systemp, velref, tempc(MAXTIME), + * pressmb(MAXTIME),relhum(MAXTIME),wind(MAXTIME), + * winddir(MAXTIME),wvr(MAXTIME) double precision antpos(3*MAXANT,MAXCONFG), ddec(MAXSRC), * decapp(MAXSRC), decepo(MAXSRC), dra(MAXSRC), * epoch(MAXSRC), eq, freqoff(MAXSRC*MAXIF), @@ -34,14 +37,15 @@ c----------------------------------------------------------------------- * raapp(MAXSRC), raepo(MAXSRC), restfreq(MAXSRC*MAXIF), * sdf(MAXIF*MAXFREQ), sfreq(MAXIF*MAXFREQ), * timeoff(MAXCONFG), timeref, tprev, veldop(MAXSRC), - * height(MAXCONFG) + * height(MAXCONFG),timewth(MAXTIME) character observer*16, source(MAXSRC)*20, telescop*16 common /tables/ raepo, decepo, raapp, decapp, dra, ddec, sfreq, * freqoff, restfreq, veldop, antpos, timeoff, freqref, - * epoch, lat, long, height, tprev, timeref, eq, sdf, dnu, - * evec, systemp, jyperk, velref, nsrc, nif, nchan, nfreq, - * nconfig, nants, srcids, freqids, srcid, freqid, srcidx, - * freqidx, sindx, findx, mount, velsys, config, mosaic, - * velcomp, llok, emok, systok, jok, inited + * epoch, lat, long, height, tprev, timeref, eq, sdf, + * timewth, tempc, pressmb, relhum, wind, winddir, wvr, + * evec, systemp, jyperk, velref, nsrc, nif, nchan, nfreq, + * dnu, nconfig, nants, srcids, freqids, srcid, freqid, + * srcidx, freqidx, sindx, findx, mount, velsys, config, + * mosaic, velcomp, llok, emok, systok, jok, inited, nwth common /tablesc/ observer, source, telescop diff --git a/src/prog/fits.for b/src/prog/fits.for index 592322fa..6233b6b1 100644 --- a/src/prog/fits.for +++ b/src/prog/fits.for @@ -486,6 +486,10 @@ c blfound = .false. if (dobl) call BlInit(lu,nif,blfound) c +c Load any WX tables (weather) +c + call WxLoad(lu,wxfound) +c c Give a summary about various tables. c call tabinfo(lu,blfound) @@ -928,6 +932,64 @@ c end +c*********************************************************************** + + subroutine WxLoad(lu,wxfound) + + integer lu + logical wxfound +c----------------------------------------------------------------------- +c Load any AIPS WX tables. +c----------------------------------------------------------------------- + include 'fits.h' + integer nval + character type*1,units*16 + real dewc(MAXTIME),b,c,t0, prfac + parameter (b=17.625, c=243.04) +c----------------------------------------------------------------------- +c Look for WX tables. + ntab = 0 + call ftabLoc(lu,'AIPS WX',wxfound) + if (wxfound) then + nval = 0 + call ftabInfo(lu,'TIME',type,units,nwth,nval) + if (nwth.gt.MAXTIME) call bug('f','Too many times in WX table') + if (nval.ne.1 .or. type.ne.'D') + * call bug('f','Something screwy with WX table') + call output(' Using weather table information') + call ftabGetd(lu,'TIME',0,timewth) +c call ftabGetr(lu,'TIME_INTERVAL',0,interval) + +c call ftabGeti(lu,'ANTENNA_NO',0,antno) +c call ftabGeti(lu,'SUBARRAY',0,subarr) + call ftabGetr(lu,'TEMPERATURE',0,tempc) + call ftabGetr(lu,'PRESSURE',0,pressmb) + call ftabGetr(lu,'DEWPOINT',0,dewc) + call ftabGetr(lu,'WIND_VELOCITY',0,wind) + call ftabGetr(lu,'WIND_DIRECTION',0,winddir) + call ftabGetr(lu,'WVR_H2O',0,wvr) +c call ftabGetr(lu,'IONOS_ELECTRON',0,elec) +c +c Convert dew point to relative humidity +c Convert sea level pressure to actual +c Convert wind in m/s to km/h +c Convert times to JD, t0 is start of ut day +c + prfac = 1.0 + if (telescop.eq.'atca') prfac = 0.975 + t0 = int(timeref - 0.5) + 0.5 + do i = 1, nwth + dewc(i) = max(-50.0,dewc(i)) + relhum(i) = 100 * exp(b*dewc(i)/(c+dewc(i)))/ + * exp(b*tempc(i)/(c+tempc(i))) + pressmb(i) = prfac * pressmb(i) + wind(i) = wind(i) * 3.6 + timewth(i) = timewth(i) + t0 + enddo + call output(' Saving weather data') + endif + + end c*********************************************************************** subroutine BlInit(lu,nif1,blfound) @@ -1401,7 +1463,7 @@ c Give information on the tables in the file. c c----------------------------------------------------------------------- integer NTAB - parameter (NTAB=6) + parameter (NTAB=7) character tabs(NTAB)*8 logical found,givecal character ename*16 @@ -1410,7 +1472,8 @@ c Externals. integer binsrcha data tabs/'AIPS AN ','AIPS CH ','AIPS FG ', - * 'AIPS FQ ','AIPS OB ','AIPS SU '/ + * 'AIPS FQ ','AIPS OB ','AIPS SU ', + * 'AIPS WX '/ c----------------------------------------------------------------------- call ftabloc(lu,' ',found) if (.not.found) call bug('f','Something is screwy') @@ -1438,8 +1501,6 @@ c call output(' ... it is assumed FILLM applied these.') else if (ename.eq.'AIPS PO') then call output(' Ingoring AIPS planetary ephemeris table') - else if (ename.eq.'AIPS WX') then - call output(' Ignoring AIPS weather table') else if (ename.eq.'AIPS TY') then call output(' Ignoring AIPS system flux cal table values') call output(' ... it is assumed FILLM applied these.') @@ -2150,7 +2211,7 @@ c s = itoaf(nconfig) ls = len1(s) write(line,'(a,a,a,f6.2,a)') - * ' Decrementing times for configration ',s(1:ls),' by', + * ' Decrementing times for configuration ',s(1:ls),' by', * -ut1utc,' seconds (UTC-UT1).' call output(line) endif @@ -2474,6 +2535,8 @@ c config = -1 Tprev = -1 inited = .true. + wthid = 0 + T0 = -1 end @@ -2498,18 +2561,20 @@ c the value corrected for clock differences. c----------------------------------------------------------------------- integer LSRRADIO parameter (LSRRADIO=257) + double precision mjd0 + parameter(mjd0= 2400000.5d0) include 'fits.h' include 'mirconst.h' - integer i,j,k,l + integer i,j,k,l,wid logical newsrc,newfreq,newconfg,newlst,newchi,newvel,neweq - logical neweph + logical neweph, newwth real chi,chi2,dT double precision lst,vel,az,el double precision sfreq0(MAXIF),sdf0(MAXIF),rfreq0(MAXIF) c c Externals. c - integer binsrchi,len1 + integer binsrchi,len1, binsrchd double precision eqeq c if (.not.inited) call TabInit(tno) @@ -2535,6 +2600,23 @@ c newconfg = config.ne.confg config = confg if (config.gt.nconfig) config = 1 + +c get latest weather + if (time.gt.tprev) then + wid = binsrchd(time,timewth,nwth) + newwth = wthid.ne.wid + if (newwth) then + wthid = wid + call uvputvrr(tno,'airtemp',tempc(wthid),1) + call uvputvrr(tno,'pressmb',pressmb(wthid),1) + call uvputvrr(tno,'relhumid',relhum(wthid),1) + call uvputvrr(tno,'wind',wind(wthid),1) + call uvputvrr(tno,'winddir',winddir(wthid),1) +c call uvputvrr(tno,'precipmm',rainmm(wthid),1) +c call uvputvrr(tno,'smonrms',real(wvr(wthid)),1) + endif + endif + c c Correct the time c diff --git a/src/prog/varlist.for b/src/prog/varlist.for index d66b6c5a..d88e8f6b 100644 --- a/src/prog/varlist.for +++ b/src/prog/varlist.for @@ -2,9 +2,9 @@ c*********************************************************************** program varlist implicit none C -C List all variable names, types, and lengths in a u,v data set +C List all variable names, types, and lengths in a u,v data set C -C user inputs: dataset - name of u,v data set +C user inputs: dataset - name of u,v data set C outfile - output file for listing; default = logfile C c= varlist - List all variables in dataset @@ -27,7 +27,7 @@ c rjs 7nov89 Some standardising and cosmetic changes. c lgm 12nov89 Fix so printed variable lengths are non-zero c pjt 30jun93 Wow, 2.5 years of bugfree riding, but now added MAXCHAN c rjs 16sep93 Call logclose. -c rjs 27apr95 Distinguish between zero-length and unset +c rjs 27apr95 Distinguish between zero-length and unset c vjm 29mar12 Handle longer dataset names c ToDo c * fix questionable practice to find all uv vars (at most 300 now) @@ -89,7 +89,7 @@ C if(option(1:3) .eq. 'nam') then do iv=1,nvar,5 jv = min(nvar,iv+4) - write(line,'(10x,5(a8,3x))') + write(line,'(10x,5(a8,3x))') * (var(ivar)(3:10),ivar=iv,jv) call LogWrite(line,more) enddo diff --git a/src/subs/bsrch.for b/src/subs/bsrch.for index 63a48df0..57204263 100644 --- a/src/subs/bsrch.for +++ b/src/subs/bsrch.for @@ -88,3 +88,49 @@ c endif enddo end +c************************************************************************ +c*binsrchd -- Search a sorted list of doubles for insertion place. +c:search,binary-search +c+ + integer function binsrchd(key,keys,nkeys) +c + implicit none + integer nkeys + double precision key,keys(nkeys) +c +c Search for interval where given key would be inserted in a list. Return the +c index of the lower bound. If key less than all, return 0. +c A binary search is used. +c +c Input: +c key The double to search for. +c keys A list of doubles. These are assumed to be in +c order to allow a binary search. +c nkeys The number of integers. +c +c Output: +c binsrchd Either the index (if interval found in the +c list), or zero. +c-- +c------------------------------------------------------------------------ + integer j,k,l +c + k = 1 + l = nkeys + binsrchd = 0 + do while(k.le.l) + j = (k+l)/2 + if(key.lt.keys(j))then + l = j - 1 + else if(key.gt.keys(j))then + k = j + 1 + else + binsrchd = j + return + endif + enddo +c If not found, l will be index of last smaller element or zero + binsrchd = l + return + + end From cc20f98e404c7b52d9aaeb36954817b05b35f5b3 Mon Sep 17 00:00:00 2001 From: Mark Wieringa Date: Mon, 10 Nov 2025 13:20:35 +1100 Subject: [PATCH 2/2] review updates --- src/prog/fits.for | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/prog/fits.for b/src/prog/fits.for index 6233b6b1..04a48e1f 100644 --- a/src/prog/fits.for +++ b/src/prog/fits.for @@ -971,7 +971,7 @@ c call ftabGeti(lu,'SUBARRAY',0,subarr) c call ftabGetr(lu,'IONOS_ELECTRON',0,elec) c c Convert dew point to relative humidity -c Convert sea level pressure to actual +c Convert sea level pressure to actual surface pressure c Convert wind in m/s to km/h c Convert times to JD, t0 is start of ut day c @@ -980,8 +980,8 @@ c t0 = int(timeref - 0.5) + 0.5 do i = 1, nwth dewc(i) = max(-50.0,dewc(i)) - relhum(i) = 100 * exp(b*dewc(i)/(c+dewc(i)))/ - * exp(b*tempc(i)/(c+tempc(i))) + relhum(i) = 100 * exp(b*dewc(i)/(c+dewc(i))- + * b*tempc(i)/(c+tempc(i))) pressmb(i) = prfac * pressmb(i) wind(i) = wind(i) * 3.6 timewth(i) = timewth(i) + t0 @@ -2613,7 +2613,6 @@ c get latest weather call uvputvrr(tno,'wind',wind(wthid),1) call uvputvrr(tno,'winddir',winddir(wthid),1) c call uvputvrr(tno,'precipmm',rainmm(wthid),1) -c call uvputvrr(tno,'smonrms',real(wvr(wthid)),1) endif endif