C+ C NAME: C find_theta_rot C PURPOSE: C find the angle of intersection between a line on the ccd of constant thetax, and a line of longitude on the globe C CATEGORY: C Data processing C CALLING SEQUENCE: C call find_theta_rot(mode,ic,q,ra,dec,dr,thetax,thetarot,badangle) C INPUTS: C mode camera binning style; mode=N -> NxN binning C ic camera number(0-4); 0 -> TMO camera C q(4) quaternion for coordinate rotation C ra/dec coordinates of star C thetarot sum of intersection angle, to be averaged C xra/ydec ra/dec to be used as cartesian coordinates in finding angle of intersection C dr number of pixels cra/cdec is from r0 C PROCEDURE: C Find two points on ccd of constant thetax, rotate those coordinates so that cra, cdec become [0,0] C then create a line from those two points and find the angle between it and y-axis (dec-axis) C Final thetarot should be clockwise angle of rotation for fishtail to be pointed north C MODIFICATION HISTORY: C Aug, 2003 Aaron Smith (CASS) C Aug, 2004 Aaron Smith (CASS) C- subroutine find_theta_rot(mode,ic,q,ra,dec,dr,thetax,thetarot,image1,image2,image3, & imagecnt1,imagecnt2,imagecnt3,npolarimage,spolarimage,npolarimagecnt,spolarimagecnt) implicit none real*8 aq, bq, cq, xccd, yccd, thetarot, xra(2), ydec(2) real*8 rq, dr, thetax, thetay, cx, cy, cz, crsq, ra, dec, arg real*8 r0, tempx, phi, rlat, rccd, q(4), cra1, cra2, cra3 real*8 alfa, beta, gamma, xtri, ytri, hyptri, slope integer*4 ii, jj, mode, ic, badangle, xarg, yarg real*4 image1(120,120), image2(120,120), image3(120,120) real*4 npolarimage(160,160), spolarimage(160,160) integer*4 imagecnt1(120,120), imagecnt2(120,120), imagecnt3(120,120) integer*4 spolarimagecnt(160,160), npolarimagecnt(160,160) integer*2 k, level, zz, ibad if(ic.eq.1) r0=1193.5D0 if(ic.eq.2) r0=1188.5D0 if(ic.eq.3) r0=1198.5D0 badangle=0 do ii=1,2 if(ii.eq.1) rq=dr+6.D0 if(ii.eq.2) rq=dr-6.D0 if(rq.gt.32.4D0.or.rq.lt.-27.D0) then badangle=1 !!dont use a line that is outside the ROI goto 976 endif !! quadratic coeffitients from Andy's July 2001 memo aq = 0.1664D0 - (0.00001*thetax*thetax) bq = (0.000036D0*thetax*thetax) - 18.0226D0 cq = (0.000936D0*thetax*thetax) - rq thetay = ( -bq-sqrt((bq*bq)-(4.D0*aq*cq)) ) / ( 2.D0*aq ) cx = dsind(thetax) cy = dsind(thetay) crsq = (cx*cx)+(cy*cy) if(crsq.gt.0.D0) then phi = datan2d(cy,cx) else phi = 0.D0 endif cz = dsqrt(1.D0-(crsq)) rlat = dasind(cz) call rotateq(q, phi, rlat) !!rotate to standard system using original qq quaternion xra(ii)=phi ydec(ii)=rlat enddo alfa=ra !!line up with x-axis beta=-dec !!rotate about new y-axis so centroid now at [0,0] gamma=0.D0 do ii=1,2 phi=xra(ii) rlat=ydec(ii) call rotate8(alfa,beta,gamma,phi,rlat) if( phi.gt.180.D0 ) phi=phi-360.D0 if( phi.lt.-180.D0 ) phi=phi+360.D0 xra(ii)=phi ydec(ii)=rlat enddo !!Find angle between y-axis (dec-axis) and line of constant thetax xtri=xra(1)-xra(2) ytri=dabs(ydec(1)-ydec(2)) hyptri=dsqrt( xtri*xtri + ytri*ytri ) thetarot=dacosd(ytri/hyptri) slope=(ydec(1)-ydec(2))/(xra(1)-xra(2)) ccc Final thetarot should be clockwise angle of rotation for fishtail (bunny ears) ccc to be pointed north...note that RA/DEC form left handed coordinates when looking out ccc from Earth! if(slope.gt.0.D0)then if(xra(1).lt.xra(2)) thetarot=thetarot if(xra(1).gt.xra(2)) thetarot=180.D0+thetarot else if(xra(1).lt.xra(2)) thetarot=180.D0-thetarot !270.D0-thetarot if(xra(1).gt.xra(2)) thetarot=360.D0-thetarot endif !! 0 < thetarot < 360 if(thetarot.lt.0.D0) thetarot=360.D0+thetarot cccccccccccccccccccccccccccc Fill Coarse Rotator Maps ccccccccccccccccccccccccccccccccccccccc cra1=60.D0 cra2=180.D0 cra3=300.D0 c c Make first RA/dec plot (0-120 in ra) c ibad=0 xarg = nint( -(ra-cra1) + 60.5D0 ) yarg = nint( dec + 60.5D0 ) if(xarg.lt.1.or.xarg.gt.120.or.yarg.lt.1.or.yarg.gt.120)ibad=1 if(ibad.eq.0)then !!At the border between 0 and 360 degrees averaging gives bad values, !!so only use one of the values if(imagecnt1(xarg,yarg).ge.1)then arg=image1(xarg,yarg)/float(imagecnt1(xarg,yarg)) if(arg.gt.355.00.and.thetarot.lt.180.00) go to 976 if(arg.lt.5.00.and.thetarot.gt.180.00) go to 976 endif image1(xarg,yarg) = image1(xarg,yarg) + thetarot imagecnt1(xarg,yarg) = imagecnt1(xarg,yarg) + 1 endif c c Make second RA/dec plot (120-240 in ra) c ibad=0 xarg = nint( -(ra-cra2) + 60.5D0 ) yarg = nint( dec + 60.5D0 ) if(xarg.lt.1.or.xarg.gt.120.or.yarg.lt.1.or.yarg.gt.120)ibad=1 if(ibad.eq.0)then !!At the border between 0 and 360 degrees averaging gives bad values, !!so only use one of the values if(imagecnt2(xarg,yarg).ge.1)then arg=image2(xarg,yarg)/float(imagecnt2(xarg,yarg)) if(arg.gt.355.00.and.thetarot.lt.180.00) go to 976 if(arg.lt.5.00.and.thetarot.gt.180.00) go to 976 endif image2(xarg,yarg) = image2(xarg,yarg) + thetarot imagecnt2(xarg,yarg) = imagecnt2(xarg,yarg) + 1 endif c c Make third RA/dec plot (240-360 in ra) c ibad=0 xarg = nint( -(ra-cra3) + 60.5D0 ) yarg = nint( dec + 60.5D0 ) if(xarg.lt.1.or.xarg.gt.120.or.yarg.lt.1.or.yarg.gt.120)ibad=1 if(ibad.eq.0)then !!At the border between 0 and 360 degrees averaging gives bad values, !!so only use one of the values if(imagecnt3(xarg,yarg).ge.1)then arg=image3(xarg,yarg)/float(imagecnt3(xarg,yarg)) if(arg.gt.355.00.and.thetarot.lt.180.00) go to 976 if(arg.lt.5.00.and.thetarot.gt.180.00) go to 976 endif image3(xarg,yarg) = image3(xarg,yarg) + thetarot imagecnt3(xarg,yarg) = imagecnt3(xarg,yarg) + 1 endif c c Make North Polar plot c if(dec.gt.33.D0)then ibad=0 xarg = nint( 2.D0*(90.D0-dec)*cosd(-ra) + 80.5D0) yarg = nint( 2.D0*(90.D0-dec)*sind(-ra) + 80.5D0) if(xarg.lt.1.or.xarg.gt.160.or.yarg.lt.1.or.yarg.gt.160)ibad=1 if(ibad.eq.0)then !! Convert eq angle to north polar angle thetarot=360.D0-(thetarot+90.D0-ra) if(thetarot.gt.360.0) then thetarot=thetarot-360.D0 elseif(thetarot.lt.0.0)then thetarot=thetarot+360.D0 endif !!At the border between 0 and 360 degrees averaging gives bad values, !!so only use one of the values if(npolarimagecnt(xarg,yarg).ge.1)then arg=npolarimage(xarg,yarg)/float(npolarimagecnt(xarg,yarg)) if(arg.gt.355.00.and.thetarot.lt.180.00) go to 976 if(arg.lt.5.00.and.thetarot.gt.180.00) go to 976 endif npolarimage(xarg,yarg) = npolarimage(xarg,yarg) + thetarot npolarimagecnt(xarg,yarg) = npolarimagecnt(xarg,yarg) + 1 endif endif c c Make South Polar plot c if(dec.lt.-33.D0)then ibad=0 xarg = nint( -2.D0*(90.D0+dec)*cosd(ra) + 80.5D0) yarg = nint( -2.D0*(90.D0+dec)*sind(ra) + 80.5D0) if(xarg.lt.1.or.xarg.gt.160.or.yarg.lt.1.or.yarg.gt.160)ibad=1 if(ibad.eq.0)then !! Convert eq angle to south polar angle thetarot=360.D0-(thetarot+ra-270.D0) if(thetarot.gt.360.0) then thetarot=thetarot-360.D0 elseif(thetarot.lt.0.0)then thetarot=thetarot+360.D0 endif !!At the border between 0 and 360 degrees averaging gives bad values, !!so only use one of the values if(spolarimagecnt(xarg,yarg).ge.1)then arg=spolarimage(xarg,yarg)/float(spolarimagecnt(xarg,yarg)) if(arg.gt.355.00.and.thetarot.lt.180.00) go to 976 if(arg.lt.5.00.and.thetarot.gt.180.00) go to 976 endif spolarimage(xarg,yarg) = spolarimage(xarg,yarg) + thetarot spolarimagecnt(xarg,yarg) = spolarimagecnt(xarg,yarg) + 1 endif endif 976 return end