character*30 f1,cc*1 real a_b(512*512) integer*2 b(512*512) integer*2 a(512*512*2) ! -32 change to 16 when cal. character*80 head(72) real*8 a8(8) common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut pi=4.*atan(1.0) cx=pi/12. cy=pi/180. phi=(40.+23/60.+36/3600.)*pi/180. ! XL latitude sphi=sin(phi) cphi=cos(phi) k=iargc() if(k.eq.0)then write(*,*) write(*,*)' usage: cc512 filename [!]' write(*,*)' for WJH prism 3 peaks image' write(*,*)' position stand for red peak' write(*,*)' jiang 2006,1,6' write(*,*) write(*,*) stop endif c check file call getarg(1,f1) if(inqh(f1,n1,n2,head,ibit).eq.0)stop 'DATA not found' if(n1.ne.512.or.n2.ne.512)stop 'not 512*512 CCD' call getpar(head) ! get s0,iy,im,id,ut idone=0 if(indexpos(head,'A81 ').ne.73)then idone=1 ipos=indexpos(head,'A81 ')-1 do 11 i=1,8 11 read(head(ipos+i)(11:),*)a8(i) endif if(k.ge.2)then call getarg(2,cc) if(cc.eq.'!')idone=0 endif if(idone.eq.1)stop 'coord already done!' if(kgetade(head).ne.0)stop 'bad coordination of fits_head' c read data 10 read(head(3)(65:),'(i5)')ii ii=ii/2+1 call readata(a,b,f1,ibit,ii) do 12 i=1,n1*n2 12 a_b(i)=a(i) iturn=0 if((iy-1900)*10+im.ge.945.or.iy.lt.1990)iturn=1 c turn image m1=512 m2=512 ipos=indexpos(head,'INSTRUME') if(ipos.ne.73)then do i=12,20 if(head(ipos)(i:i).eq.'3')iturn=0 if(head(ipos)(i:i).eq.'4')iturn=4 enddo endif if(iturn.eq.0)call turnx2(b,m1,m2,-180) if(iturn.eq.1)call turnx2(b,m1,m2,-360) if(iturn.eq.4)call turnx2(b,m1,m2,180) c main_job call main_job(head,a,a(512*512+1),a8,b,f1,iturn, * match,rms,seeing,white,a_b) c write something to fits head call put_a8(head,a8,match,rms,seeing,white) ! 45-52 write down a8 call put_gr(head) ! 53-57 write gal_co. redden call put_xmass(head) ! 15-17, 58 write center,h_angle,airmass call put_aa(head) ! 59,60 write Azimuth, Altitude call put_moon(head,a8) ! 61-65 phase, position, pos_angle close (1) c write(*,*) ' re_write fits head !' write(*,'(a,"GSC:",i4," RMS:",f6.2)')f1,match,rms i=lnblnk(f1)+1 f1(i:i)=char(0) call tvcoordc(f1,head) end function kgetade(head) character*80 head(72) character*1 c character*20 line common /platec/alpha_c,delta_c,epoch kgetade=1 ipos=indexpos(head,'RA ') if(ipos.eq.73)return line=head(ipos)(12:31) line(20:20)=char(0) do 15 i=1,19 if(line(i:i).eq.':')line(i:i)=',' 15 if(ichar(line(i:i)).eq.39)line(i:i)=' ' read(line(1:),*,err=10)i,j,x alpha_c=i+j/60.+x/3600. c write(*,*)'alpha',alpha_c c i=19 c if(head(ipos)(28:28).eq.'.')i=20 c if(head(ipos)(29:29).eq.'.')i=21 c read(head(ipos)(i:),'(i2,1x,i2,1x,f4.1)',err=10)i,j,x c alpha_c=i+j/60.+x/3600. c write(*,*)'alpha',alpha_c ipos=indexpos(head,'DEC ') if(ipos.eq.73)return line=head(ipos)(12:31) line(20:20)=char(0) c=' ' do 20 i=1,19 if(line(i:i).eq.'-')c='-' if(line(i:i).eq.'-')line(i:i)=' ' if(line(i:i).eq.':')line(i:i)=',' 20 if(ichar(line(i:i)).eq.39)line(i:i)=' ' read(line(1:),*,err=10)i,j,x delta_c=i+j/60.+x/3600. if(c.eq.'-')delta_c=-delta_c c write(*,*)'delta',delta_c c i=19 c if(head(ipos)(29:29).eq.'.')i=20 c if(head(ipos)(i+1:i+1).eq.'-')then c read(head(ipos)(i+1:),'(a,i1,1x,i2,1x,f4.1)',err=10)c,i,j,x c else c read(head(ipos)(i:),'(a,i2,1x,i2,1x,f4.1)',err=10)c,i,j,x c endif c delta_c=i+j/60.+x/3600. c if(c.eq.'-')delta_c=-delta_c c write(*,*)'delta',delta_c ipos=indexpos(head,'EPOCH ') if(ipos.eq.73)write(*,*)' no keywords EPOCH =' if(ipos.eq.73)then epoch=1995.1 else i=24 if(head(ipos)(29:29).eq.'.')i=25 read(head(ipos)(i:),'(f6.1)',err=10)epoch if(epoch.lt.1950.or.epoch.gt.2050.)return endif kgetade=0 return 10 write(*,*)' wrong R.A. or DEC data' end subroutine getstar(a,b,iturn,starx,stary,starv,nstar,white,seeing) integer*2 a(512,512),b(512*512) real starx(1),stary(1),starv(1) real see400(400) do 30 i=1,nstar see400(i)=see(b,512,512,starx(i),stary(i),white) starx(i)=513.-starx(i) 30 stary(i)=513.-stary(i) do 120 i=1,nstar-1 do 120 j=i+1,nstar if(starv(i).ge.starv(j))goto 120 x=starx(i) starx(i)=starx(j) starx(j)=x x=stary(i) stary(i)=stary(j) stary(j)=x x=starv(i) starv(i)=starv(j) starv(j)=x 120 continue seeing=xmedian(see400,nstar) end function see(map,n1,n2,x,y,white) c (a81+a84)/2 *180*3600/pi = arcsec / pixel c schmidt is 1.709"/pixel c seeing= 2.355*sigma*1.709 integer*2 map(n1,n2) integer b(113) ! (15-1)*8+1 ix=x+0.5 iy=y+0.5 j=1 do 10 i=ix-7,ix+7 b(j)=map(i,iy)-white 10 j=j+8 do 20 j=1,112,8 20 b(j+4)=(b(j)+b(j+8))*0.5+0.5 do 30 j=1,112,4 30 b(j+2)=(b(j)+b(j+4))*0.5+0.5 do 40 j=1,112,2 40 b(j+1)=(b(j)+b(j+2))*0.5+0.5 call histat(b,mode,maxh,xmean,xpeak,fsigma1,gsigma,113) j=1 do 110 i=iy-7,iy+7 b(j)=map(ix,i)-white 110 j=j+8 do 120 j=1,112,8 120 b(j+4)=(b(j)+b(j+8))*0.5+0.5 do 130 j=1,112,4 130 b(j+2)=(b(j)+b(j+4))*0.5+0.5 do 140 j=1,112,2 140 b(j+1)=(b(j)+b(j+2))*0.5+0.5 call histat(b,mode,maxh,xmean,xpeak,fsigma2,gsigma,113) c write(*,'(a,f9.2,a)')'seeing:',(fsigma1+fsigma2)/16.*2.355*1.709,'"' see=(fsigma1+fsigma2)/16.*2.355*1.709 end subroutine getmxy(a8,x,y,z,phase,pl,pb,xx) real*8 a8(8),b8(12),jd common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut call xytoad(a8,b8) j=(iy+((im-1)*30.4+id)/365.)*10.+0.5 epoch_g=j*0.1 x=a8(7)/cx y=a8(8)/cy call astprs(x,y,2000.,alpha_g,delta_g,epoch_g) call ad_azl(alpha_g,delta_g,a1,d1) alpha_g=alpha_g*cx delta_g=delta_g*cy call azl_ad(a1,d1-0.1,x,y) call standc(alpha_g,delta_g,x*cx,y*cy,xi,xn) a3=b8(1)*xi+b8(3)*xn+b8(5) d3=b8(2)*xi+b8(4)*xn+b8(6) call azl_ad(a1,d1+0.1,x,y) call standc(alpha_g,delta_g,x*cx,y*cy,xi,xn) a2=b8(1)*xi+b8(3)*xn+b8(5) d2=b8(2)*xi+b8(4)*xn+b8(6) z=atan2(a2-a3,d2-d3)/cy call fjd(jd,iy,im,id) jd=jd+ut/24. call moon2(jd,a2,d2) ! get monn l,b call lb_ad(a2,d2,x,y) ! get moon alpha delta phase=amod(a2-sunlamda(jd)+360.,360.)/360.*29.53 ! get moon phase call ad_azl(x,y,a2,d2) ! get moon A, H pl=a2 pb=d2 x=sin(d1*cy)*sin(d2*cy)+cos(d1*cy)*cos(d2*cy)*cos(a1*cy-a2*cy) y=acos(x) xx=y/cy cosA=(sin(d2*cy)-sin(d1*cy)*x)/cos(d1*cy)/sin(y) if(cosa.gt.0.999999)cosa=0.999999 if(cosa.lt.-.999999)cosa=-.999999 y=acos(cosA)/cy x=a2-a1 if(abs(x).gt.180.)x=a1-a2 if(x.lt.0.)y=-y z=270.+z-y z=amod(z,360.) if(z.lt.0.)z=z+360. x=cos(z*cy) y=sin(z*cy) end subroutine azl_ad(a,h,alpha,delta) !dd->hd common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut aa=a*cy hh=h*cy sina=sin(aa) cosa=cos(aa) sinh=sin(hh) cosh=cos(hh) sinde=sphi*sinh-cphi*cosa*cosh cosde=sqrt(1.-sinde*sinde) sint=cosh*sina/cosde cost=(cosa*cosh+cphi*sinde)/(sphi*cosde) delta=atan(sinde/cosde)/cy t=atan(sint/cost) if(cost.lt.0.)t=t+pi alpha=s0-t/cx if(alpha.lt.0.)alpha=alpha+24. if(alpha.gt.24.)alpha=alpha-24. end subroutine ad_azl(alpha,delta,a,h) ! hd,dd common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut hour=(s0-alpha)*cx delt=delta*cy s_delta=sin(delt) c_delta=cos(delt) s_hour=sin(hour) c_hour=cos(hour) sinh=sphi*s_delta+cphi*c_delta*c_hour cosh=sqrt(1.-sinh*sinh) sinA=c_delta*s_hour/cosh cosA=(sphi*c_delta*c_hour-cphi*s_delta)/cosh h=atan(sinh/cosh)/cy a=atan(sinA/cosA) if(cosA.lt.0.)a=a+pi if(a.lt.0.)a=a+pi+pi if(a.gt.pi+pi)a=a-pi-pi a=a/cy end subroutine put_aa(head) ! 59,60 character*80 head(72) common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut ipos=indexpos(head,'HA ') read(head(ipos)(19:),'(i2,1x,i2,1x,f5.2)')i,j,x t=i+j/60.+x/3600. ipos=indexpos(head,'DEC ') read(head(ipos)(20:),'(i2,1x,i2,1x,f4.1)')i,j,x d=i+j/60.+x/3600. if(head(ipos)(19:19).eq.'-')d=-d call ad_azl(s0-t,d,a,h) ! hd,dd ipos=indexpos(head,'A81 ')+14 ! 59 row head(ipos)(1:32)='AZIMUTH = /' head(ipos+1)(1:32)='ALTITUDE= /' head(ipos)(34:76)='Azimuth is measured from south through west' head(ipos+1)(34:54)='both units in degrees' write(head(ipos)(21:30),'(f10.2)')a ! in degree write(head(ipos+1)(21:30),'(f10.2)')h end subroutine put_moon(head,a8) ! 61-65 character*80 head(72) real*8 a8(8) common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut ipos=indexpos(head,'A81 ')+16 ! 61 row head(ipos)(1:32)='MPHASE = /' head(ipos+1)(1:32)='MAZIMUTH= /' head(ipos+2)(1:32)='MALITIUD= /' head(ipos+3)(1:32)='MANGLE = /' head(ipos+4)(1:32)='MDIRECT = /' head(ipos)(34:71)='Define the length of Moon_Month 29.53d' head(ipos+1)(34:78)='Moon Azimuth is measured from south thr. west' head(ipos+2)(34:54)='both units in degrees' head(ipos+3)(34:75)='Position angle of Moon to CCD_Field center' head(ipos+4)(34:77)='Moon direction in the X_Y plane of CCD frame' call getmxy(a8,x,y,z,phase,pl,pb,xx) write(head(ipos)(21:30),'(f10.1)')phase write(head(ipos+1)(21:30),'(f10.2)')pl write(head(ipos+2)(21:30),'(f10.2)')pb write(head(ipos+3)(21:30),'(f10.2)')xx write(head(ipos+4)(21:30),'(f10.2)')z end function sunlamda(jd) ! in degree c sun lamda in degree, beta==0 real*8 jd,n,l,g,x n=jd-2451545.d0 l=280.460d0+0.9856474d0*n g=357.528d0+0.9856993d0*n x=dmod(g,360d0)*datan(1d0)/45d0 x=l+1.915d0*dsin(x)+0.020d0*dsin(x+x) x=dmod(x,360d0) if(x.lt.0d0)x=x+360d0 sunlamda=x end subroutine moon2(jd,pl,pb) !in degree real*8 jd,t,g,r,x1,x2,x3,x4,x5,x6 g=datan(1d0)/45d0 t=(jd-2451545.d0)/36525.d0 r=218.32d0+dmod(481267.833d0*t,360d0) x1=(134.9d0+dmod(477198.85d0*t,360d0))*g x2=(259.2d0-dmod(413335.38d0*t,360d0))*g x3=(235.7d0+dmod(890534.23d0*t,360d0))*g x4=(269.9d0+dmod(954397.70d0*t,360d0))*g x5=(357.5d0+dmod( 35999.05d0*t,360d0))*g x6=(186.6d0+dmod(966404.05d0*t,360d0))*g r=r+6.29d0*dsin(x1)-1.27d0*dsin(x2)+0.66d0*dsin(x3) * +0.21d0*dsin(x4)-0.19d0*dsin(x5)-0.11d0*dsin(x6) r=dmod(r,360d0) if(r.lt.0d0)r=r+360d0 pl=r x1=( 93.3d0+dmod(483202.03d0*t,360d0))*g x2=(228.2d0+dmod(960400.87d0*t,360d0))*g x3=(318.3d0+dmod( 6003.18d0*t,360d0))*g x4=(217.6d0-dmod(407332.20d0*t,360d0))*g r=5.13d0*dsin(x1)+0.28d0*dsin(x2)-0.28d0*sin(x3)-0.17d0*dsin(x4) r=dmod(r,360d0) if(r.gt.90d0)r=r-360d0 pb=r end subroutine lb_ad(pl,pb,alpha,delta) !dd-->hd cy=atan(1.)/45. xl=cos(pb*cy)*cos(pl*cy) xm=0.9175*cos(pb*cy)*sin(pl*cy)-0.3978*sin(pb*cy) xn=0.3978*cos(pb*cy)*sin(pl*cy)+0.9175*sin(pb*cy) alpha=atan2(xm,xl)/cy/15. if(alpha.lt.0.)alpha=alpha+24. delta=asin(xn)/cy end subroutine ad_lb(al,de,pl,pb) !hd-->dd cy=atan(1.)/45. cx=cy*15. xl=cos(de*cy)*cos(al*cx) xm=0.9175*cos(de*cy)*sin(al*cx)+0.3978*sin(de*cy) xn=-.3978*cos(de*cy)*sin(al*cx)+0.9175*sin(de*cy) pl=atan2(xm,xl)/cy if(pl.lt.0.)pl=pl+360. pb=asin(xn)/cy end subroutine put_a8(head,a8,match,rms,seeing,white) !45-52 character*80 head(72),ac*12,dc*12 real*8 a8(8) common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut ipos=indexpos(head,'VOLT3 ') if(ipos.ne.73)then write(head(ipos)(24:30),'(f7.2)')seeing head(ipos)(34:49)='SEEING (arcsec)' endif ipos=indexpos(head,'VOLT2 ') if(ipos.ne.73)then write(head(ipos)(22:),"(f9.2,' / SKY ADU/PIXEL')")white endif ipos=indexpos(head,'A81 ') if(ipos.eq.73)ipos=indexpos(head,' ') if(ipos.gt.51)ipos=indexpos(head,'COMMEMT ') if(ipos.gt.51)ipos=51 head(ipos)(1:46)='A81 = / 6 plate coef.' head(ipos+1)(1:46)='A82 = / matched star:' head(ipos+2)(1:46)='A83 = / RMS (arcsec):' head(ipos+3)(1:32)='A84 = /' head(ipos+4)(1:32)='A85 = /' head(ipos+5)(1:32)='A86 = /' head(ipos+6)(1:46)='A87 = / p_c (2000.0) ' head(ipos+7)(1:46)='A88 = / in rad. ' do 11 i=1,6 11 write(head(ipos-1+i)(11:30),'(E20.12)')a8(i) write(head(ipos+6)(11:30),'(f20.7)')a8(7) write(head(ipos+7)(11:30),'(f20.7)')a8(8) write(head(ipos+1)(47:52),'(i6)')match write(head(ipos+2)(47:52),'(f6.2)')rms x=a8(7)/cx ! put lamda beta y=a8(8)/cy call ad_lb(x,y,alpha,delta) call toms(alpha,ac,2) call toms(delta,dc,0) head(ipos+3)(57:63)='l_sun =' head(ipos+4)(57:63)='b_sun =' head(ipos+3)(65:73)=ac(1:9) head(ipos+4)(65:73)=dc(1:9) end subroutine put_gr(head) !51:53-57 character*80 head(72) character*12 ac,dc common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut ipos=indexpos(head,'A81 ')+6 read(head(ipos)(21:),*)xa ! center plate in rad. read(head(ipos+1)(21:),*)xd xa=xa/cx xd=xd/cy ! to h,d call toms(xa,ac,1) call toms(xd,dc,0) head(ipos+2)(1:46)="RA2000 = ' ' / plate center " head(ipos+3)(1:32)="DEC2000 = ' ' /" head(ipos+2)(18:29)=ac head(ipos+3)(19:29)=dc(1:11) call astprs(xa,xd,2000.,x,y,1950.) call galactic(x*cx,y*cy,xa,xd) head(ipos+4)(1:32)='GALLONG = /' head(ipos+4)(34:62)='galactic coordinates (1950.0)' head(ipos+5)(1:32)='GALLATI = /' head(ipos+5)(34:54)='both units in degrees' write(head(ipos+4)(21:30),'(f10.2)')xa ! in degree write(head(ipos+5)(21:30),'(f10.2)')xd head(ipos+6)(1:32)='EXTINCT = /' head(ipos+6)(34:63)='BH Bmag, EXTINCTION = 4*E(B-V)' xa=redden(xa,xd) write(head(ipos+6)(21:30),'(f10.3)')xa end subroutine put_xmass(head) !15,16,17,58 character*80 head(72) character*12 ac,dc common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut ipos=indexpos(head,'A81 ')+6 read(head(ipos)(21:),*)xa ! center plate in rad. read(head(ipos+1)(21:),*)xd xa=xa/cx xd=xd/cy ! to h,d j=(iy+((im-1)*30.4+id)/365.)*10.+0.5 epoch=j*0.1 call astprs(xa,xd,2000.,x,y,epoch) call toms(x,ac,1) call toms(y,dc,0) ipos=indexpos(head,'RA ') head(ipos)(11:30)="' '" head(ipos)(18:29)=ac write(head(ipos)(50:),"('(',f6.1,')')")epoch ipos=indexpos(head,'EPOCH ') if(ipos.ne.73)head(ipos)(11:30)="' '" if(ipos.ne.73)write(head(ipos)(24:29),'(f6.1)')epoch ipos=indexpos(head,'DEC ') head(ipos)(11:30)="' '" head(ipos)(19:29)=dc(1:11) write(head(ipos)(50:),"('(',f6.1,')')")epoch h=amod(s0-x,24.) ! get HA, here I wast if(h.lt.0.)h=h+24. ! more time, x->xa call toms(h,ac,1) ipos=indexpos(head,'HA ') head(ipos)(11:30)="' '" head(ipos)(18:29)=ac xmm=xmass(h*cx,xd*cy) ! middle time ipos=indexpos(head,'EXPOSURE') read(head(ipos)(20:),*)k ! exp time xmb=xmass((h-k/7200.)*cx,xd*cy) ! begin time xme=xmass((h+k/7200.)*cx,xd*cy) ! end time x=(xmb+4.*xmm+xme)/6. ipos=indexpos(head,'A81 ')+13 head(ipos)(1:32)='AIRMASS = /' write(head(ipos)(21:30),'(f10.3)')x end subroutine getpar(head) c get s0 iy im id ut character*80 head(72),c9*9 real*8 d1,d common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut ipos=indexpos(head,'DATE-OBS') read(head(ipos)(22:),'(i2,1x,i2,1x,i2)')id,im,iy ! read obs date iy=iy+1900 if(iy.lt.1950)iy=iy+100 ipos=indexpos(head,'TIME ') read(head(ipos)(20:),'(i2,1x,i2,1x,f4.1)')i,j,x ipos=indexpos(head,'EXPOSURE') read(head(ipos)(20:),*)k ! exp time ut=i+j/60.+x/3600.+k/7200. ! center UT time, in h call fjd(d1,1992,12,31) ! 1992,12,31 d2=(6+38/60.+40.1954/3600.)/24. ! 6 38 40.1954 x1=7+50/60.+18.344/3600. ! 7:50:18.344 lamda of XL call fjd(d,iy,im,id) d=(d-d1)*1.0027379093d0+d2 d=dmod(d,1d0) if(d.lt.0d0)d=d+1d0 s0=d*24d0 s0=s0+x1+ut*1.002738 end subroutine fjd(d,iy,im,id) real*8 d i=(im+9)*0.09 i=(i+iy-1900)*1.75+0.01 j=30.56*im j=(iy-1950)*367+id+j-i d=2433338.5d0+j end function xmass(angle,delta) common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut x=sin(delta)*sphi+cos(delta)*cos(angle)*cphi c z=90.-atan2(x,sqrt(1.-x*x))*180./pi ! Zenth x=1./x-1. xmass=1.+x-x*(0.0018167+x*(0.002875+x*0.0008083)) end subroutine put_a2(head) character*80 head(72) character*35 dis common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut ipos=indexpos(head,'A81 ') if(ipos.eq.73)ipos=indexpos(head,' ') if(ipos.gt.51)ipos=indexpos(head,'COMMEMT ') if(ipos.gt.51)ipos=51 head(ipos)(1:32)='A81 = /' head(ipos)(34:53)='6 plate coefficients' head(ipos+1)(1:32)='A82 = /' head(ipos+2)(1:32)='A83 = /' head(ipos+3)(1:32)='A84 = /' head(ipos+4)(1:32)='A85 = /' head(ipos+5)(1:32)='A86 = /' head(ipos+6)(1:46)='A87 = / p_c (2000.0) ' head(ipos+7)(1:46)='A88 = / in radians ' 10 write(*,*)'input center R.A, Decl, [1950.0]: ex. 3:46:02.1 23:0' read(*,'(a)')dis if(itohd(dis,x,y,epoch).eq.0)goto 10 if(epoch.eq.0.)epoch=1950. call astprs(x,y,epoch,x,y,2000.) write(head(ipos+6)(11:30),'(f20.7)')x*cx write(head(ipos+7)(11:30),'(f20.7)')y*cy end function redden(al,b) c this subroutine written by burstein. INTEGER*4 IRED(1200),IHI(201) XRAD = 57.2958 ALIN = AL/0.3 + 0.5 ILIN = ALIN + 0.51 BB = ABS(B) redden= -0.99 if(bb.lt.10.) return OPEN(UNIT=11,FILE='/EOD/bur/redsouth.dat', * STATUS='OLD',ACCESS='DIRECT',recl=4800) OPEN(UNIT=12,FILE='/EOD/bur/rednorth.dat', * STATUS='OLD',ACCESS='DIRECT',recl=4800) OPEN(UNIT=13,FILE='/EOD/bur/hinorth.dat', * STATUS='OLD',ACCESS='DIRECT',recl=804) OPEN(UNIT=14,FILE='/EOD/bur/hisouth.dat', * STATUS='OLD',ACCESS='DIRECT',recl=804) IF (B) 11,11,18 11 IF (B + 62.) 13,12,12 12 AREC = -(B+10.)/0.6 + 1. KZ = AREC + 0.51 READ (11,rec=KZ) IRED call swap4(ired,4800) KQ = KZ -1 IEBV = IRED(ILIN) GO TO 19 13 ALC = AL/XRAD AREC = 101. + SIN(ALC)*(90.+B)/0.3 ALIN = 101. + COS(ALC)*(90.+B)/0.3 KW = AREC + 0.51 ILIN = ALIN + 0.51 READ (14,rec=KW) IHI call swap4(ihi,804) AHI = IHI(ILIN) BMV = -0.0372 + 0.357*AHI/10000. IEBV = BMV*1000. + 0.5 KQ = KW - 1 GO TO 19 18 IF (B - 62.) 17,17,16 17 AREC = (B-10.)/0.6 + 1. KY = AREC + 0.51 READ (12,rec=KY) IRED call swap4(ired,4800) KQ = KY - 1 IEBV = IRED(ILIN) GO TO 19 16 ALC = AL/XRAD AREC = 101. + SIN(ALC)*(90.-B)/0.3 ALIN = 101. + COS(ALC)*(90.-B)/0.3 KX = AREC + 0.51 ILIN = ALIN + 0.51 READ (13,rec=KX) IHI call swap4(ihi,804) AHI = IHI(ILIN) BMV = -0.0372 + 0.357*AHI/10000. IEBV = BMV*1000. + 0.5 KQ = KX - 1 19 ABV = IEBV redden= ABV*4./1000. + 0.005 close(11) close(12) close(13) close(14) return END subroutine galactic(alpha,delta,xl,xb) common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut a0=282.25*cy theta=62.6*cy Xl0=33.0*cy x1=cos(delta)*cos(alpha-a0) x2=cos(delta)*sin(alpha-a0)*cos(theta)+sin(delta)*sin(theta) xl=atan2(x2,x1)+xl0 if(xl.lt.0.)xl=xl+2.0*pi xb=asin(sin(delta)*cos(theta)- * cos(delta)*sin(alpha-a0)*sin(theta)) xl=xl/cy xb=xb/cy end subroutine readata(a,b,f1,ibit,ii) integer c(720),a(1),b(1) character*30 f1,ch*2880 equivalence (c,ch) nn=720 open(1,file=f1,status='old',access='direct',recl=nn*4) read(1,rec=1)ch n1=512 n2=512 n=1 55 k=index(ch,'END ') if(k.eq.0)then n=n+1 read(1,rec=n)ch goto 55 endif c write(*,"(' reading data ')") i=0 if(ibit.lt.0)ibit=-ibit k1=n+1 k2=n+(n1*n2*ibit/32-1)/nn do 60 k=k1,k2 c if(mod(k,120).eq.0)call dispdot() read(1,rec=k)ch do 60 j=1,nn i=i+1 60 a(i)=c(j) close (1) if(ibit.eq.32)then call swap4(a,n1*n2*4) call shrink4(a,n1,n2,a,ii) else call swap2(a,n1*n2*2) endif call shrink2(a,n1,n2,b) ! 2048*2048-->512*512 c write(*,*) end subroutine shrink2(a,n1,n2,b) integer*2 a(n1,n2) integer*2 b(512*512) k=0 m=1 do 10 j=1,n2,m do 10 i=1,n1,m k=k+1 10 b(k)=a(i,j) end subroutine shrink4(a,n1,n2,b,ii) real a(n1*n2) integer*2 b(n1*n2) do 10 i=1,n1*n2 x=a(i)/ii if(x.gt.32767.)x=32767. if(x.lt.0.)x=0. 10 b(i)=x end function inqh(f1,n1,n2,head,ibit) character*(*) f1 character*80 head(72) logical inqh k=index(f1,'.') if(k.eq.0)f1=f1(1:lnblnk(f1))//'.fit' inquire(file=f1,exist=inqh) if(.not.inqh)return open(1,file=f1,status='old', * access='direct',recl=5760) read(1,rec=1)head read(head(2)(22:),*)ibit read(head(4)(22:),*)n1 !NAXIS1 read(head(5)(22:),*)n2 !NAXIS2 c ipos=indexpos(head,'OBJECT ') c write(*,'(1x,a)')head(ipos)(1:70) ! object ipos=indexpos(head,'DATE-OBS') head(ipos)(63:72)='(dd/mm/yy)' c do 10 i=ipos,ipos+5 ! da.ti.ex.ra.de.ha c10 write(*,'(1x,a)')head(i)(1:78) c ipos=indexpos(head,'EPOCH ') c if(ipos.ne.73)write(*,'(1x,a)')head(ipos)(1:78) ipos=indexpos(head,'TELESCOP') if(ipos.ne.73)head(ipos)(12:28)='BAO 60-90 SCHIMDT' ipos=indexpos(head,'BZERO ') if(ipos.ne.73)head(ipos)(1:1)='P' ipos=indexpos(head,'BSCALE ') if(ipos.ne.73)head(ipos)(1:1)='P' ipos=indexpos(head,'TIME ') if(ipos.ne.73)head(ipos)(42:57)='OF START OF OBS.' ipos=indexpos(head,'HA ') if(ipos.ne.73)head(ipos)(45:61)='OF MIDDLE OF OBS.' close (1) end subroutine main_job(head,map,map1,a8,ma,f1,iturn, * match,rms,seeing,white,a_b) real a_b(512*512) character*80 head(72) integer*2 map(512,512),ma(512,512) ! map buffer integer*2 map1(512,512) ! 2nd buffer real*8 a8(8) character*12 ac,dc,c3*3 ! alpha delta display line character*50 ns,dis*60 ! for display real x9(99),y9(99) ! al_del of starard posi star real sx(4),sy(4) ! scale on right_bottom of pic real*8 adcoef(2,6),xycoef(2,6),dummy8 ! 2 set 6 plte coef real starx(400),stary(400),starv(400) ! take ccd stars byte starm(400) parameter (ngsc=3999) ! No. of gsc stars. real aagsc(ngsc),ddgsc(ngsc),mmgsc(ngsc) character f1*30 common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut common /platec/alpha_c,delta_c,epoch call whitexblack(ma,512,512,white,x) black=white+20.*x n1=512 n2=512 ir=4 ******** auto match call findstar(head,a_b,starx,stary,starv,nstar) call sortxyv(starx,stary,starv,nstar) if(nstar.gt.200)nstar=200 call sortxyv(stary,starv,starx,nstar) call couple(starx,stary,starv,starm,nstar) call getstar(map1,map,iturn,starx,stary,starv,nstar,white,seeing) ipp=0 c******** take out gsc stars epoch1=epoch alpha1=alpha_c delta1=delta_c 1011 call fgsc(alpha_c,delta_c,epoch,20.,aagsc,ddgsc,mmgsc,ng) call astprs(alpha_c,delta_c,epoch,xa,xd,2000.) epoch=2000. xycoef(1,1)=-0.83e-5 xycoef(2,1)= 0.10e-6 xycoef(1,2)=-0.10e-6 xycoef(2,2)=-0.83e-5 xycoef(1,3)= 0.21e-2 xycoef(2,3)= 0.21e-2 alpha_c=xa*cx delta_c=xd*cy if(ng.gt.30)ng=30 rms=100. nst=nstar if(nst.gt.30)nst=30 match=kgsc(aagsc,ddgsc,mmgsc,ng, * starx,stary,starv,nst, * xycoef,alpha_c,delta_c,rms) if(match.le.5)then xycoef(1,1)=-0.83e-5 xycoef(2,1)= 0.10e-6 xycoef(1,2)=-0.10e-6 xycoef(2,2)=-0.83e-5 xycoef(1,3)= 0.21e-2 xycoef(2,3)= 0.21e-2 call shift_c(starx,stary,aagsc,ddgsc,ng,xycoef,alpha_c,delta_c) xycoef(1,1)=-0.83e-5 xycoef(2,1)= 0.10e-6 xycoef(1,2)=-0.10e-6 xycoef(2,2)=-0.83e-5 xycoef(1,3)= 0.21e-2 xycoef(2,3)= 0.21e-2 match=kgsc(aagsc,ddgsc,mmgsc,ng, * starx,stary,starv,nst, * xycoef,alpha_c,delta_c,rms) endif call fgsc(alpha_c/cx,delta_c/cy,epoch,20.,aagsc,ddgsc,mmgsc,ng) rms=5. match=kgsc(aagsc,ddgsc,mmgsc,ng, * starx,stary,starv,nstar, * xycoef,alpha_c,delta_c,rms) if(match.eq.0)stop 'match fail (star too faint) !' ********* save a8 a8(1)=xycoef(1,1) a8(2)=xycoef(2,1) a8(3)=xycoef(1,2) a8(4)=xycoef(2,2) a8(5)=xycoef(1,3) a8(6)=xycoef(2,3) ******95,12,19 add if(alpha_c.lt.0.)alpha_c=alpha_c+pi+pi a8(7)=alpha_c a8(8)=delta_c ! 2000 epoch ********* save a8 end end subroutine couple(x,y,v,m,n) real x(1),y(1),v(1) byte m(1) do 10 i=1,n 10 m(i)=0 kk=0 15 k=0 20 k=k+1 if(k.ge.n)goto 100 if(m(k).ne.0)goto 20 i1=k 30 k=k+1 if(k.gt.n)goto 100 if(m(k).ne.0)goto 30 i2=k z=y(i1)-y(i2) u=v(i1)/v(i2) if(x(i1)-x(i2).lt.1. .and. z.gt.16.5 .and. z.lt.17.4 c .and. u.gt.0.5 .and. u.lt.3.)then m(i1)=1 m(i2)=2 else k=i1 endif goto 20 100 kk=kk+1 if(kk.eq.1)goto 15 c do 110 i=1,nstar c110 write(*,*)x(i),y(i),v(i),m(i) k=0 do 200 i=1,n if(m(i).ne.1)goto 200 k=k+1 x(k)=x(i) y(k)=y(i) v(k)=v(i) 200 continue n=k end subroutine sortxyv(x,y,v,n) real x(1),y(1),v(1) do 10 i=1,n-1 do 10 j=i+1,n if(v(i).lt.v(j))then z=v(i) v(i)=v(j) v(j)=z z=x(i) x(i)=x(j) x(j)=z z=y(i) y(i)=y(j) y(j)=z endif 10 continue end subroutine xy_rad(x,y,a8,alpha,delta) c 2048*2048 ccd's X Y (in pixel) -------------> alpha delta (in h,d) real x,y,alpha,delta real*8 a8(8) ! from fits_head "A81"--"A88" pi=4.*atan(1.) xi=a8(1)*x+a8(3)*y+a8(5) xn=a8(2)*x+a8(4)*y+a8(6) rc=a8(7) ! center of ccd (in rad.) dc=a8(8) cd=cos(dc) td=tan(dc) alpha=atan(xi/cd/(1.-xn*td)) delta=atan((xn+td)*cos(alpha)/(1.-xn*td)) alpha=alpha+rc if(abs(dc-delta).gt.2.)then delta=-delta alpha=alpha+pi endif alpha=alpha*12./pi delta=delta*180./pi end function indexpos(head,f1) character*80 head(72),f1*8 do 10 indexpos=1,72 10 if(head(indexpos)(1:8).eq.f1)return end function kgsc(aa,dd,gg,ng,pxx,pyy,pzz,nc,xycoef,alpha_c,delta_c,rms) parameter (ngsc=400,nstar=300) real aa(1),dd(1),gg(1),pxx(1),pyy(1),pzz(1) real xx(nstar),yy(nstar),zz(nstar) c real*8 cf8(8) real wa(ngsc),wd(ngsc),gr(nstar),gd(nstar) character*12 ac,dc c character*12 ac1,dc1,ns*70,ch*1 byte xyw(nstar) real*8 xixn(2,nstar),xy(2,nstar),xycoef(2,6),adcoef(2,6) do 111 i=1,nc xx(i)=pxx(i) yy(i)=pyy(i) 111 zz(i)=pzz(i) n12=12 if(nc.le.30)n12=6 pi=4.*atan(1.0) cx=pi/12. cy=pi/180. epoch=2000. n1=512 n2=512 call xytoad(xycoef,adcoef) do 100 i=1,ng xa=aa(i)*cx xd=dd(i)*cy call standc(alpha_c,delta_c,xa,xd,xi,xn) wa(i)=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3) 100 wd(i)=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3) c100 write(*,*)aa(i),dd(i),wa(i),wd(i) c write(*,*)xycoef,alpha_c,delta_c c call pgbegin(0,'/xw',1,1) c call pgenv(1.,512.,512.,1.,1,0) c do 101 i=1,ng c k=36.-gg(i) c101 call pgpoint(1,wa(i),wd(i),k) c call pgsci(3) c do 102 i=1,nc c n9=alog(zz(i)/4000.)+20 c if(n9.lt.20)n9=20 c if(n9.gt.27)n9=27 c102 call pgpoint(1,xx(i),yy(i),n9) c c read(*,*) do 105 i=1,nc 105 xyw(i)=0 c m=0 do 110 i=1,ng x=rms do 120 j=1,nc z1=wa(i)-xx(j) z2=wd(i)-yy(j) z=abs(z1)+abs(z2) if(z.lt.x)then k=j x=z endif 120 continue if(x.ne.rms)then xyw(k)=1 gr(k)=aa(i) gd(k)=dd(i) c else c m=m+1 c call toms(aa(i),ac1,1) c call toms(dd(i),dc1,0) c write(*,'(i5,a,2f7.1,2(2x,a))')m, c * ' a gsc star lost',wa(i),wd(i),ac1,dc1 endif 110 continue j=0 do 152 i=1,nc if(xyw(i).ne.0)then j=j+1 xx(j)=xx(i) yy(j)=yy(i) xy(1,j)=xx(i) xy(2,j)=yy(i) gr(j)=gr(i) gd(j)=gd(i) xyw(j)=1 xa=gr(j)*cx xd=gd(j)*cy call standc(alpha_c,delta_c,xa,xd,xi,xn) xixn(1,j)=xi xixn(2,j)=xn endif 152 continue n=j k=0 kgsc=0 if(n.lt.n12)return ! Fail, too less star match 151 k=k+1 c if(n-k.lt.n12)return ! Fail, too less star left call plate(xixn,xy,xyw,n,xycoef,3) call toms(alpha_c/cx,ac,1) call toms(delta_c/cy,dc,0) c write(ns(1:),"('Center: ',2a,' (',f6.1,')')")ac,dc,epoch c call pgsch(.7) c call pgtext(0.,-45.,ns) c write(*,2)epoch,ac,dc c2 format(' **Frame Center: (',f6.1,')',/5x,2a) c write(*,*)(xycoef(1,i),i=1,3) c write(*,*)(xycoef(2,i),i=1,3) sig=0. m=0. z=0. c call pgsci(6) do 150 i=1,n xi=xycoef(1,1)*xy(1,i)+xycoef(1,2)*xy(2,i)+xycoef(1,3) xn=xycoef(2,1)*xy(1,i)+xycoef(2,2)*xy(2,i)+xycoef(2,3) call astand(alpha_c,delta_c,xi,xn,xa,xd) xa=xa/cx xd=xd/cy sigma=((gr(i)-xa)*15.*cos(delta_c))**2+(gd(i)-xd)**2 sigma=sqrt(sigma)*3600. if(xyw(i).eq.1)then if(sigma.gt.z)then z=sigma j=i endif sig=sig+sigma m=m+1 endif c call toms(gr(i),ac1,1) c call toms(gd(i),dc1,0) c call toms(xa,ac,1) c call toms(xd,dc,0) c if(i.eq.1)write(*,449)epoch,epoch c449 format(9x,'GSC (',f6.1,')',12x,'Result (',f6.1,')' c * /2x,24('-'),3x,24('-')) c ch=' ' c if(xyw(i).eq.1)ch='*' c write(ns(1:),'(i2)')i c call pgsci(2) c call pgsch(.7) c call pgtext(xx(i)-15.,yy(i)-15.,ns(1:2)) c call pgsch(1.) c if(xyw(i).ne.1)call pgsci(3) c call pgpoint(1,xx(i),yy(i),24) c write(*,'(i3,2a,2x,2a,g12.2,,4x,a)')i,ac1,dc1,ac,dc,sigma,ch 150 continue x=sig/m c write(*,'(62x,g12.2)')x c 3 choice here, after rather good result, must adjest plate center, use 0 c60 write(*,*)'1,2....n delete or recover this star to re_cal 6 coef.' c write(*,*)' 0 calculate center of frame' c write(*,*)'-1 quit' c read(*,*,err=60)k if(x.gt.6.0.and.k.gt.9)return ! Fail, too big error if(z.gt.2.0)then xyw(j)=-xyw(j) goto 151 endif k=-k xa=(n1+1)*0.5 xd=(n2+1)*0.5 xi=xycoef(1,1)*xa+xycoef(1,2)*xd+xycoef(1,3) xn=xycoef(2,1)*xa+xycoef(2,2)*xd+xycoef(2,3) call astand(alpha_c,delta_c,xi,xn,xa,xd) alpha_c=xa delta_c=xd do 141 i=1,n xa=gr(i)*cx xd=gd(i)*cy call standc(alpha_c,delta_c,xa,xd,xi,xn) xixn(1,i)=xi 141 xixn(2,i)=xn if(k.lt.0)goto 151 c cf8(1)=xycoef(1,1) c cf8(2)=xycoef(2,1) c cf8(3)=xycoef(1,2) c cf8(4)=xycoef(2,2) c cf8(5)=xycoef(1,3) c cf8(6)=xycoef(2,3) c epoch is 2000 c cf8(7)=alpha_c c cf8(8)=delta_c c write(*,*)'GSC:',ng,' CCD:',nc,' Matched:',n c if(ng.gt.30)write(*,3)n-k,x c3 format(9x,'By using',i4,' stars to calculate 8 ceof, RMS:', c * f5.2,' (arcsec)') kgsc=n-k rms=x return end subroutine fgsc(a_c,d_c,e_p,s_z,aa,dd,xmag,n) parameter (ngsc=3999) ! No. of local area stars. character*6 dir,ns*30 ! gsc directory name character*8 file1,file2,file3,file4 ! gsc file name real aa(1),dd(1),xmag(1) ! alpha delta for gsc integer*2 ee(ngsc),ingsc(ngsc),e2 byte w910(2) character col(ngsc),color,color1 logical logi equivalence (w910,e2) nn=ngsc pi=4.*atan(1.) cx=pi/12. cy=pi/180. ns='/EOD/GSC/regions.bin' inquire(file=ns,exist=logi) if(.not.logi)stop 'GSC cataloge not found !' alph_a=a_c delt_a=d_c if(delta.gt. 89.5)delta= 89.5 ! near 90 is too large cal.err if(delta.lt.-89.5)delta=-89.5 call astprs(alph_a,delt_a,e_p,alpha,delta,2000.) c get 4 corner position of pic. unit is hour & degree hw_de=s_z/120. up_d=delta+hw_de dn_d=delta-hw_de hs=hw_de*cy hs2=511./hs*0.5 x=min(cos(up_d*cy),cos(dn_d*cy)) hw_al=hw_de/x/15. if(hw_al.gt.12..or.x.lt.0.)hw_al=12. up_a=alpha+hw_al dn_a=alpha-hw_al c for R.A. to joint 24 hour and 0 hour c iraf=0 normal, iraf=1 center R.A. is large, iraf=-1 center R.A. near 0 hour iaf=0 if(up_a.ge.24.) iaf=1 if(up_a.ge.24.) up_a=up_a-24. if(dn_a.lt.0.) iaf=-1 if(dn_a.lt.0.) dn_a=dn_a+24. **************************FIND STAR*********************** n=0 c if delta>83.5 use all N8230 directory files if(abs(delta).gt.83.5)then if(delta.gt.0.)then dir='N8230/' do 20 i=4615,4662 ! file name is 4615.BAO~~4662.BAO write(file1(1:),"(i4,'.BAO')")i 20 call find(dir,file1,aa,dd,ee,col,nn,n,up_d,dn_d,up_a,dn_a,iaf) else dir='S8230/' do 21 i=9490,9537 ! as north write(file1(1:),"(i4,'.BAO')")i 21 call find(dir,file1,aa,dd,ee,col,nn,n,up_d,dn_d,up_a,dn_a,iaf) endif else c if delta <83.5 use 4 files or lease. c according to 4 corner to index which file should be read. if(jindex(up_a,up_d,dir,file1).ne.0)goto 30 call find(dir,file1,aa,dd,ee,col,nn,n,up_d,dn_d,up_a,dn_a,iaf) c write(*,'(a)')file1 30 if(jindex(up_a,dn_d,dir,file2).ne.0)goto 31 if(file2.eq.file1)goto 31 call find(dir,file2,aa,dd,ee,col,nn,n,up_d,dn_d,up_a,dn_a,iaf) c write(*,'(a)')file2 31 if(jindex(dn_a,up_d,dir,file3).ne.0)goto 32 if(file3.eq.file1.or.file3.eq.file2)goto 32 call find(dir,file3,aa,dd,ee,col,nn,n,up_d,dn_d,up_a,dn_a,iaf) c write(*,'(a)')file3 32 if(jindex(dn_a,dn_d,dir,file4).ne.0)goto 33 if(file4.eq.file1.or.file4.eq.file2.or.file4.eq.file3)goto 33 call find(dir,file4,aa,dd,ee,col,nn,n,up_d,dn_d,up_a,dn_a,iaf) c write(*,'(a)')file4 33 continue endif c find star end. but in some case a star appear several time in c different files, find it, keep only one, ignore repeat one. c GSC says, a star at more appears 7 times in different pltes. c according to positions' error to detect the repeat one. nn=n c 1. determine how many stars cosz=15.*cos(delta*cy) n=1 do 34 i=1,nn-1 ingsc(i)=n x=((aa(i)-aa(i+1))*cosz)**2 y=(dd(i)-dd(i+1))**2 34 if(x+y.gt.5e-6)n=n+1 c accordind to test, 1e-6 or 1e-5 all OK, so choice middle value ingsc(nn)=n c 2. acording color, change positions of alpha and delta color1='V' i=0 3410 i=i+1 if(i.ge.nn)goto 3430 if(ingsc(i).eq.ingsc(i+1))then if(col(i).eq.color1)then j=i+1 x=aa(j) y=dd(j) k=ee(j) color=col(j) aa(j)=aa(i) dd(j)=dd(i) ee(j)=ee(i) col(j)=col(i) aa(i)=x dd(i)=y ee(i)=k col(i)=color endif 3420 if(i.ge.nn-1)goto 3430 if(ingsc(i+1).eq.ingsc(i+2))then i=i+1 goto 3420 endif endif goto 3410 3430 continue c 3. put each stars in front, mult stars in last do 341 i=1,n if(ingsc(i).ne.i)then do 342 j=i+1,nn if(ingsc(j).eq.i)goto 343 342 continue 343 x=aa(j) y=dd(j) k=ee(j) c color=col(j) aa(j)=aa(i) dd(j)=dd(i) ee(j)=ee(i) c col(j)=col(i) ingsc(j)=ingsc(i) aa(i)=x dd(i)=y ee(i)=k c col(i)=color ingsc(i)=i endif 341 continue c 4. delete out of frame star k=0 alph_a=alpha*cx ! center of pic. delt_a=delta*cy do 344 i=1,n xa=aa(i)*cx xd=dd(i)*cy call standc(alph_a,delt_a,xa,xd,xi,xn) ix=hs2*(hs-xi)+1.5 iy=hs2*(xn-hs)+512.5 if(inpic(ix,iy).eq.0)goto 344 c release last 2 byte message e2=ee(i) mag=w910(1) if(mag.lt.0)x=mag*0.1-5. if(mag.gt.0)x=mag*0.1+5. if(mag.lt.0)x=-x c if(mag.lt.0)col(i)=col(i)+32 k=k+1 aa(k)=aa(i) dd(k)=dd(i) xmag(k)=x c col(k)=col(i) c ingsc(k)=ingsc(i) 344 continue c write(*,*)nn,k c 5. 2nd delete same position star n=1 do 400 i=2,k do 410 j=1,n x=((aa(i)-aa(j))*cosz)**2 y=(dd(i)-dd(j))**2 if(x+y.lt.5e-6)then ! keep bright star if(xmag(j).le.xmag(i))goto 400 xmag(j)=xmag(i) aa(j)=aa(i) dd(j)=dd(i) goto 400 endif 410 continue n=n+1 aa(n)=aa(i) dd(n)=dd(i) xmag(n)=xmag(i) 400 continue c write(*,*)nn,n do 500 i=1,n-1 do 500 j=i+1,n if(xmag(i).lt.xmag(j))goto 500 x=xmag(i) xmag(i)=xmag(j) xmag(j)=x x=aa(i) aa(i)=aa(j) aa(j)=x x=dd(i) dd(i)=dd(j) dd(j)=x 500 continue if(n.gt.400)n=400 ! control gsc number 400, star 300 end subroutine find(dir,f1,aa,dd,ee,col,nn,n,up_d,dn_d,up_a,dn_a,iaf) c find gsc c dir file: is cataloge of star c aa dd ee: 10 byte of gsc star c nn n : max array and now size of array c up.. dn..: 4 corner of chart, in hour & degree c iaf : R.A. indicator, for join 24~0 hour data real aa(nn),dd(nn) integer*2 ee(nn) character*1 col(nn),color,dir1 character dir*6,f1*8 integer*2 k512,mag byte u(512),v(10,51),w(10) equivalence (u,v),(u(511),k512),(alpha,w),(delta,w(5)),(mag,w(9)) c 1 2 3 4 5 6 7 8 9 10 c ------------------ ------------------ - --- -- -- c alpha delta c mag pe me c c~classify in first bit, 0~stellar 1~non-stellar c mag 7 bit, (mag-5)*10 c pe~ position error in ", pe=pe*10, if pe > 15 then pe=0, up 4 bit c me~ magnitude error in m, me=me*10, if me > 15 then me=0, low 4 bit open(1,file='/EOD/GSC/'//dir//f1,status='old', * access='direct',recl=512) nrec=0 10 nrec=nrec+1 read(1,rec=nrec,err=100)u do 20 i=1,51 do 30 j=1,10 30 w(j)=v(j,i) call swap4(w,8) if(w(9).eq.0)goto 100 c add color identfy color='V' if(dir1.eq.'S')color='B' if(alpha.lt.0.)then alpha=-alpha color='B' if(dir1.eq.'S')color='V' endif if(delta.gt.90.)then delta=delta-90. color='R' endif if(delta.gt.up_d.or.delta.lt.dn_d)goto 20 if(iaf) 40,50,60 c iaf=-1, ex. center in 1 hour, dn_a in 23 hour. 24 23...........3 2 1 40 if(alpha.gt.up_a.and.alpha.lt.dn_a)goto 20 if(alpha.gt.dn_a)alpha=alpha-24. ! value around 0 pi goto 70 c iaf=0 normal ex. .........8 7 6 5 ................. 50 if(alpha.gt.up_a.or.alpha.lt.dn_a)goto 20 goto 70 c iaf=1 ex, center in 23 hour, up_a in 1 hour. 23 22 21.........1 0 60 if(alpha.gt.up_a.and.alpha.lt.dn_a)goto 20 if(alpha.lt.up_a)alpha=alpha+24. ! around 2*pi c bug 96,12,11, fixed c70 n=n+1 c if(n.gt.nn)return 70 if(n+1.gt.nn)return n=n+1 aa(n)=alpha dd(n)=delta ee(n)=mag col(n)=color 20 continue goto 10 100 close (1) end function jindex(alpha,delta,dir,file) c given a alpha delta, output a dir & file name c in this file, the given star must be there. character dir*6,file*8 c 1991,1,5 real a1(9537),a2(9537),d1(9537),d2(9537) character*5 cc(24) integer*2 ii(24) data ii/ 593,1177,1728,2258,2780,3245,3651,4013, * 4293,4491,4614,4662,5259,5837,6411,6988, * 7522,8021,8463,8839,9133,9345,9489,9537/ data ip/0/ data cc/'N0000','N0730','N1500','N2230','N3000','N3730', * 'N4500','N5230','N6000','N6730','N7500','N8230', * 'S0000','S0730','S1500','S2230','S3000','S3730', * 'S4500','S5230','S6000','S6730','S7500','S8230'/ if(ip.ne.0)goto 99 c read index file of gsc only once c because slow operate at optical disk, so search file first at /usr/local open(1,file='/EOD/GSC/regions.bin',status='old',form='unformatted') read(1)a1,a2,d1,d2 close (1) call swap4(a1,9537*4) call swap4(a2,9537*4) call swap4(d1,9537*4) call swap4(d2,9537*4) 99 ip=1 jindex=0 do 20 i=1,9537 if(alpha.lt.a1(i).or.alpha.gt.a2(i))goto 20 dd1=d1(i) dd2=d2(i) if(delta.lt.0.)dd1=dd2 if(delta.lt.0.)dd2=d1(i) if(delta.lt.dd1.or.delta.gt.dd2)goto 20 goto 30 20 continue jindex=1 return 30 write(file(1:),"(i4.4,'.BAO')")i do 40 j=1,24 if(i.le.ii(j))goto 50 40 continue 50 dir=cc(j)//'/' end function inpic(ix,iy) c judge ix,iy in frame or not inpic=0 if(ix.lt.1.or.ix.gt.512)return if(iy.lt.1.or.iy.gt.512)return inpic=1 end subroutine shift_c(gx,gy,aa,dd,n,xycoef,alpha_c,delta_c) c n===30 real gx(n),gy(n),aa(n),dd(n) real cx(30),cy(30) real czr1(30),czr2(30) integer*2 w(60) integer cz1(30),cz2(30) real*8 xycoef(2,6),adcoef(2,6) pi=4.*atan(1.0) x1=pi/12. x2=pi/180. call xytoad(xycoef,adcoef) do 100 i=1,n xa=aa(i)*x1 xd=dd(i)*x2 call standc(alpha_c,delta_c,xa,xd,xi,xn) cx(i)=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3) 100 cy(i)=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3) c call pgbegin(0,'/xw',1,1) c call pgenv(0.,512.,512.,0.,1,0) c call pgpoint(n,cx,cy,5) c call pgsci(3) c call pgpoint(n,gx,gy,6) m=0. do 30 i=1,n r=500 cz1(i)=0 do 40 j=1,n rr=sqrt((cx(i)-gx(j))**2+(cy(i)-gy(j))**2) if(rr.lt.r)then r=rr cz1(i)=j czr1(i)=r endif 40 continue r=500 cz2(i)=0 do 50 j=1,n rr=sqrt((cx(i)-gx(j))**2+(cy(i)-gy(j))**2) if(rr.lt.r.and.rr.ne.czr1(i))then r=rr cz2(i)=j czr2(i)=r endif 50 continue if(cz1(i).ne.0)then m=m+1 w(m)=czr1(i)*10 endif if(cz2(i).ne.0)then m=m+1 w(m)=czr2(i)*10 endif c write(*,*)i,czr1(i),czr2(i) 30 continue call whitexblack(w,m,1,peak,s) peak=peak*0.1 s=3 c write(*,*)'peak: ',peak m=0 do 60 i=1,n if(cz1(i).eq.0.and.cz2(i).eq.0)goto 60 70 if(cz1(i).eq.0)then if(abs(czr2(i)-peak).gt.s)goto 60 j=cz2(i) c write(*,*)i,cx(i)-gx(j),cy(i)-gy(j) m=m+1 czr1(m)=cx(i)-gx(j) czr2(m)=cy(i)-gy(j) goto 60 endif if(cz2(i).eq.0)then if(abs(czr1(i)-peak).gt.s)goto 60 j=cz1(i) c write(*,*)i,cx(i)-gx(j),cy(i)-gy(j) m=m+1 czr1(m)=cx(i)-gx(j) czr2(m)=cy(i)-gy(j) goto 60 endif x1=abs(czr1(i)-peak) x2=abs(czr2(i)-peak) if(x1.gt.x2)cz1(i)=0 if(x2.gt.x1)cz2(i)=0 goto 70 60 continue c do 80 i=1,m c80 write(*,*)i,czr1(i),czr2(i) x1=xmedian(czr1,m) x2=xmedian(czr2,m) s=-sqrt(x1*x1+x2*x2)*xycoef(1,1)*10800/pi write(*,11)s,x1,x2 11 format(' Auto shift',f5.1,' arc_minutes, ',2f6.1) c call pgend alpha_c=alpha_c+1.*x1*xycoef(1,1) delta_c=delta_c+1.*x2*xycoef(1,1) end function xmedian(a,n) real a(n) do 10 i=1,n-1 do 10 j=i+1,n if(a(i).lt.a(j))goto 10 x=a(i) a(i)=a(j) a(j)=x 10 continue m=n/2 xmedian=a(m+1) if(m+m.ne.n)return xmedian=(a(m)+a(m+1))/2 end subroutine findstar(head,a_b,starx,stary,starv,nstar) parameter (n512=512) parameter (npos=5) PARAMETER (MAXEXP=10, MAXPAR=6, MAXBOX=13, MAXSKY=10000) PARAMETER (NOPT=20, NCMD=23, MAXPSF=207) parameter (MAXPIC=3*n512*MAXBOX) common /SIZE/ ncol,nrow common /EDGE/ n12,ix1,ix2,iy1,iy2,seeing real w1(MAXPIC) real w2(MAXBOX*MAXBOX*2) real opt(nopt) real a_b(n512*n512) real jnk(n512*n512) character*80 head(72) character*60 f1,f3 real starx(400),stary(400),starv(400) seeing=2.0 ! pixel ipos=indexpos(head,'NAXIS ') if(ipos.lt.72)read(head(ipos)(65:77),'(i5,i8)')n12,n49 if(n12.ne.0)then ipos=indexpos(head,'NAXIS1 ') read(head(ipos)(65:74),'(2i5)')ix1,ix2 ipos=indexpos(head,'NAXIS2 ') read(head(ipos)(65:74),'(2i5)')iy1,iy2 b_high=490000. if(n49.ne.0)b_high=n49 else ipos=indexpos(head,'TTIME ') b_high=28000. endif ipos=indexpos(head,'VOLT3 ') ! SEEING if(ipos.lt.72)read(head(ipos)(25:),*)seeingo if(seeingo.gt.3.)seeing=seeingo/1.67 ipos=indexpos(head,'NAXIS1 ') read(head(ipos)(15:),*)ncol ipos=indexpos(head,'NAXIS2 ') read(head(ipos)(15:),*)nrow c READ NOISE (ADU; 1 frame) = 29.50 / 1 c GAIN (e-/ADU; 1 frame) = 3.30 / 2 c LOW GOOD DATUM (in sigmas) = 5.00 / 3 c HIGH GOOD DATUM (in ADU) = 28000.00 / 4 c FWHM OF OBJECT = 3.37 / 5 c THRESHOLD (in sigmas) = 3.50 / 6 c LS (LOW SHARPNESS CUTOFF) = 0.40 / 7 c HS (HIGH SHARPNESS CUTOFF) = 1.20 / 8 c LR (LOW ROUNDNESS CUTOFF) = -1.15 / 9 c HR (HIGH ROUNDNESS CUTOFF) = 1.25 /10 c WATCH PROGRESS = -2.00 /11 c FITTING RADIUS = 3.00 /12 c PSF RADIUS = 6.00 /13 c VARIABLE PSF = 2.00 /14 c FRACTIONAL-PIXEL EXPANSION = 0.00 /15 c ANALYTIC MODEL PSF = 3.00 /16 c EXTRA PSF CLEANING PASSES = 3.00 /17 c USE SATURATED PSF STARS = 0.00 /18 c PERCENT ERROR (in %) = 0.75 /19 c PROFILE ERROR (in %) = 5.00 /20 opt( 1)= 29.50 opt( 2)= 3.5 opt( 3)= 5.6 opt( 4)=28800. opt( 5)= 3.37 opt( 6)= 3.50 opt( 7)= 0.40 opt( 8)= 1.20 opt( 9)= -1.15 opt(10)= 1.25 opt(11)= -2.00 opt(12)= 3.00 opt(13)= 6.00 opt(14)= 2.00 opt(15)= 0.00 opt(16)= 3.00 opt(17)= 3.00 opt(18)= 0.00 opt(19)= 0.75 opt(20)= 5.00 if(opt(4).gt.b_high)b_high=opt(4) opt(4)=b_high if(seeingo.gt.5..and.opt(5).lt.5.)then opt(7)=opt(7)-(seeingo-5.)*0.1 endif if(n12.eq.0)then n12=1 ix1=1 iy1=1 ix2=ncol iy2=nrow endif ix1=ix1+npos iy1=iy1+npos ix2=ix2-npos iy2=iy2-npos MAX = MAXPIC/3 MAXCOL = MAX/MAXBOX CALL FINDS(w1,w1(MAX+1),w1(2*MAX+1),w2,w2(MAXBOX*MAXBOX+1), * MAX, MAXBOX, MAXCOL, MAXSKY, OPT, NOPT,a_b,jnk, * starx,stary,starv,nstar) c open(3,file="1.coo",status='UNKNOWN') c do 330 i=1,nstar c330 WRITE (3,331)i,starx(i),stary(i),starv(i) c331 FORMAT(I6, 2F9.3, f9.0) c close(3) end SUBROUTINE FINDS(D, H, JCYLN, G, SKIP, MAX, MAXBOX, MAXCOL, . MAXSKY, OPT, NOPT, a_b,jnk, * starx,stary,starv,nstar) real starx(1),stary(1),starv(1) INTEGER MAX, MAXBOX, MAXCOL, MAXSKY, NOPT, n12,ix1,ix2,iy1,iy2 REAL D(MAXCOL,MAXBOX), H(MAXCOL,MAXBOX), DATA(2), OPT(NOPT) REAL G(MAXBOX,MAXBOX), AMAX1 INTEGER JCYLN(MAX) LOGICAL SKIP(MAXBOX,MAXBOX) CHARACTER LINE*5 REAL PIXELS, RADIUS, FWHM, SIGSQ, RSQ, RELERR, SKYLVL, TEMP REAL HMIN, LOBAD, HIBAD, P, DATUM, HEIGHT, DENOM, SGOP REAL SHARP, ROUND, SHRPLO, SHRPHI, RNDLO, RNDHI REAL SUMG, SUMGSQ, SUMGD, SUMD, SG, SGSQ, SGD, SD, WT, HX, HY REAL DGDX, SDGDX, SDGDXS, SDDGDX, SGDGDX, SUM2, SUM4 REAL XCEN, YCEN, DX, DY, PHPADU, READNS, SKYMOD INTEGER NHALF, NBOX, MIDDLE, LASTCL, LASTRO, NCOL, NROW, JSQ INTEGER ISTAT, NROWS, NSTAR INTEGER I, J, K, N, IX, IY, JX, JY, KX, LX, LY LOGICAL SATUR8 real a_b(512*512) real jnk(512*512) common /SIZE/ ncol,nrow common /EDGE/ n12,ix1,ix2,iy1,iy2,seeing HIBAD = OPT(4) FWHM = OPT(5) if(seeing.ne.0)FWHM=seeing SHRPLO= OPT(7) SHRPHI= OPT(8) RNDLO = OPT(9) RNDHI = OPT(10) C RADIUS=AMAX1(2.001, 0.637*FWHM) NHALF=MIN0((MAXBOX-1)/2, INT(RADIUS)) NBOX=2*NHALF+1 ! Length of the side of the subarray MIDDLE=NHALF+1 C LASTRO=NROW-NHALF LASTCL=NCOL-NHALF C SIGSQ=(FWHM/2.35482)**2 RADIUS=RADIUS**2 C SUMG=0.0 SUMGSQ=0.0 PIXELS=0.0 DO J=1,NBOX JSQ=(J-MIDDLE)**2 C DO I=1,NBOX RSQ=FLOAT((I-MIDDLE)**2+JSQ) G(I,J)=EXP(-0.5*RSQ/SIGSQ) IF (RSQ .LE. RADIUS) THEN SKIP(I,J)=.FALSE. SUMG=SUMG+G(I,J) SUMGSQ=SUMGSQ+G(I,J)**2 PIXELS=PIXELS+1.0 ELSE SKIP(I,J)=.TRUE. END IF END DO END DO DENOM=SUMGSQ-(SUMG**2)/PIXELS SGOP=SUMG/PIXELS C RELERR=1.0/DENOM RELERR=SQRT(RELERR) CALL SKY (D, H, JCYLN, MIN0(MAX,MAXSKY), OPT(1), HIBAD, READNS, . PHPADU, SKYMOD, IX, a_b) data(1)=1 data(2)=n12 IF ((DATA(1) .LT. 0.5) .OR. (DATA(2) .LT. 0.5)) RETURN READNS = OPT(1)**2*DATA(2)/DATA(1) PHPADU = OPT(2)*DATA(1) HMIN = SQRT(READNS + AMAX1(0.,SKYMOD)/PHPADU) LOBAD = 0.1*NINT(10.*(SKYMOD-OPT(3)*HMIN)) HMIN = 0.01*NINT(100.*OPT(6)*RELERR*HMIN) READNS = SQRT(READNS) C DO JY=1,MIDDLE-1 DO IX=1,NCOL D(IX,JY) = -1.1E38 END DO END DO C LX = 1 LY = 1 NROWS = NHALF CALL RDARAY (a_b,LX,LY,NCOL,NROWS,MAXCOL,D(1,MIDDLE),ISTAT) NROWS=1 C JY=0 2020 JY=JY+1 ! Increment image-row pointer IF (JY .GT. NROW) GO TO 2100 ! Have we reached the bottom? C DO J=1,NBOX IY = JY + (J - MIDDLE) C I = IY + NHALF C JCYLN(J) = MOD(I-1,NBOX) + 1 END DO C LY = JY+NHALF IF (LY .LE. NROW) THEN CALL RDARAY (a_b, LX, LY, NCOL, NROWS, MAXCOL, . D(1,JCYLN(NBOX)), ISTAT) ELSE K = JCYLN(NBOX) DO IX=1,NCOL D(IX,K) = -1.1E38 END DO END IF C DO 2050 JX=1,NCOL C SGD=0. SD=0. SGSQ=SUMGSQ SG=SUMG P=PIXELS C DO 2040 IX = JX-NHALF,JX+NHALF I = MIDDLE + (IX-JX) DO 2040 J=1,NBOX K = JCYLN(J) IF (SKIP(I,J)) GO TO 2040 IF ((IX .GE. 1) .AND. (IX .LE. NCOL)) THEN DATUM = D(IX,K) IF ((DATUM .GE. LOBAD).AND.(DATUM .LE. HIBAD)) THEN SGD = SGD+G(I,J)*DATUM SD = SD+DATUM GO TO 2040 END IF END IF SGSQ = SGSQ-G(I,J)**2 SG = SG-G(I,J) P = P-1. 2040 CONTINUE C IF (P .GT. 1.5) THEN IF (P .LT. PIXELS) THEN SGSQ = SGSQ-(SG**2)/P IF (SGSQ .NE. 0.) THEN SGD = (SGD-SG*SD/P)/SGSQ ELSE SGD = 0. END IF ELSE SGD = (SGD-SGOP*SD)/DENOM END IF ELSE SGD = 0. END IF H(JX,2) = SGD 2050 CONTINUE CALL WRARAY (jnk, LX, JY, NCOL, NROWS, MAXCOL, H(1,2), ISTAT) GO TO 2020 2100 CONTINUE SKIP(MIDDLE,MIDDLE) = .TRUE. 3000 CONTINUE C DO JY=1,MIDDLE-1 DO IX=1,NCOL D(IX,JY) = -1.1E38 H(IX,JY) = 0. END DO END DO C LX = 1 LY = 1 NROWS = NHALF CALL RDARAY (a_b, LX, LY, NCOL, NROWS, MAXCOL, . D(1,MIDDLE), ISTAT) CALL RDARAY (jnk, LX, LY, NCOL, NROWS, MAXCOL, . H(1,MIDDLE), ISTAT) NROWS = 1 C NSTAR = 0 JY = 0 3020 JY = JY+1 C IF (JY .GT. NROW) return C DO J=1,NBOX IY = JY+(J-MIDDLE)+NHALF JCYLN(J) = MOD(IY-1,NBOX) + 1 END DO C LY = JY+NHALF IF (LY .LE. NROW) THEN CALL RDARAY (a_b, LX, LY, NCOL, NROWS, MAXCOL, . D(1,JCYLN(NBOX)), ISTAT) CALL RDARAY (jnk, LX, LY, NCOL, NROWS, MAXCOL, . H(1,JCYLN(NBOX)), ISTAT) ELSE K = JCYLN(NBOX) DO IX=1,NCOL D(IX,K) = -1.1E38 H(IX,K) = 0.0 END DO END IF C JX = 1 3040 HEIGHT = H(JX,JCYLN(MIDDLE)) C IF (HEIGHT .LT. HMIN) GO TO 3300 DO 3051 IX=JX-NHALF,JX+NHALF IF ((IX .LT. 1) .OR. (IX .GT. NCOL)) GO TO 3051 I = MIDDLE + (IX-JX) DO 3050 J=1,NBOX K = JCYLN(J) IF (SKIP(I,J)) GO TO 3050 IF (HEIGHT .LT. H(IX,K)) GO TO 3300 3050 CONTINUE 3051 CONTINUE C SHARP=0. ROUND=0. DATUM=D(JX,JCYLN(MIDDLE)) IF ((DATUM .LT. LOBAD) .OR. (DATUM .GT. HIBAD)) GO TO 3068 P=0. DO 3061 IX=JX-NHALF,JX+NHALF IF ((IX .LT. 1) .OR. (IX .GT. NCOL)) GO TO 3061 I = MIDDLE + (JX-IX) TEMP = 0.0 DO 3060 J=1,NBOX K=JCYLN(J) IF ((IX.EQ.JX).AND.(K.EQ.MIDDLE)) GO TO 3060 DATUM=D(IX,K) IF ((DATUM .GE. LOBAD) .AND. (DATUM .LE. HIBAD)) THEN TEMP = TEMP+(DATUM-SKYMOD) P = P + 1. END IF 3060 CONTINUE SHARP = SHARP+TEMP 3061 CONTINUE C SHARP=(D(JX,JCYLN(MIDDLE))-SKYMOD-SHARP/P)/HEIGHT IF ((SHARP .LT. SHRPLO) .OR. (SHARP .GT. SHRPHI)) GO TO 3200 3068 CONTINUE C IF ((JX .LT. MIDDLE) .OR. (JX .GT. LASTCL) .OR. . (JY .LT. MIDDLE) .OR. (JY .GT. LASTRO)) THEN XCEN = REAL(JX) YCEN = REAL(JY) GO TO 3190 END IF C start SUM2 = 0. SUM4 = 0. DO I=0,NHALF DO J=1,NHALF c if (jy .ge. 77) type 6667, h(jx-i,jcyln(middle-j)), c . h(jx+i,jcyln(middle+j)), c . h(jx-j,jcyln(middle+i)), h(jx+j,jcyln(middle-i)) c6667 format (4f6.0) SUM2 = SUM2 + . H(JX-I,JCYLN(MIDDLE-J)) + H(JX+I,JCYLN(MIDDLE+J)) - . H(JX-J,JCYLN(MIDDLE+I)) - H(JX+J,JCYLN(MIDDLE-I)) SUM4 = SUM4 + . ABS(H(JX-I,JCYLN(MIDDLE-J))) + . ABS(H(JX+I,JCYLN(MIDDLE+J))) + . ABS(H(JX-J,JCYLN(MIDDLE+I))) + . ABS(H(JX+J,JCYLN(MIDDLE-I))) END DO END DO ROUND = 2.*SUM2/SUM4 IF ((ROUND .LT. RNDLO) .OR. (ROUND .GT. RNDHI)) GO TO 3200 IX = JX-MIDDLE C SUMGD=0.0 SUMGSQ=0.0 SUMG=0.0 SUMD=0.0 SDGDX=0.0 SDGDXS=0.0 SDDGDX=0.0 SGDGDX=0.0 P=0. N=0 DO 3073 I=1,NBOX SG=0. SD=0. KX = IX+I DO 3070 J=1,NBOX WT=FLOAT(MIDDLE-ABS(J-MIDDLE)) K=JCYLN(J) DATUM=D(KX,K) IF ((DATUM .GE. LOBAD) .AND. (DATUM .LE. HIBAD)) THEN SD=SD+(DATUM-SKYMOD)*WT SG=SG+G(I,J)*WT END IF 3070 CONTINUE IF (SG .GT. 0.0) THEN WT=FLOAT(MIDDLE-ABS(I-MIDDLE)) SUMGD=SUMGD+WT*SG*SD SUMGSQ=SUMGSQ+WT*SG**2 SUMG=SUMG+WT*SG SUMD=SUMD+WT*SD P=P+WT N=N+1 DGDX=SG*(MIDDLE-I) SDGDXS=SDGDXS+WT*DGDX**2 SDGDX=SDGDX+WT*DGDX SDDGDX=SDDGDX+WT*SD*DGDX SGDGDX=SGDGDX+WT*SG*DGDX END IF 3073 CONTINUE C IF (N .LE. 2) GO TO 3200 HX=(SUMGD-SUMG*SUMD/P)/(SUMGSQ-(SUMG**2)/P) C IF (HX .LE. 0.) GO TO 3200 C SKYLVL=(SUMD-HX*SUMG)/P DX=(SGDGDX-(SDDGDX-SDGDX*(HX*SUMG+SKYLVL*P)))/(HX*SDGDXS/SIGSQ) XCEN=JX+DX/(1.+ABS(DX)) C IF ((XCEN .LT. 0.5) .OR. (XCEN .GT. NCOL+0.5)) GO TO 3200 C SUMGD=0. SUMGSQ=0. SUMG=0. SUMD=0. SDGDX=0. SDGDXS=0. SDDGDX=0. SGDGDX=0. P=0. N=0 SATUR8 = .FALSE. DO 3078 J=1,NBOX K=JCYLN(J) SG=0. SD=0. DO 3076 I=1,NBOX WT=FLOAT(MIDDLE-ABS(I-MIDDLE)) KX=IX+I DATUM=D(KX,K) IF (DATUM .LE. HIBAD) THEN IF (DATUM .GE. LOBAD) THEN SD=SD+(DATUM-SKYMOD)*WT SG=SG+G(I,J)*WT END IF ELSE IF (.NOT. SKIP(I,J)) SATUR8 = .TRUE. END IF 3076 CONTINUE C IF (SG .GT. 0.0) THEN WT=FLOAT(MIDDLE-ABS(J-MIDDLE)) SUMGD=SUMGD+WT*SG*SD SUMGSQ=SUMGSQ+WT*SG**2 SUMG=SUMG+WT*SG SUMD=SUMD+WT*SD P=P+WT DGDX=SG*(MIDDLE-J) SDGDX=SDGDX+WT*DGDX SDGDXS=SDGDXS+WT*DGDX**2 SDDGDX=SDDGDX+WT*SD*DGDX SGDGDX=SGDGDX+WT*SG*DGDX N=N+1 END IF C 3078 CONTINUE C IF (N .LE. 2) GO TO 3200 HY=(SUMGD-SUMG*SUMD/P)/(SUMGSQ-(SUMG**2)/P) IF (HY .LE. 0.0) GO TO 3200 SKYLVL=(SUMD-HY*SUMG)/P DY=(SGDGDX-(SDDGDX-SDGDX*(HY*SUMG+SKYLVL*P)))/(HY*SDGDXS/SIGSQ) YCEN=JY+DY/(1.+ABS(DY)) IF ((YCEN .LT. 0.5) .OR. (YCEN .GT. NROW+0.5)) GO TO 3200 DY = 2.*(HX-HY)/(HX+HY) C IF (.NOT. SATUR8) THEN IF ((DY .LT. RNDLO) .OR. (DY .GT. RNDHI)) GO TO 3200 END IF C 3190 if(xcen.lt.ix1)goto 3200 if(xcen.gt.ix2)goto 3200 if(ycen.lt.iy1)goto 3200 if(ycen.gt.iy2)goto 3200 if(nstar.lt.400)NSTAR=NSTAR+1 ! 51228 starx(nstar)=XCEN stary(nstar)=YCEN starv(nstar)=(hx+hy)/2. 3200 CONTINUE C JX = JX+NHALF 3300 JX = JX+1 C IF (JX .LE. NCOL) GO TO 3040 GO TO 3020 END! SUBROUTINE SKY (D, S, INDEX, MAX, READNS, HIBAD, SKYMN, SKYMED, . SKYMOD, N, a_b) C IMPLICIT NONE INTEGER MAX C REAL S(MAX), D(MAX) INTEGER INDEX(MAX) C REAL READNS, HIBAD, SKYMN, SKYMED, SKYMOD, SKYSIG, SKYSKW INTEGER NCOL, NROW, ISTEP, LX, LY, NX, NY, IROW, I, N INTEGER ISTAT, IFIRST real a_b(512*512) common /SIZE/ ncol,nrow C ISTEP = NCOL*NROW/MAX+1 C LX = 1 NX = NCOL NY = 1 IFIRST = 0 N = 0 DO IROW=1,NROW LY = IROW CALL RDARAY (a_b, LX, LY, NX, NY, MAX, D, ISTAT) IF (ISTAT .NE. 0) RETURN IFIRST = IFIRST + 1 IF (IFIRST .GT. ISTEP) IFIRST = IFIRST - ISTEP I = IFIRST 1010 IF (ABS(D(I)) .LE. HIBAD) THEN N = N+1 S(N) = D(I) IF (N .EQ. MAX) GO TO 1100 I = I + ISTEP ELSE I = I+1 END IF IF (I .LE. NCOL) GO TO 1010 END DO C 1100 CONTINUE CALL QUICK (S, N, INDEX) CALL MMM (S, N, HIBAD, READNS, SKYMN, SKYMED, SKYMOD, SKYSIG, . SKYSKW) RETURN END! SUBROUTINE QUICK (DATUM, N, INDEX) IMPLICIT NONE INTEGER MAXSTK, N PARAMETER (MAXSTK=28) REAL DATUM(N) INTEGER INDEX(N), STKLO(MAXSTK), STKHI(MAXSTK) C REAL DKEY INTEGER I, HI, LO, NSTAK, LIMLO, LIMHI, IKEY C DO I=1,N INDEX(I)=I END DO C NSTAK=0 LIMLO=1 LIMHI=N C 100 DKEY=DATUM(LIMLO) IKEY=INDEX(LIMLO) C LO=LIMLO HI=LIMHI 101 CONTINUE C IF (LO .EQ. HI)GO TO 200 C IF (DATUM(HI) .LE. DKEY) GO TO 109 HI=HI-1 C GO TO 101 C 109 DATUM(LO)=DATUM(HI) INDEX(LO)=INDEX(HI) LO=LO+1 110 CONTINUE C IF (LO .EQ. HI) GO TO 200 C IF (DATUM(LO) .GE. DKEY) GO TO 119 C LO=LO+1 GO TO 110 C 119 DATUM(HI)=DATUM(LO) INDEX(HI)=INDEX(LO) HI=HI-1 C GO TO 101 C 200 CONTINUE C DATUM(LO)=DKEY INDEX(LO)=IKEY C IF (LIMHI-LO .GT. LO-LIMLO) GO TO 300 C IF (LO-LIMLO .LE. 1) GO TO 400 C IF (LIMHI-LO .GE. 2) GO TO 250 C LIMHI=LO-1 GO TO 100 C 250 CONTINUE C NSTAK=NSTAK+1 STKLO(NSTAK)=LIMLO STKHI(NSTAK)=LO-1 LIMLO=LO+1 GO TO 100 C 300 CONTINUE C IF (LIMHI-LO .LE. 1) GO TO 400 C IF (LO-LIMLO .GE. 2) GO TO 350 C LIMLO=LO+1 GO TO 100 C 350 CONTINUE C NSTAK=NSTAK+1 STKLO(NSTAK)=LO+1 STKHI(NSTAK)=LIMHI LIMHI=LO-1 GO TO 100 C 400 CONTINUE C IF (NSTAK .LE. 0) THEN RETURN ! Normal return END IF LIMLO=STKLO(NSTAK) LIMHI=STKHI(NSTAK) NSTAK=NSTAK-1 GO TO 100 C END! C SUBROUTINE MMM (SKY, NSKY, HIBAD, READNS, SKYMN, SKYMED, . SKYMOD, SIGMA, SKEW) C IMPLICIT NONE INTEGER NSKY REAL SKY(NSKY) C DOUBLE PRECISION DSQRT, DBLE REAL ALOG10, AMIN1, AMAX1 C DOUBLE PRECISION SUM,SUMSQ REAL CUT, CUT1, CUT2, DELTA, SKYMID, SKYMED, SKYMN, SKYMOD REAL SIGMA, SKEW, R, SIGN, HIBAD, CENTER, SIDE, READNS REAL DMOD, OLD, CLAMP INTEGER I, J, K, L, M INTEGER MINIMM, MAXIMM, NITER, ISTEP, MAXIT, MINSKY, JSTEP LOGICAL REDO DATA MAXIT / 30 /, MINSKY / 20 / C IF (NSKY .LE. 0) THEN GO TO 9900 END IF SKYMID=0.5*(SKY((NSKY+1)/2)+SKY(NSKY/2+1)) C SUM=0.D0 SUMSQ=0.D0 CUT1=AMIN1(SKYMID-SKY(1), SKY(NSKY)-SKYMID, HIBAD-SKYMID) C CUT2=SKYMID + CUT1 CUT1=SKYMID - CUT1 C MINIMM=0 DO 1010 I=1,NSKY IF (SKY(I) .LT. CUT1) THEN MINIMM=I GO TO 1010 END IF IF (SKY(I) .GT. CUT2) GO TO 1020 DELTA=SKY(I)-SKYMID SUM=SUM+DELTA SUMSQ=SUMSQ+DELTA**2 MAXIMM=I 1010 CONTINUE C 1020 CONTINUE SKYMED=0.5*(SKY((MINIMM+MAXIMM+1)/2)+SKY((MINIMM+MAXIMM)/2+1)) SKYMN=SUM/DBLE(MAXIMM-MINIMM) SIGMA=DSQRT(SUMSQ/DBLE(MAXIMM-MINIMM)-SKYMN**2) SKYMN=SKYMN+SKYMID C SKYMOD=SKYMN IF (SKYMED .LT. SKYMN) SKYMOD=3.*SKYMED-2.*SKYMN C NITER=0 OLD = 0. CLAMP = 1. 2000 NITER=NITER+1 IF ((NITER .GT. MAXIT) .OR. (MAXIMM-MINIMM .LT. MINSKY)) THEN GO TO 9900 END IF C R=ALOG10(FLOAT(MAXIMM-MINIMM)) R=AMAX1(2., (-.1042*R+1.1695)*R+.8895) C CUT=R*SIGMA+0.5*ABS(SKYMN-SKYMOD) CUT=AMAX1(1.5,CUT) CUT1=SKYMOD-CUT CUT2=SKYMOD+CUT C REDO=.FALSE. C ISTEP=INT(SIGN(1.0001, CUT1-SKY(MINIMM+1))) JSTEP=(ISTEP+1)/2 C IF (ISTEP .GT. 0) GO TO 2120 2100 IF ((ISTEP .LT. 0) .AND. (MINIMM .LE. 0)) GO TO 2150 C IF ((SKY(MINIMM) .LE. CUT1) .AND. (SKY(MINIMM+1) .GE. CUT1)) . GO TO 2150 C 2120 CONTINUE DELTA=SKY(MINIMM+JSTEP)-SKYMID SUM=SUM-REAL(ISTEP)*DELTA SUMSQ=SUMSQ-REAL(ISTEP)*DELTA**2 MINIMM=MINIMM+ISTEP REDO=.TRUE. ! A change has occured GO TO 2100 C 2150 CONTINUE C ISTEP=INT(SIGN(1.0001, CUT2-SKY(MAXIMM))) JSTEP=(ISTEP+1)/2 C IF (ISTEP .LT. 0) GO TO 2220 2200 IF ((ISTEP .GT. 0) .AND. (MAXIMM .GE. NSKY)) GO TO 2250 C IF ((SKY(MAXIMM) .LE. CUT2) .AND. (SKY(MAXIMM+1) .GE. CUT2)) . GO TO 2250 C 2220 DELTA=SKY(MAXIMM+JSTEP)-SKYMID SUM=SUM+REAL(ISTEP)*DELTA SUMSQ=SUMSQ+REAL(ISTEP)*DELTA**2 MAXIMM=MAXIMM+ISTEP REDO=.TRUE. ! A change has occured GO TO 2200 C 2250 CONTINUE C SKYMN=SUM/DBLE(MAXIMM-MINIMM) SIGMA=DSQRT(SUMSQ/DBLE(MAXIMM-MINIMM)-SKYMN**2) SKYMN=SKYMN+SKYMID C SKYMED=0.0 CENTER = REAL(MINIMM+1 + MAXIMM)/2. SIDE = REAL(NINT(0.2*REAL(MAXIMM-MINIMM)))/2. + 0.25 J = NINT(CENTER-SIDE) K = NINT(CENTER+SIDE) L = NINT(CENTER-0.25) M = NINT(CENTER+0.25) R = 0.25*READNS 2305 CONTINUE IF ((J .GT. 1) .AND. (K .LT. NSKY) .AND. ( . (SKY(L)-SKY(J) .LT. R) .OR. (SKY(K)-SKY(M) .LT. R) )) THEN J = J-1 K = K+1 GO TO 2305 END IF C DO 2310 I=J,K 2310 SKYMED=SKYMED+SKY(I) C SKYMED=SKYMED/REAL(K-J+1) IF (SKYMED .LT. SKYMN) THEN DMOD=3.*SKYMED-2.*SKYMN - SKYMOD ELSE DMOD=SKYMN - SKYMOD END IF C IF (DMOD*OLD .LT. 0.) CLAMP = 0.5*CLAMP SKYMOD = SKYMOD + CLAMP*DMOD OLD = DMOD IF (REDO) GO TO 2000 C SKEW=(SKYMN-SKYMOD)/AMAX1(1., SIGMA) NSKY=MAXIMM-MINIMM RETURN C 9900 SIGMA=-1.0 SKEW=0.0 RETURN C END! subroutine rdaray (id, lx, ly, mx, my, nx, func, ier) real func(nx,*),id(512*512) common /SIZE/ ncol,nrow mx = lx+mx-1 my = ly+my-1 lx = max(1,lx) ly = max(1,ly) mx = min(ncol,mx) my = min(nrow,my) my = my-ly+1 do 10 j=1,my jy = ly+j-1 10 call imgs2r (id, func(1,j), lx, mx, jy, jy, ier, ncol) mx = mx-lx+1 end subroutine wraray (id, lx, ly, mx, my, maxx, func, ier) real id(512*512) real func(maxx,*) common /size/ ncol, nrow mx = lx+mx-1 my = ly+my-1 lx = max(1,lx) ly = max(1,ly) mx = min(ncol, mx) my = min(nrow, my) nx = mx-lx+1 ny = my-ly+1 do 10 j=1,ny jy = ly+j-1 10 call imps2r (id, func(1,j), lx, mx, jy, jy, ier, ncol) mx = nx my = ny end subroutine imps2r(id,buf,i1,i2,j1,j2,ier,ncol) real id(ncol,*),buf(1) k=0 do 10 j=j1,j2 do 10 i=i1,i2 k=k+1 10 id(i,j)=buf(k) ier=0 end subroutine imgs2r(id,buf,i1,i2,j1,j2,ier,ncol) real id(ncol,*),buf(1) k=0 do 10 j=j1,j2 do 10 i=i1,i2 k=k+1 10 buf(k)=id(i,j) ier=0 end