program noorbweight c c This program determines weights for the orbit program using the new c scheme involving separation, method (mcode), and observer (ocode) c [ orbits.3 --> orbits.4 ] c c 040218 : set maximum value for wt1 for each method c : automatically reduce wt3 for HJ measures to 0.1 c : reduce wt3 for uncertain/estimated measures or measures with c theta determined from quadrant (i.e., "U" or "L" codes) from c 0.5 to 0.1 c 040322 : use new set of observer weights (observer.weight.2004) c : remove automatic HJ reduction, as is now incorporated in c observer.weight.2004 c : use four "generic" visual weights (VI1 - VI4) based on date c 110708 : added new method codes (EYE, PTG) to weighting c 120709 : modified for new WDS format c character*3 ref1,method,mcode(500),ocode1(500),ocode2(500) character*8 ref character*1 tcode,rcode,c95,c96,tecode,fcode character*9 adate character*6 atheta character*7 arho character*100 title,data1,data2 dimension w2(500) c open(10,file='observer.weight.2004',status='UNKNOWN') open(11,file='noorb.3',status='UNKNOWN') open(12,file='noorb.4',status='UNKNOWN') c c read in observer weights c nobs=394 do 100 n=1,nobs 100 read(10,901) mcode(n),ocode1(n),ocode2(n),w2(n) 901 format(t7,a3,2x,a3,2x,a3,t99,f9.2) c c read and write wds title and orbital elements lines c 200 read(11,907,end=999,err=999) title 907 format(a100) write(12,907) title read(11,907) title write(12,907) title c c read data and determine weight c 400 read(11,908,end=999,err=999) data1,wt,data2,date,tcode,theta, $ rcode,rho,method,filt,fcode,tel,tecode,enn,ref1, $ ref,c95,c96 908 format(a52,f5.1,a40,t7,f10.5,t19,a1,f7.3,t32,a1,f9.5,t50,a3, $ t63,f4.0,t71,a1,t73,f5.1,a1,f3.0,1x,a3,t83,a8,t95,2a1) if (date .eq. 0.0) go to 600 c c throw out erroneous or unresolved measures c if ((c95 .eq. 'X') .or. (c96 .eq. 'X')) go to 400 if ((c95 .eq. 'F') .or. (c96 .eq. 'F')) go to 400 c c adjust rho for mas, arcmin, or degree measures c if (rcode .eq. 'D') rho=rho*3600. if (rcode .eq. 'M') rho=rho*60. if (rcode .eq. 'm') rho=rho/1000. c c adjust filter wavelength c if (filt .eq. 0.0) filt=550. if (fcode .eq. 'B') filt=440. if (fcode .eq. 'R') filt=700. if (fcode .eq. 'K') filt=2200. if (fcode .eq. 'n') filt=550. if (fcode .eq. 'x') filt=1. if (fcode .eq. 'u') filt=filt*1.0e3 if (fcode .eq. 'm') filt=filt*1.0e6 if (fcode .eq. 'c') filt=filt*1.0e7 if (fcode .eq. 'M') filt=filt*1.0e9 c c adjust telescope aperture, etc c if (tel .eq. 0.0) tel=0.1 if (tecode .eq. 'k') tel=tel*1.0e3 if (tecode .eq. 'M') tel=tel*1.0e6 c if (enn .eq. 0.0) enn=1.0 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c determine weighting factor = sqrt(n) c wt4=sqrt(enn) c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c give half weight to measures flagged as uncertain or estimated c wt3=1.0 if ((tcode .eq. ':') .or. (rcode .eq. 'L')) wt3=0.5 if (rcode .eq. ':') wt3=0.5 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c determine separation (x) in units of Rayleigh limit c use fits to (rms relative rho) vs (rho/rholim) for c appropriate method to determine weighting factor wt1 c c rholim = 1.2 * lambda/D c for rholim in radians and lambda and D in same units c rholim = lambda/(4040*D) c for lambda in nm, D in m, rholim in arcsec c rholim=filt/(4040.*tel) x=rho/rholim c imethod=3 if (method .eq. 'PTG') imethod=2 if (method .eq. 'VIS') imethod=1 if (method .eq. 'HIP') imethod=3 if (method .eq. 'TYC') imethod=3 c if ((imethod .eq. 1) .and. (x .lt. 1.4)) $ w1 = 0.710442 -0.807388*x +0.3711700*x**2 -0.05581570*x**3 if ((imethod .eq. 1) .and. (x .ge. 1.4)) $ w1 = 0.190430 -0.028076*x +0.0021362*x**2 -0.00005746*x**3 if ((imethod .eq. 1) .and. (x .ge. 7.0)) $ w1 = 0.103455 -0.003657*x +0.0000500*x**2 -0.00000022*x**3 if ((imethod .eq. 1) .and. (x .ge. 20.0)) $ w1=0.0485 c if ((imethod .eq. 2) .and. (x .lt. 20.)) $ w1 = 0.110840 -0.011873*x +0.0005499*x**2 -0.00000980*x**3 if ((imethod .eq. 2) .and. (x .ge. 20.)) $ w1 = 0.018280 -0.000184*x +0.0000008*x**2 if ((imethod .eq. 2) .and. (x .ge. 100.)) $ w1 = 0.007500 c if ((imethod .eq. 3) .and. (x .lt. 3.)) $ w1 = 0.397269 -0.298306*x +0.0908508*x**2 -0.00958455*x**3 if ((imethod .eq. 3) .and. (x .ge. 3.)) $ w1 = 0.111693 -0.022882*x +0.0018639*x**2 -0.00004715*x**3 if ((imethod .eq. 3) .and. (x .ge. 10.)) $ w1 = 0.022200 -0.000020*x c wt1=0.01/w1/w1 wt1a=wt1 if ((imethod .eq. 2) .and. (wt1 .gt. 50.)) wt1=50. if ((imethod .eq. 3) .and. (wt1 .gt. 20.)) wt1=20. c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c determine proper observer code and corresponding weight factor c based on "observer quality" c iobs=0 do 500 n=1,nobs if ((method .eq. mcode(n)) .and. (ref1 .eq. ocode1(n))) iobs=n if ((method .eq. mcode(n)) .and. (ref1 .eq. ocode2(n))) iobs=n 500 continue if ((method .eq. 'VIS') .and. (iobs .eq. 0) .and. $ (date .le. 1850.)) iobs=270 if ((method .eq. 'VIS') .and. (iobs .eq. 0) .and. $ (date .gt. 1850.)) iobs=271 if ((method .eq. 'VIS') .and. (iobs .eq. 0) .and. $ (date .gt. 1900.)) iobs=272 if ((method .eq. 'VIS') .and. (iobs .eq. 0) .and. $ (date .gt. 1950.)) iobs=273 if ((method .eq. 'PTG') .and. (iobs .eq. 0)) iobs=328 if ((method .eq. 'EYE') .and. (iobs .eq. 0)) iobs=394 if ((method .eq. 'SPE') .and. (iobs .eq. 0)) iobs=389 if ((method .eq. 'WSI') .and. (iobs .eq. 0)) iobs=389 if ((method .eq. 'HIP') .and. (iobs .eq. 0)) iobs=7 if ((method .eq. 'TYC') .and. (iobs .eq. 0)) iobs=8 if ((method .eq. 'CSP') .and. (iobs .eq. 0)) iobs=1 c c temporarily set observer weights = 1.0 for redetermining these values c wt2=w2(iobs) c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c combine all factors to determine overall weight c wt=wt1*wt2*wt3*wt4 if ((wt .lt. 0.1) .and. (wt .ne. 0.0)) wt= 0.1 if (wt .gt. 99.9) wt=99.9 write(12,909) data1,wt,data2 909 format(a52,f5.1,a40) go to 400 c 600 write(12,910) 910 format(95x,' ') go to 200 999 stop 'obsweighter program finished' end