program obsweighter3 c c This program determines weights for the orbit program c using the new scheme involving separation, method (mcode), c and observer (ocode) c [ orbits.3* --> orbits.4* ] c c Modifications made 2/18/2004: c (1) set maximum value for wt1 for each method c (2) automatically reduce wt3 for HJ measures to 0.1 c (3) reduce wt3 for uncertain/estimated measures or measures with theta c determined from quadrant (i.e., "U" or "L" codes) from 0.5 to 0.1 c c Modifications made 3/22/2004: c (1) use new set of observer weights (observer.weight.2004) c (2) remove automatic HJ reduction, as is now incorporated in observer.weight.2004 c (3) use four "generic" visual weights (VI1 - VI4) based on date c c Modifications made 7/8/2011: c (1) added new method codes (EYE, PTG) to weighting c character*3 ref1,method,mcode(500),ocode1(500),ocode2(500) character*8 ref character*1 c77,c78,c79 character*20 hcode character*9 adate character*6 atheta character*7 arho character*100 title dimension w2(500) c open(10,file='observer.weight.2004',status='UNKNOWN') open(11,file='orbits.3',status='UNKNOWN') open(12,file='orbits.4',status='UNKNOWN') c c read in observer weights c c nobs=169 c do 100 n=1,nobs c 100 read(10,901) mcode(n),ocode1(n),ocode2(n),w2(n) c 901 format(t3,a3,2x,a3,2x,a3,t61,f9.2) 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) date,theta,rho,method,wt, $ hcode,enn,ref1,ref,tel,c77,c78,c79 908 format(6x,f9.4,4x,f6.2,6x,f11.6,7x,a3,f5.1,5x,a20,t63, $ f2.0,t66,a3,t66,a8,f3.0,3a1) if (date .eq. 0.0) go to 600 c c throw out erroneous or unresolved measures c if ((c78 .eq. 'D') .or. (c79 .eq. 'D')) go to 400 if ((c78 .eq. 'E') .or. (c79 .eq. 'E')) go to 400 if ((c78 .eq. 'F') .or. (c79 .eq. 'F')) go to 400 if ((c78 .eq. 'H') .or. (c79 .eq. 'H')) go to 400 if ((c78 .eq. 'J') .or. (c79 .eq. 'J')) go to 400 if ((c78 .eq. 'S') .or. (c79 .eq. 'S')) go to 400 if ((c78 .eq. 'T') .or. (c79 .eq. 'T')) go to 400 c c adjust rho for mas, arcmin, or degree measures c if ((c78 .eq. 'A') .or. (c79 .eq. 'A')) rho=rho/1000. if ((c78 .eq. 'C') .or. (c79 .eq. 'C')) rho=rho*60. if ((c78 .eq. 'H') .or. (c79 .eq. 'H')) rho=rho*3600. c c adjust telescope aperture c if (ref1 .eq. 'AnJ') tel=240. if (ref1 .eq. 'Mrr') tel=240. if (ref1 .eq. 'MKT') tel=23.6*39.37 if (ref1 .eq. 'MkT') tel=23.6*39.37 if (ref1 .eq. 'NOI') tel=37.5*39.37 if (ref .eq. 'Lab1974 ') tel= 200. if (ref .eq. 'Koe1983 ') tel=20.0*39.37 if (tel .eq. 0.0) tel=3.0 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 ((c78 .eq. 'U') .or. (c79 .eq. 'U')) wt3=0.5 if ((c78 .eq. 'L') .or. (c79 .eq. 'L')) 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 rholim=5.3596/tel x=rho/rholim c imethod=3 if ((method .eq. 'VPH') .or. (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. 'VPH') .and. (iobs .eq. 0)) iobs=328 if ((method .eq. 'PTG') .and. (iobs .eq. 0)) iobs=328 if ((method .eq. 'XVI') .and. (iobs .eq. 0)) iobs=394 if ((method .eq. 'EYE') .and. (iobs .eq. 0)) iobs=394 if ((method .eq. 'XSP') .and. (iobs .eq. 0)) iobs=389 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) date,theta,rho,method,wt,hcode 909 format(6x,f9.4,4x,f6.2,7x,f10.6,7x,a3,f5.1,5x,a20) go to 400 c 600 write(12,910) 910 format(95x,' ') go to 200 999 stop 'obsweighter program finished' end