program obsweighter2
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
	character*3 ref1,method,mcode(200),ocode1(200),ocode2(200)
	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(200)
c
	open(10,file='observer.weight',status='UNKNOWN')
	open(11,file='orbits.3',status='UNKNOWN')
	open(12,file='orbits.4',status='UNKNOWN')
c
c	read in observer weights 
c
	nobs=169
	do 100 n=1,nobs
  100	read(10,901) mcode(n),ocode1(n),ocode2(n),w2(n)
  901	format(t3,a3,2x,a3,2x,a3,t61,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 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	adjust for uncertain measures
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	adjust rho for minutes or milliarcseconds codes
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
	if ((c78 .eq. 'C') .or. (c79 .eq. 'C')) rho=rho*60.
	if ((c78 .eq. 'A') .or. (c79 .eq. 'A')) rho=rho/1000.
c
	rholim=5.3596/tel
	x=rho/rholim
c
	imethod=3
	if (method .eq. 'VPH') imethod=2
	if (method .eq. 'VIS') imethod=1
	if (ref1   .eq. 'HIP') imethod=3
	if (ref1   .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
	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))       iobs=1
	if ((method .eq. 'VPH'     .and. iobs .eq. 0))       iobs=2
	if ((method .eq. 'XVI'     .and. iobs .eq. 0))       iobs=3
	if ((method .eq. 'XSP'     .and. iobs .eq. 0))       iobs=4
	if ((method .eq. 'HIP'     .and. iobs .eq. 0))       iobs=148
	if ((method .eq. 'TYC'     .and. iobs .eq. 0))       iobs=155
	if ((method .eq. 'CSP'     .and. iobs .eq. 0))       iobs=5
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