Index: trunk/tools/eam/pscoords/distortion.sh =================================================================== --- trunk/tools/eam/pscoords/distortion.sh (revision 41267) +++ trunk/tools/eam/pscoords/distortion.sh (revision 41267) @@ -0,0 +1,681 @@ +#!/usr/bin/env mana + +if (not($?MAKEPLOT)) set MAKEPLOT = 0 +if (not($?TEMPDIR)) + $TEMPDIR = . +end + +$TIMEFORMAT = mjd +$TIMEREF = 2000/01/01,00:00:00 + +# GPC1 / GPC2 sizes +$NXch = 4846 +$NYch = 4868 + +macro choose.smfs + if ($0 != 2) + echo "USAGE: choose.smfs (version)" + break + end + + local i version + $version = $1 + + mkdir $version + $TEMPDIR = $version + + load.pscoords.corners + + file smflist.$version.txt isFound + if (not($isFound)) + echo "invalid version $version : smflist.$version.txt not found" + break + end + + list names -x "cat smflist.$version.txt" + for i 1 $names:n + list word -split $names:$i + $expname = $word:0 + + $nebpath = $word:1\.smf + + ## NOTE: assignment from backticks does not set $STATUS to FALSE if command exits with non-zero status + ## break -auto off + ## $diskpath = `neb-locate -p $nebpath` + ## $myStatus = $STATUS + ## echo $STATUS $myStatus "...$diskpath..." + ## break -auto on + ## if (not($myStatus)) + ## echo "skipping $nebpath" + ## continue + ## end + + $diskpath = `neb-locate -p $nebpath` + if ("$diskpath" == "") + echo "skipping $nebpath" + next + end + + echo $diskpath $expname + compare.smf.pscoords $diskpath $version/$expname + end +end + +macro compare.smf.pscoords + if ($0 != 3) + echo "USAGE: compare.smf.pscoords (smffile) (outroot)" + break + end + + local smffile outroot + $smffile = $1 + $outroot = $2 + + make.fpa.centers $smffile + + set Lo_um = +10 * Lo_pix + set Mo_um = -10 * Mo_pix + + foreach ver 00 10 01 11 + set L$ver\_um = +10 * L$ver\_pix + set M$ver\_um = -10 * M$ver\_pix + end + + # compare IPP prediction with mechanical model: + + # fit FP positions + fit2d Lo_mat Mo_mat Lo_um 1 + applyfit2d Lo_mat Mo_mat Lo_fit + set dLoR = Lo_um - Lo_fit + $C00_L = $CX0Y0 + $C10_L = $CX1Y0 + $C01_L = $CX0Y1 + + foreach type 00 01 10 11 + applyfit2d L$type\_mat M$type\_mat L$type\_fit + end + + fit2d Lo_mat Mo_mat Mo_um 1 + applyfit2d Lo_mat Mo_mat Mo_fit + set dMoR = Mo_um - Mo_fit + $C00_M = $CX0Y0 + $C10_M = $CX1Y0 + $C01_M = $CX0Y1 + + foreach type 00 01 10 11 + applyfit2d L$type\_mat M$type\_mat M$type\_fit + end + + # C10_L C01_L + # C10_M C01_M + + $Scale1 = sqrt($C10_L^2 + $C10_M^2) + $Scale2 = sqrt($C01_L^2 + $C01_M^2) + + $Rot1 = datan2($C10_M, $C10_L) + $Rot2 = datan2($C01_L, $C01_M) + + echo Scale1, Scale2: $Scale1, $Scale2 + echo Rot1, Rot2: $Rot1, $Rot2 + + vstat -q dLoR; $SIGMA_Lo = $SIGMA + vstat -q dMoR; $SIGMA_Mo = $SIGMA + + $BOXST = -tickpad 0.1 -xpad 3.5 -labelpadx 3 + + foreach type 00 01 10 11 + # the vector diff for a corner is a measurement of the rotation + set dL$type = L$type\_um - L$type\_fit + set dM$type = M$type\_um - M$type\_fit + + # rotation vector + set drL$type = dL$type - dLoR + set drM$type = dM$type - dMoR + set dr$type = sqrt(drL$type^2 + drM$type^2) + + # XXX pre-calculate this: + set dRL$type = L$type\_mat - Lo_mat + set dRM$type = M$type\_mat - Mo_mat + set R$type = sqrt(dRL$type^2 + dRM$type^2) + + set phi$type = 1000*asin(dr$type / R$type) +# set flipSign = ((drL$type < 0) && (otaID_mat < 40)) || ((drL$type > 0) && (otaID_mat >=40)) + set flipSign = ((drL00 < 0) && (otaID_mat < 40)) || ((drL00 > 0) && (otaID_mat >=40)) + set phi$type = flipSign ? phi$type : 0 - phi$type + end + + set phi_ave = 0.25*(phi00 + phi01 + phi10 + phi11) + set phi_std = sqrt((4/3)*(0.25*(phi00^2 + phi01^2 + phi10^2 + phi11^2) - phi_ave^2)) + + if ($MAKEPLOT) + clear -s -n 0 + resize 1800 1800; + label -fn helvetica 18 + + section a00 0.0 0.0 0.5 0.5 + section a10 0.5 0.0 0.5 0.5 + section a01 0.0 0.5 0.5 0.5 + section a11 0.5 0.5 0.5 0.5 + + sprint lineL "dL stdev = %5.1f" $SIGMA_Lo + sprint lineM "dM stdev = %5.1f" $SIGMA_Mo + sprint lineR "Rot1: %6.2f, Rot2: %6.2f" $Rot1 $Rot2 + sprint lineS "Sca1: %6.2f, Sca2: %6.2f" {100*$Scale1} {100*$Scale2} + + section a00; lim Lo_mat dLoR; label -y dL -x L; box $BOXST; plot Lo_mat dLoR; label +x "$lineL, $lineM" + section a11; lim Mo_mat dMoR; label -y dM -x M; box $BOXST; plot Mo_mat dMoR; label +x "$lineR" + section a01; lim Lo_mat dMoR; label -y dM -x L; box $BOXST; plot Lo_mat dMoR; label +x "$lineS" + section a10; lim Mo_mat dLoR; label -y dL -x M; box $BOXST; plot Mo_mat dLoR + antialias 0.4 + + png -name $outroot.dLdM.png + + clear -s -n 1 + lim Lo_mat Mo_mat; clear; resize 1800 1800; + box; plot Lo_mat Mo_mat -pt ocir -c blue -sz 0.5; needles Lo_mat Mo_mat dLoR dMoR -c red -scale 100 + + antialias 0.4 + png -name $outroot.LM.png + +# lim -n 0 otaID_mat -10 10; clear; box; plot otaID_mat phi00 -pt ocir; plot -c indigo otaID_mat phi11 -pt circ; plot -c blue otaID_mat phi10 -pt plus; plot -c red otaID_mat phi01 -pt x + lim -n 0 otaID_mat -10 10; clear; box; plot otaID_mat phi_ave -dy phi_std -pt cir -c blue + +# NOTE: this test is for o8635g0206o.1479168.cm.2139537.smf + data rot.psc.dat + read otaID_psc:int 1 phi_psc 2 +# lim -n 1 phi_psc -5 5; clear; box; plot phi_psc phi00 -pt ocir; plot -c indigo phi_psc phi11 -pt circ; plot -c blue phi_psc phi10 -pt plus; plot -c red phi_psc phi01 -pt x + lim -n 2 phi_psc -5 5; clear; box; plot phi_psc phi_ave -dy phi_std -pt cir -c blue + end + + mcreate phu 0 0 + + keyword phu ROT1 -wf $Rot1 + keyword phu ROT2 -wf $Rot2 + keyword phu SCALE1 -wf $Scale1 + keyword phu SCALE2 -wf $Scale2 + + keyword phu ROTANGLE -wf $ROTANGLE + keyword phu POSANGLE -wf $POSANGLE + keyword phu ALT -wf $ALT + keyword phu AZ -wf $AZ + keyword phu FILTERID -w $FILTERID + keyword phu CRVAL1 -wf $CRVAL1_PHU + keyword phu CRVAL2 -wf $CRVAL2_PHU + keyword phu DETTEM -wf $DETTEM + keyword phu M2Z -wf $M2Z + keyword phu MJDOBS -wf $MJDOBS + + wd phu $outroot.psc.fits + write -fits DATA -append $outroot.psc.fits otaID_mat Lo_mat Mo_mat Lo_um Mo_um dLoR dMoR phi_ave phi_std +end + +# I need to calculate the coordinates of the center pixels in the focal plane. +# pscoords wants (dx, dy, theta) in (um, um, radian) + +macro make.fpa.centers + if ($0 != 2) + echo "USAGE: make.fpa.centers (smffile)" + break + end + + local smffile chipname i j + $smffile = $1 + + rd -head smfphu $smffile + + keyword smfphu ROT ROTANGLE + keyword smfphu POSANGLE POSANGLE + keyword smfphu ALT ALT + keyword smfphu AZ AZ + keyword smfphu FILTERID FILTERID + keyword smfphu CRVAL1 CRVAL1_PHU + keyword smfphu CRVAL2 CRVAL2_PHU + keyword smfphu DETTEM DETTEM + keyword smfphu M2Z M2Z + keyword smfphu MJD-OBS MJDOBS + + delete -q Lo_pix Mo_pix L00_pix M00_pix L10_pix M10_pix L01_pix M01_pix L11_pix M11_pix + + # generate a list of the OTA IDs in order + exec ./get.chip.IDs $smffile > $TEMPDIR/tmp.smf.dat + data $TEMPDIR/tmp.smf.dat + read Nastro 1 otaID_all:int 2 + + # we need to track the chip sequence number in the file, but then + # only capture a subset of the entries (Nastro >= 1) + + set keep = (Nastro >= 1) + subset otaID_smf = otaID_all where keep + + for i 0 otaID_all[] + if (not(keep[$i])) continue + make.chip.points $smffile $i + concat fpa_x[0] Lo_pix + concat fpa_y[0] Mo_pix + concat fpa_x[1] L00_pix + concat fpa_y[1] M00_pix + concat fpa_x[2] L10_pix + concat fpa_y[2] M10_pix + concat fpa_x[3] L11_pix + concat fpa_y[3] M11_pix + concat fpa_x[4] L01_pix + concat fpa_y[4] M01_pix + end + + match.pscoords.to.smf + + echo "done with make.fpa.centers" +end + +######## pscoords ops ########### + +# generate a file with all chips (00 - 77) and entries for center and all 4 corners (00, 10, 01, 11) +macro make.pscoords.corners + + exec rm -f pscoords.corners.dat + + local i j + + delete -q otaID + for i 0 8 + for j 0 8 + concat $i\$j otaID + end + end + + set Xc = zero(otaID) + $NXch/2 + set Yc = zero(otaID) + $NYch/2 + set X0 = zero(otaID) + 0 ; ## use this for Y0 as well + set X1 = zero(otaID) + $NXch + set Y1 = zero(otaID) + $NYch + + write $TEMPDIR/tmp.ota.CC.dat otaID Xc Yc + write $TEMPDIR/tmp.ota.00.dat otaID X0 X0 + write $TEMPDIR/tmp.ota.01.dat otaID X0 Y1 + write $TEMPDIR/tmp.ota.10.dat otaID X1 X0 + write $TEMPDIR/tmp.ota.11.dat otaID X1 Y1 + + foreach ver CC 00 01 10 11 + # chip centers, with offsets + exec pscoords -nooff in=ota out=fp < $TEMPDIR/tmp.ota.$ver.dat > $TEMPDIR/tmp.fp.$ver.dat + end + + write $TEMPDIR/tmp.ota.IDs.dat otaID + + exec paste $TEMPDIR/tmp.ota.IDs.dat $TEMPDIR/tmp.fp.CC.dat | paste - $TEMPDIR/tmp.fp.00.dat | paste - $TEMPDIR/tmp.fp.10.dat | paste - $TEMPDIR/tmp.fp.01.dat | paste - $TEMPDIR/tmp.fp.11.dat > pscoords.corners.dat +end + +macro load.pscoords.corners + + # read the pscoords data into vectors with ext _ps + data pscoords.corners.dat + read otaID_ps:int 1 Lo_ps 2 Mo_ps 3 L00_ps 4 M00_ps 5 L10_ps 6 M10_ps 7 L01_ps 8 M01_ps 9 L11_ps 10 M11_ps 11 +end + +macro match.pscoords.to.smf + # pscoords vectors have extension _ps + # generate a matched set by joining on the otaID values + + join otaID_smf otaID_ps + + foreach field otaID Lo Mo L00 M00 L10 M10 L01 M01 L11 M11 + reindex $field\_mat = $field\_ps using index2 + end +end + +# result: (dL, dM) for chips with Nastro > 0 +macro make.pscoords.rottest + + # confirm that we understand the structure in libpscoords.c + data pscoords.gpc1.v3.dat + read dLx 1 dMx 2 rotrad 3 otaID:int 5 + set Nastro = zero(otaID) + 1 + + # Can I generate the FP corner map from pscoords? + + exec rm -f $TEMPDIR/tmp.ota.CC.dat $TEMPDIR/tmp.ota.00.dat $TEMPDIR/tmp.ota.01.dat $TEMPDIR/tmp.ota.10.dat $TEMPDIR/tmp.ota.11.dat + + delete -q otaIDx Xc Yc X00 X01 X10 X11 Y00 Y01 Y10 Y11 + for i 0 otaID[] + if (Nastro[$i] < 1) continue + concat otaID[$i] otaIDx + end + + set Xc = zero(otaIDx) + $NXch/2 + set Yc = zero(otaIDx) + $NYch/2 + set X0 = zero(otaIDx) + 0 ; ## use this for Y0 as well + set X1 = zero(otaIDx) + $NXch + set Y1 = zero(otaIDx) + $NYch + + write $TEMPDIR/tmp.ota.CC.dat otaIDx Xc Yc + write $TEMPDIR/tmp.ota.00.dat otaIDx X0 X0 + write $TEMPDIR/tmp.ota.01.dat otaIDx X0 Y1 + write $TEMPDIR/tmp.ota.10.dat otaIDx X1 X0 + write $TEMPDIR/tmp.ota.11.dat otaIDx X1 Y1 + + # chip centers, with offsets + exec pscoords -dooff in=ota out=fp < $TEMPDIR/tmp.ota.CC.dat > $TEMPDIR/tmp.fp.dat + data $TEMPDIR/tmp.fp.dat + read Lx 1 Mx 2 + + # chip centers, sans offsets + exec pscoords -nooff in=ota out=fp < $TEMPDIR/tmp.ota.CC.dat > $TEMPDIR/tmp.fp.dat + data $TEMPDIR/tmp.fp.dat + read Lo 1 Mo 2 + + set dL = Lx - Lo + set dM = Mx - Mo + + foreach type 00 01 10 11 + # chip corner, with offsets + exec pscoords -dooff in=ota out=fp < $TEMPDIR/tmp.ota.$type.dat > $TEMPDIR/tmp.fp.dat + data $TEMPDIR/tmp.fp.dat + read Ltmpx 1 Mtmpx 2 + + # chip corner, sans offsets + exec pscoords -nooff in=ota out=fp < $TEMPDIR/tmp.ota.$type.dat > $TEMPDIR/tmp.fp.dat + data $TEMPDIR/tmp.fp.dat + read Ltmpo 1 Mtmpo 2 + + # the vector diff for a corner is a measurement of the rotation + set dL$type = Ltmpx - Ltmpo + set dM$type = Mtmpx - Mtmpo + + set dRL$type = Ltmpo - Lo + set dRM$type = Mtmpo - Mo + + set drL$type = dL$type - dL + set drM$type = dM$type - dM + + set dr$type = sqrt(drL$type^2 + drM$type^2) + set R$type = sqrt(dRL$type^2 + dRM$type^2) + + set phi$type = 1000*asin(dr$type / R$type) +# set flipSign = ((drL$type < 0) && (otaIDx < 40)) || ((drL$type > 0) && (otaIDx >=40)) + set flipSign = ((drL00 < 0) && (otaIDx < 40)) || ((drL00 > 0) && (otaIDx >=40)) + set phi$type = flipSign ? phi$type : 0 - phi$type + end + + lim -n 0 otaID -10 10; clear; box; plot otaID phi00 -pt ocir; plot -c indigo otaID phi11 -pt circ; plot -c blue otaID phi10 -pt plus; plot -c red otaID phi01 -pt x + lim -n 1 -10 10 -10 10; clear; box; plot rotrad phi00 -pt ocir; plot -c indigo rotrad phi11 -pt circ; plot -c blue rotrad phi10 -pt plus; plot -c red rotrad phi01 -pt x + +# lim -n 0 Lo Mo; clear; box; plot -pt cir -c blue -sz 1 Lo Mo; needles Lo Mo dL dM -c red -scale 100 + +# confirm that we understand the structure in libpscoords.c +# data pscoords.gpc1.v3.dat +# read dLx 1 dMx 2 rotrad 3 + +# dLx, dMx are offsets for the chip centers +# lim -n 1 dL dLx; clear; box; plot dL dLx + +# lim -n 2 dL dM; clear; box; plot dL dM -pt ocir -sz 2.0 -c red; plot dLx dMx -pt plus -sz 2.0 -c blue +end + +# result: (dL, dM) for chips with Nastro > 0 +macro make.pscoords.test + + # Can I generate the FP corner map from pscoords? + + exec rm -f $TEMPDIR/tmp.ota.CC.dat + output $TEMPDIR/tmp.ota.CC.dat + for i 0 otaID[] + if (Nastro[$i] < 1) continue + echo otaID[$i] {$NXch/2} {$NYch/2} + end + output stdout + + exec pscoords -dooff in=ota out=fp < $TEMPDIR/tmp.ota.CC.dat > $TEMPDIR/tmp.fp.dat + + data $TEMPDIR/tmp.fp.dat + read Lx 1 Mx 2 + + exec pscoords -nooff in=ota out=fp < $TEMPDIR/tmp.ota.CC.dat > $TEMPDIR/tmp.fp.dat + + data $TEMPDIR/tmp.fp.dat + read Lo 1 Mo 2 + +# set dL = Lx - Lo +# set dM = Mx - Mo +# lim -n 0 Lo Mo; clear; box; plot -pt cir -c blue -sz 1 Lo Mo; needles Lo Mo dL dM -c red -scale 100 + +# confirm that we understand the structure in libpscoords.c +# data pscoords.gpc1.v3.dat +# read dLx 1 dMx 2 rotrad 3 + +# dLx, dMx are offsets for the chip centers +# lim -n 1 dL dLx; clear; box; plot dL dLx + +# lim -n 2 dL dM; clear; box; plot dL dM -pt ocir -sz 2.0 -c red; plot dLx dMx -pt plus -sz 2.0 -c blue +end + +######## smf WCS ops ########### + +macro make.chip.points + if ($0 != 3) + echo "USAGE: make.chip.points (smffile) (chipnum)" + break + end + + local smffile chipname + $smffile = $1 + $chipnum = $2 + + rd -x {1 + 3*$chipnum} chp -head $smffile + + vlist chip_x {$NXch/2} 0 $NXch $NXch 0 + vlist chip_y {$NYch/2} 0 0 $NYch $NYch + set fpa_x_nan = zero(chip_x) + NAN + set fpa_y_nan = zero(chip_y) + NAN + + keyword chp NASTRO NASTRO + if ($NASTRO == 0) + echo "chip has no astrometry solution, skipping" + set fpa_x = fpa_x_nan + set fpa_y = fpa_y_nan + return + end + + keyword chp CRVAL1 CRVAL1 + keyword chp CRVAL2 CRVAL2 + keyword chp CRPIX1 CRPIX1 + keyword chp CRPIX2 CRPIX2 + keyword chp CRVAL1 CRVAL1 + keyword chp CRVAL2 CRVAL2 + keyword chp CDELT1 CDELT1 + keyword chp CDELT2 CDELT2 + + break -auto off + keyword chp NPLYTERM NPLYTERM + $myStatus = $STATUS + break -auto off + if ($myStatus == 0) + $NPLYTERM = 1 + end + + keyword chp PC001001 PC001001 + keyword chp PC001002 PC001002 + keyword chp PC002001 PC002001 + keyword chp PC002002 PC002002 + if ($NPLYTERM > 1) + keyword chp PCA1X0Y2 PCA1X0Y2 + keyword chp PCA2X0Y2 PCA2X0Y2 + keyword chp PCA1X1Y1 PCA1X1Y1 + keyword chp PCA2X1Y1 PCA2X1Y1 + keyword chp PCA1X2Y0 PCA1X2Y0 + keyword chp PCA2X2Y0 PCA2X2Y0 + if ($NPLYTERM > 2) + keyword chp PCA1X0Y3 PCA1X0Y3 + keyword chp PCA2X0Y3 PCA2X0Y3 + keyword chp PCA1X1Y2 PCA1X1Y2 + keyword chp PCA2X1Y2 PCA2X1Y2 + keyword chp PCA1X2Y1 PCA1X2Y1 + keyword chp PCA2X2Y1 PCA2X2Y1 + keyword chp PCA1X3Y0 PCA1X3Y0 + keyword chp PCA2X3Y0 PCA2X3Y0 + end + end + # echo $CRVAL1 $CRVAL2 + + set tmp_x = $CDELT1*(chip_x - $CRPIX1) + set tmp_y = $CDELT2*(chip_y - $CRPIX2) + + set fpa_x = $CRVAL1 + $PC001001*tmp_x + $PC001002*tmp_y + set fpa_y = $CRVAL2 + $PC002001*tmp_x + $PC002002*tmp_y + + if ($NPLYTERM > 1) + set tmp_x2 = tmp_x^2 + set tmp_xy = tmp_x*tmp_y + set tmp_y2 = tmp_y^2 + + set fpa_x = fpa_x + $PCA1X2Y0*tmp_x2 + $PCA1X1Y1*tmp_xy + $PCA1X0Y2*tmp_y2 + set fpa_y = fpa_y + $PCA2X2Y0*tmp_x2 + $PCA2X1Y1*tmp_xy + $PCA2X0Y2*tmp_y2 + + if ($NPLYTERM > 2) + set tmp_x3 = tmp_x2*tmp_x + set tmp_y3 = tmp_y2*tmp_y + + set fpa_x = fpa_x + $PCA1X3Y0*tmp_x3 + $PCA1X2Y1*tmp_x2*tmp_y + $PCA1X1Y2*tmp_x*tmp_y2 + $PCA1X0Y3*tmp_y3 + set fpa_y = fpa_y + $PCA2X3Y0*tmp_x3 + $PCA2X2Y1*tmp_x2*tmp_y + $PCA2X1Y2*tmp_x*tmp_y2 + $PCA2X0Y3*tmp_y3 + end + end + + delete -q tmp_x tmp_y tmp_x2 tmp_y2 tmp_xy tmp_x3 tmp_y3 +end + + +############################### old code of interest ################# + +macro make.distortion.map + if ($0 != 3) + echo "USAGE: make.distortion.map (smffile) (outroot)" + break + end + + local smffile outroot + $smffile = $1 + $outroot = $2 + + # the FPA image ranges from -250,000 to +250,000 in each dimension + $NX = 1000 + $NY = 1000 + $dX = 500 + $dY = 500 + + rd phu -head $smffile + + keyword phu IMNAXIS1 IMNAXIS1 + keyword phu IMNAXIS2 IMNAXIS2 + keyword phu CRVAL1 CRVAL1 + keyword phu CRVAL2 CRVAL2 + keyword phu CRPIX1 CRPIX1 + keyword phu CRPIX2 CRPIX2 + keyword phu CDELT1 CDELT1 + keyword phu CDELT2 CDELT2 + keyword phu PC001001 PC001001 + keyword phu PC001002 PC001002 + keyword phu PC002001 PC002001 + keyword phu PC002002 PC002002 + keyword phu PCA1X0Y2 PCA1X0Y2 + keyword phu PCA2X0Y2 PCA2X0Y2 + keyword phu PCA1X0Y3 PCA1X0Y3 + keyword phu PCA2X0Y3 PCA2X0Y3 + keyword phu PCA1X1Y1 PCA1X1Y1 + keyword phu PCA2X1Y1 PCA2X1Y1 + keyword phu PCA1X1Y2 PCA1X1Y2 + keyword phu PCA2X1Y2 PCA2X1Y2 + keyword phu PCA1X2Y0 PCA1X2Y0 + keyword phu PCA2X2Y0 PCA2X2Y0 + keyword phu PCA1X2Y1 PCA1X2Y1 + keyword phu PCA2X2Y1 PCA2X2Y1 + keyword phu PCA1X3Y0 PCA1X3Y0 + keyword phu PCA2X3Y0 PCA2X3Y0 + keyword phu NPLYTERM NPLYTERM + + if ($NPLYTERM != 3) + echo "smf has NPLYTERM = $NPLYTERM, not 3. skipping" + break + end + + # make an image of a fixed size (say 1000 x 1000) + mcreate zm $NX $NY + + # generate an image of the FPA coordinates: + set fpa_xm = $dX*(xramp(zm) - $NX / 2) + set fpa_ym = $dY*(yramp(zm) - $NY / 2) + + set tmp_xm = $CDELT1*(fpa_xm - $CRPIX1) + set tmp_ym = $CDELT2*(fpa_ym - $CRPIX2) + + set tpa_xm = $PC001001*tmp_xm + $PC001002*tmp_ym + set tpa_ym = $PC002001*tmp_xm + $PC002002*tmp_ym + + set tmp_x2m = tmp_xm^2 + set tmp_xym = tmp_xm*tmp_ym + set tmp_y2m = tmp_ym^2 + + set tpa_xm = tpa_xm + $PCA1X2Y0*tmp_x2m + $PCA1X1Y1*tmp_xym + $PCA1X0Y2*tmp_y2m + set tpa_ym = tpa_ym + $PCA2X2Y0*tmp_x2m + $PCA2X1Y1*tmp_xym + $PCA2X0Y2*tmp_y2m + + set tmp_x3m = tmp_x2m*tmp_xm + set tmp_y3m = tmp_y2m*tmp_ym + + set tpa_xm = tpa_xm + $PCA1X3Y0*tmp_x3m + $PCA1X2Y1*tmp_x2m*tmp_ym + $PCA1X1Y2*tmp_xm*tmp_y2m + $PCA1X0Y3*tmp_y3m + set tpa_ym = tpa_ym + $PCA2X3Y0*tmp_x3m + $PCA2X2Y1*tmp_x2m*tmp_ym + $PCA2X1Y2*tmp_xm*tmp_y2m + $PCA2X0Y3*tmp_y3m + + set tpa_dx = tpa_xm - tmp_xm + set tpa_dy = tpa_ym - tmp_ym +end + +macro choose.smfs.gpc1.old + + local t0 t1 mydate + + ctimes -abs 2010/01/01,00 -var t0 + ctimes -abs now -var t1 + + $exp_last = none + + # these queries run too slowly + for time $t0 $t1 + ctimes -ref $time -var mydate + queueinit camRun + queueload camRun -x "camtool -processedexp -dbname gpc1 -label LAP.PV3.20140730.final -dateobs_begin $mydate -limit 2" + queuesize camRun -var Nlines + +# queueload camRun -x "camtool -processedexp -dbname gpc1 -dateobs_begin $mydateS -dateobs_end $mydateE -limit 2" +# echo $Nlines $mydateS $mydateE +# $mydateS = $mydateE + if ($Nlines == 0) continue + + queue2book camRun camBook + queuedelete camRun + + book getword camBook page.000 exp_name -var exp_name + book getword camBook page.000 path_base -var path_base + book delete camBook + + if ("$exp_name" == "$exp_last") + echo "skipping $mydate : $exp_name = $exp_last" + continue + end + $exp_last = $exp_name + + $nebpath = `neb-locate -p $path_base.smf` + + echo $nebpath $exp_name + end +end + +############################3 + +if ($SCRIPT) + + if ($argv:n != 1) + echo "USAGE: distortion.sh (set)" + exit 2 + end + + choose.smfs $argv:0 + if ($STATUS) exit 0 + exit 1 +end Index: trunk/tools/eam/pscoords/get.chip.IDs =================================================================== --- trunk/tools/eam/pscoords/get.chip.IDs (revision 41267) +++ trunk/tools/eam/pscoords/get.chip.IDs (revision 41267) @@ -0,0 +1,9 @@ +#!/bin/csh -f +if ($#argv != 1) then + echo "USAGE: get.chip.IDs (smffile)" + exit 2 +endif + +echo $1 | fields -n hdr NASTRO EXTNAME | awk '{print $2, substr($3, 3, 2)}' + +exit 0 Index: trunk/tools/eam/pscoords/libps2coords.c =================================================================== --- trunk/tools/eam/pscoords/libps2coords.c (revision 41267) +++ trunk/tools/eam/pscoords/libps2coords.c (revision 41267) @@ -0,0 +1,769 @@ +/* Routine to convert microns in the focal plane to OTA/cell index and posn */ +/* Syntax: pscoords [args] < infile > outfile + * args: + * in=[sky | tp | fp | ota | cell] + * out=[sky | tp | fp | ota | cell] + * + * and other arguments required for various transformations: see syntax() + */ +/* See PSDC-730-003-01 (OTA FITS) and PSDC-730-005-00 (GPC Arch) for details */ +/* 120605 JT, added real refraction calculation and updated offsets */ +/* 080930 Sidik added TEST_MAIN - must be defined if you want main() test program. */ +/* 080904 Rev 1.1, renamed from fpota.c */ +/* 071023 Rev 1.0 John Tonry */ +#include +#include +#include +#include +#include +#include "pscoords.h" + +static char +#if defined(__GNUC__) && defined(__STDC__) +__attribute__ ((__unused__)) +#endif +rcsid[] = "$Id: pscoords.c,v 1.2 2008/09/06 03:26:48 jt Exp jt $"; + +#define MRAD2 (0.001) /* Convert from mrad to rad, used for PSC_OFFROT_T */ + +/* Overall focal plane offset and rotation (um, rad) */ +PSC_OFFROT_T psc_psc_fpoff={0.,0.,0.}; + +/* Offsets and rotations of each OTA (um, millirad) */ +PSC_OFFROT_T psc_otaoff[PSC_NX*PSC_NY]={ +/* Offsets and rotations derived from ~1000 IPP smf fit parameters from 2018/04 - 2020/01. + OTAs 00, 43, 70, 77 were manually set to 0.0 due to poor fits -- 2020.02.20 EAM + */ +{ 0.0, 0.0, 0.00}, { -51.6, 11.0, -6.59}, /* OTA00 , OTA10 */ +{ -13.5, -22.1, -5.96}, { -13.2, -70.5, -5.70}, /* OTA20 , OTA30 */ +{ 0.0, 0.0, 0.00}, { -17.2, 39.8, -5.88}, /* OTA40 , OTA50 */ +{ 10.0, -126.8, -6.41}, { 0.0, 0.0, 0.00}, /* OTA60 , OTA70 */ +{ -4.0, 40.4, -6.12}, { -33.0, -33.8, -5.66}, /* OTA01 , OTA11 */ +{ -8.7, -2.5, -5.28}, { -51.1, -44.4, -5.06}, /* OTA21 , OTA31 */ +{ 21.3, 69.4, -4.99}, { 32.1, 93.4, -5.20}, /* OTA41 , OTA51 */ +{ 34.4, -34.3, -5.48}, { -16.9, -60.7, -5.95}, /* OTA61 , OTA71 */ +{ -80.3, 18.7, -5.69}, { -24.4, 28.8, -5.48}, /* OTA02 , OTA12 */ +{ 20.9, -13.5, -5.33}, { 60.7, -62.8, -5.25}, /* OTA22 , OTA32 */ +{ 0.0, 0.0, 0.00}, { 16.4, 95.1, 5.40}, /* OTA42 , OTA52 */ +{ 16.5, 73.4, -5.38}, { 53.9, 24.4, -5.48}, /* OTA62 , OTA72 */ +{ 43.6, -30.4, -5.60}, { 35.7, 34.9, -5.58}, /* OTA03 , OTA13 */ +{ 54.0, -15.9, -5.52}, { 1.4, -45.0, -5.48}, /* OTA23 , OTA33 */ +{ 0.0, 0.0, 0.00}, { 13.0, 52.5, -5.49}, /* OTA43 , OTA53 */ +{ -25.4, 7.5, -5.51}, { -33.1, -42.0, -5.50}, /* OTA63 , OTA73 */ +{ 32.5, 25.5, -5.56}, { -33.7, -61.8, -5.60}, /* OTA04 , OTA14 */ +{ 38.8, -43.8, -5.55}, { 35.8, -65.6, -5.51}, /* OTA24 , OTA34 */ +{ -29.7, 66.3, -5.60}, { -38.0, -2.2, -5.53}, /* OTA44 , OTA54 */ +{ 4.9, 9.3, -5.55}, { -71.2, 7.7, -5.51}, /* OTA64 , OTA74 */ +{ 0.0, 0.0, 0.00}, { 30.7, 22.2, -5.52}, /* OTA05 , OTA15 */ +{ 0.0, 0.0, 0.00}, { 0.0, 0.0, 0.00}, /* OTA25 , OTA35 */ +{ 23.3, 119.2, -5.23}, { 46.7, -43.6, -5.36}, /* OTA45 , OTA55 */ +{ -41.1, -13.4, -5.45}, { -10.6, -4.0, 5.96}, /* OTA65 , OTA75 */ +{ -34.2, 44.3, -6.09}, { -39.1, 8.7, -5.54}, /* OTA06 , OTA16 */ +{ 11.2, -47.9, -5.21}, { -51.9, -101.9, -5.05}, /* OTA26 , OTA36 */ +{ -14.4, 10.9, -5.07}, { -29.5, 15.5, -5.26}, /* OTA46 , OTA56 */ +{ 21.1, -68.6, -5.60}, { -39.3, -19.1, -6.07}, /* OTA66 , OTA76 */ +{ 0.0, 0.0, 0.00}, { -13.5, 33.6, -6.47}, /* OTA07 , OTA17 */ +{ 2.0, -110.9, -5.98}, { -31.0, -73.5, -5.63}, /* OTA27 , OTA37 */ +{ 32.3, -7.6, -5.67}, { 0.0, 0.0, 0.00}, /* OTA47 , OTA57 */ +{ 3.3, 44.1, -6.57}, { 0.0, 0.0, 0.00} /* OTA67 , OTA77 */ +}; + +#ifdef PS1_OFFSETS +PSC_OFFROT_T psc_otaoff[PSC_NX*PSC_NY]={ /* Offsets and rotations derived for PS1 / GPC1 (see libpscoords.h) */ + { 0.0, 0.0, 0.00}, { -12.8, -52.7, -1.37}, /* OTA00 , OTA10 */ + { 6.4,-142.0, 0.04}, { -20.8, -48.7, -1.30}, /* OTA20 , OTA30 */ + { 22.0, 64.2, -0.59}, { 11.9, 36.4, -0.97}, /* OTA40 , OTA50 */ + { 34.1, -51.8, -1.77}, { 0.0, 0.0, 0.00}, /* OTA60 , OTA70 */ + { -20.1, -30.6, 0.64}, { -58.3, -71.8, 1.06}, /* OTA01 , OTA11 */ + { 17.6, -17.7, -0.69}, { -31.4, 11.1, -0.98}, /* OTA21 , OTA31 */ + { -11.0, 38.3, 0.88}, { -8.0, 43.9, -0.30}, /* OTA41 , OTA51 */ + { -20.9, -11.0, -1.58}, { 15.9, 4.2, 0.51}, /* OTA61 , OTA71 */ + { -27.0, 32.6, 2.63}, { -0.4, -60.4, 0.51}, /* OTA02 , OTA12 */ + { -3.6, -45.7, -1.61}, { -9.0,-143.6, 1.62}, /* OTA22 , OTA32 */ + { -20.7, 15.0, -1.20}, { -13.6, -4.6, 0.03}, /* OTA42 , OTA52 */ + { 6.9, -17.8, -2.02}, { 91.8, 31.8, 1.70}, /* OTA62 , OTA72 */ + { -15.3, 22.5, 0.14}, {-111.7, -78.6, 1.89}, /* OTA03 , OTA13 */ + { 65.8, 28.9, 1.09}, { -30.0, 12.3, 0.57}, /* OTA23 , OTA33 */ + { -31.8, 43.3, -0.69}, { 12.9, 19.7, 0.41}, /* OTA43 , OTA53 */ + { 17.8, -38.0, 1.39}, { 23.0, -73.5, -0.84}, /* OTA63 , OTA73 */ + { -38.7, 14.7, 0.13}, { -18.6, 2.5, 0.16}, /* OTA04 , OTA14 */ + { 31.3, 215.1, -0.35}, { -83.3, -99.4, 0.90}, /* OTA24 , OTA34 */ + { 5.8, -13.9, -0.16}, { -18.9, -11.7, -0.15}, /* OTA44 , OTA54 */ + { 11.6, 45.0, 1.22}, { 82.9, -7.2, 0.35}, /* OTA64 , OTA74 */ + { -24.8, 58.9, 2.63}, { 24.2, 15.6, -0.65}, /* OTA05 , OTA15 */ + { -51.2, 7.5, 1.67}, { -33.9, -38.8, -0.11}, /* OTA25 , OTA35 */ + { -44.7, 13.3, -0.76}, { 32.9, 72.2, 0.56}, /* OTA45 , OTA55 */ + { -42.0, 58.6, -0.44}, { 79.8, -32.2, 0.63}, /* OTA65 , OTA75 */ + { 7.5, 75.4, 0.07}, { -5.0, 39.4, -0.91}, /* OTA06 , OTA16 */ + { -77.8, -2.6, 1.29}, { 7.7, -12.9, 0.47}, /* OTA26 , OTA36 */ + { -4.7, 33.2, -1.39}, { 28.4, 41.2, 0.20}, /* OTA46 , OTA56 */ + { 20.0, 10.3, -1.60}, { 48.3, 47.6, -0.36}, /* OTA66 , OTA76 */ + { 0.0, 0.0, 0.00}, { 39.5, 37.5, 0.10}, /* OTA07 , OTA17 */ + { 78.8, 28.1, -0.71}, { -25.5, -3.3, 0.62}, /* OTA27 , OTA37 */ + { -57.3, 11.3, -0.98}, { -31.7, -0.5, -2.00}, /* OTA47 , OTA57 */ + { -22.3, 52.6, -0.34}, { 0.0, 0.0, 0.00} /* OTA67 , OTA77 */ +}; +#endif + +/* OTA offset reported by astrometry.sh between ROT0, ROTM, and ROT1 */ +PSC_ROTWEAK_T psc_otarot[PSC_NX*PSC_NY]={ /* Rotator-induced offsets not measured for PS2 / GPC2 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA00 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA10 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA20 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA30 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA40 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA50 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA60 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA70 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA01 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA11 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA21 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA31 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA41 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA51 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA61 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA71 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA02 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA12 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA22 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA32 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA42 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA52 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA62 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA72 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA03 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA13 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA23 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA33 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA43 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA53 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA63 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA73 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA04 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA14 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA24 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA34 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA44 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA54 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA64 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA74 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA05 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA15 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA25 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA35 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA45 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA55 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA65 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA75 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA06 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA16 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA26 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA36 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA46 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA56 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA66 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA76 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA07 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA17 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA27 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA37 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA47 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA57 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, /* OTA67 */ + { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}}; /* OTA77 */ + +/* Do offsets or no? + * & 0x01 => apply static OTA offsets + * & 0x02 => apply rotator offsets +*/ +static int DOOFFSET=1; + +#if 0 +#define IRAS /* Gnomonic? irsa.ipac.edu/IRASdocs/surveys/coordproj.html */ +#endif + +/* Refraction at observatory at wavelength, default is nominal Haleakala */ +/* Units are arcsec per tanz */ +static double psc_refract_coeff=PSC_REFRACT_CONST*PS_airdens; + +/* Refraction: rad per tanz. Set psc_refract_coeff for calculations below. */ +/* J.C. Owens (1967) from www.mpe.mpg.de/ir/instruments/sharp2+/adc/adc.html */ +int psc_refract(double T /* K */, double RH /* % */, + double P /* mbar */, double wave /* um */) +{ + double Psat, Pdry, Pwet, d1, d2, x2; + double sec2rad=atan(1.0)/(45*3600); + Psat = -10474 + 116.43*T - 0.43284*T*T + 0.0005384*T*T*T; + Pwet = (RH/100) * Psat; + Pdry = P - Pwet; + d1 = Pdry/T * (1 + Pdry*(57.9e-8 - 9.324e-4/T + 0.25844/(T*T))); + d2 = Pwet/T * (1 + Pwet*(1 + 3.7e-4*Pwet) * + (-2.37321e-3 + 2.23366/T - 710.792/(T*T) + 77514.1/(T*T*T))); + x2 = 1 / (wave*wave); + psc_refract_coeff = 1e-8/sec2rad* + (d1*(2371.34 + 683939.7/(130-x2) + 4547.3/(38.99-x2)) + + d2*(6487.31 + 58.058*x2 - 0.7115*x2*x2 + 0.08851*x2*x2*x2)); + return(0); +} + +/* Return default values for current fit */ +int psc_defaults(double *pscale, double *d2, double *d3, double *airdens) +{ + *pscale = PS_scale; + *d2 = PS_d2; + *d3 = PS_d3; + *airdens = PS_airdens; + return(0); +} + +/* Load up a new offset table? */ +int psc_load_otaoff(const char *fname) +{ + int i; + char rbuf[1024]; + FILE *fp; + if( (fp=fopen(fname, "r")) == NULL) { + fprintf(stderr, "error: Error opening OTA offset file '%s'\n", fname); + return(-1); + } + for(i=0; i 0 ? +1 : -1; + +/* Guess at which OTA this is */ + *ota_xid = NINT(PSC_NX/2-0.5-x/PSC_HMECH); + *ota_yid = NINT((y-lhs*0.5*(PSC_VMOFF+PSC_VMECH-PSC_VDIE-2*mechgap))/PSC_VMECH+(PSC_NY/2-0.5)); +/* Stop at the edges of the focal plane */ + if (*ota_xid < 0) *ota_xid = 0; + if (*ota_xid >= PSC_NX) *ota_xid = PSC_NX-1; + if (*ota_yid < 0) *ota_yid = 0; + if (*ota_yid >= PSC_NY) *ota_yid = PSC_NY-1; + idx = *ota_xid+*ota_yid*PSC_NX; + +/* Center of nominal mechanical PSC_HMECH,PSC_VMECH envelope wrt FP center */ + xc = ((PSC_NX/2-0.5)-*ota_xid) * PSC_HMECH; + yc = (*ota_yid-(PSC_NY/2-0.5)) * PSC_VMECH + lhs*0.5*(PSC_VMOFF+PSC_VMECH-PSC_VDIE-2*mechgap); + +/* Center of OTA (8x8 set of cells) within mechanical envelope wrt FP center */ + xc -= lhs*(0.5*PSC_HMECH - mechgap - PSC_LBORDER - PSC_NX/2*PSC_HCELL - (PSC_NX/2-0.5)*PSC_VSTREET); + yc -= lhs*(0.5*PSC_VMECH - mechgap - PSC_TBORDER - PSC_NY/2*PSC_VCELL - (PSC_NY/2-0.5)*PSC_HSTREET); + +/* Position wrt nominal OTA center prior to offset and rotation correction */ + xp = x - xc; + yp = y - yc; + +/* Position wrt OTA center after correction for offset and rotation */ + x = xp; + y = yp; +/* Correct for static OTA offsets? */ + if(DOOFFSET & 0x01) { + x = (xp-psc_otaoff[idx].dx)*cos(psc_otaoff[idx].rot*MRAD2) + + (yp-psc_otaoff[idx].dy)*sin(psc_otaoff[idx].rot*MRAD2); + y = -(xp-psc_otaoff[idx].dx)*sin(psc_otaoff[idx].rot*MRAD2) + + (yp-psc_otaoff[idx].dy)*cos(psc_otaoff[idx].rot*MRAD2); + } + +/* Position in pixels relative to OTA origin */ + *xota = (lhs*x + (PSC_NX/2*PSC_HCELL + (PSC_NX/2-0.5)*PSC_VSTREET)) / PSC_PIXEL; + *yota = -(lhs*y - (PSC_NY/2*PSC_VCELL + (PSC_NY/2-0.5)*PSC_HSTREET)) / PSC_PIXEL; + + return(0); +} + +/* Convert a focal plane position to OTA otax,otay, Cell cellx,celly, Pixel */ +/* This version can apply rotator-induced offsets if DOOFFSET permits */ +int psc_fp_to_pixel_rot(double xfp, double yfp, /* FP position (um) */ + int *ota_xid, int *ota_yid, /* OTA ID (0:7) */ + double *xota, double *yota, /* OTA position (pix) */ + double rotrad) /* Rotator [rad] */ +{ + double mechgap=0.5*(PSC_HMECH-PSC_HDIE); /* Assume equal L/R/T chip gap */ + int idx; /* Index of this OTA in otaoff */ + int lhs; /* Left/right hand side = +1/-1 */ + double xc, yc, x, y, xp, yp; + double rot, deg2rad=atan(1.0)/45; + +/* Position in FP after removing offset and rotation */ + if(DOOFFSET & 0x01) { + x = (xfp-psc_fpoff.dx)*cos(psc_fpoff.rot) + (yfp-psc_fpoff.dy)*sin(psc_fpoff.rot); + y = -(xfp-psc_fpoff.dx)*sin(psc_fpoff.rot) + (yfp-psc_fpoff.dy)*cos(psc_fpoff.rot); + } else { + x = xfp; + y = yfp; + } + +/* Guess at which side we are on */ + lhs = x > 0 ? +1 : -1; + +/* Guess at which OTA this is */ + *ota_xid = NINT(PSC_NX/2-0.5-x/PSC_HMECH); + *ota_yid = NINT((y-lhs*0.5*(PSC_VMOFF+PSC_VMECH-PSC_VDIE-2*mechgap))/PSC_VMECH+(PSC_NY/2-0.5)); +/* Stop at the edges of the focal plane */ + if (*ota_xid < 0) *ota_xid = 0; + if (*ota_xid >= PSC_NX) *ota_xid = PSC_NX-1; + if (*ota_yid < 0) *ota_yid = 0; + if (*ota_yid >= PSC_NY) *ota_yid = PSC_NY-1; + idx = *ota_xid+*ota_yid*PSC_NX; + +/* Center of nominal mechanical PSC_HMECH,PSC_VMECH envelope wrt FP center */ + xc = ((PSC_NX/2-0.5)-*ota_xid) * PSC_HMECH; + yc = (*ota_yid-(PSC_NY/2-0.5)) * PSC_VMECH + lhs*0.5*(PSC_VMOFF+PSC_VMECH-PSC_VDIE-2*mechgap); + +/* Center of OTA (8x8 set of cells) within mechanical envelope wrt FP center */ + xc -= lhs*(0.5*PSC_HMECH - mechgap - PSC_LBORDER - PSC_NX/2*PSC_HCELL - (PSC_NX/2-0.5)*PSC_VSTREET); + yc -= lhs*(0.5*PSC_VMECH - mechgap - PSC_TBORDER - PSC_NY/2*PSC_VCELL - (PSC_NY/2-0.5)*PSC_HSTREET); + +/* Position wrt nominal OTA center prior to offset and rotation correction */ + xp = x - xc; + yp = y - yc; + +/* Position wrt OTA center after correction for offset and rotation */ + x = xp; + y = yp; +/* Correct for static OTA offsets? */ + if(DOOFFSET & 0x01) { + x = (xp-psc_otaoff[idx].dx)*cos(psc_otaoff[idx].rot*MRAD2) + + (yp-psc_otaoff[idx].dy)*sin(psc_otaoff[idx].rot*MRAD2); + y = -(xp-psc_otaoff[idx].dx)*sin(psc_otaoff[idx].rot*MRAD2) + + (yp-psc_otaoff[idx].dy)*cos(psc_otaoff[idx].rot*MRAD2); + } + +/* Correct for rotator induced OTA offsets? */ + if(DOOFFSET & 0x02) { + rot = rotrad / deg2rad; + if(rot < ROTM) { + x -= psc_otarot[idx].dxm + (rot-ROTM)/(ROT0-ROTM) * + (psc_otarot[idx].dx0 - psc_otarot[idx].dxm); + y -= psc_otarot[idx].dym + (rot-ROTM)/(ROT0-ROTM) * + (psc_otarot[idx].dy0 - psc_otarot[idx].dym); + } else { + x -= psc_otarot[idx].dxm + (rot-ROTM)/(ROT1-ROTM) * + (psc_otarot[idx].dx1 - psc_otarot[idx].dxm); + y -= psc_otarot[idx].dym + (rot-ROTM)/(ROT1-ROTM) * + (psc_otarot[idx].dy1 - psc_otarot[idx].dym); + } + } + +/* Position in pixels relative to OTA origin */ + *xota = (lhs*x + (PSC_NX/2*PSC_HCELL + (PSC_NX/2-0.5)*PSC_VSTREET)) / PSC_PIXEL; + *yota = -(lhs*y - (PSC_NY/2*PSC_VCELL + (PSC_NY/2-0.5)*PSC_HSTREET)) / PSC_PIXEL; + + return(0); +} + +/* Convert a pixel position in OTA ota_xid,ota_yid to FP */ +int psc_pixel_to_fp(int ota_xid, int ota_yid, /* OTA ID (0:7) */ + double xota, double yota, /* OTA position (pix) */ + double *xfp, double *yfp) /* FP position (um) */ +{ + int idx=ota_xid+ota_yid*PSC_NX; /* Index of this OTA in otaoff */ + double mechgap=0.5*(PSC_HMECH-PSC_HDIE); /* Assume equal L/R/T chip gap */ + int lhs=ota_xid < PSC_NX/2 ? +1 : -1; /* Left/right hand side = +1/-1 */ + double xc, yc, x, y, xp, yp; + +/* Pixel position in um relative to OTA center */ + x = lhs * ( xota*PSC_PIXEL - (PSC_NX/2*PSC_HCELL + (PSC_NX/2-0.5)*PSC_VSTREET)); + y = lhs * (-yota*PSC_PIXEL + (PSC_NY/2*PSC_VCELL + (PSC_NY/2-0.5)*PSC_HSTREET)); + +/* Position wrt nominal OTA center after its offset and rotation */ + xp = x; + yp = y; +/* Correct for static OTA offsets? */ + if(DOOFFSET & 0x01) { + xp = psc_otaoff[idx].dx + x*cos(psc_otaoff[idx].rot*MRAD2) - y*sin(psc_otaoff[idx].rot*MRAD2); + yp = psc_otaoff[idx].dy + x*sin(psc_otaoff[idx].rot*MRAD2) + y*cos(psc_otaoff[idx].rot*MRAD2); + } + +/* Center of nominal mechanical PSC_HMECH,PSC_VMECH envelope wrt FP center */ + xc = ((PSC_NX/2-0.5)-ota_xid) * PSC_HMECH; + yc = (ota_yid-(PSC_NY/2-0.5)) * PSC_VMECH + lhs*0.5*(PSC_VMOFF+PSC_VMECH-PSC_VDIE-2*mechgap); + +/* Center of OTA (8x8 set of cells) within mechanical envelope wrt FP center */ + xc -= lhs*(0.5*PSC_HMECH - mechgap - PSC_LBORDER - PSC_NX/2*PSC_HCELL - (PSC_NX/2-0.5)*PSC_VSTREET); + yc -= lhs*(0.5*PSC_VMECH - mechgap - PSC_TBORDER - PSC_NY/2*PSC_VCELL - (PSC_NY/2-0.5)*PSC_HSTREET); + +/* Position within FP */ + x = xp + xc; + y = yp + yc; + +/* Final position from FP offset and rotation */ + *xfp = x; + *yfp = y; +/* Correct for static FP offset? */ + if(DOOFFSET & 0x01) { + *xfp = psc_fpoff.dx + x*cos(psc_fpoff.rot) - y*sin(psc_fpoff.rot); + *yfp = psc_fpoff.dy + x*sin(psc_fpoff.rot) + y*cos(psc_fpoff.rot); + } + + return(0); +} + +/* Convert a pixel position in OTA ota_xid,ota_yid to FP */ +int psc_pixel_to_fp_rot(int ota_xid, int ota_yid, /* OTA ID (0:7) */ + double xota, double yota, /* OTA position (pix) */ + double *xfp, double *yfp, /* FP position (um) */ + double rotrad) /* Rotator [rad] */ +{ + int idx=ota_xid+ota_yid*PSC_NX; /* Index of this OTA in otaoff */ + double mechgap=0.5*(PSC_HMECH-PSC_HDIE); /* Assume equal L/R/T chip gap */ + int lhs=ota_xid < PSC_NX/2 ? +1 : -1; /* Left/right hand side = +1/-1 */ + double xc, yc, x, y, xp, yp; + double rot, deg2rad=atan(1.0)/45; + +/* Pixel position in um relative to OTA center */ + x = lhs * ( xota*PSC_PIXEL - (PSC_NX/2*PSC_HCELL + (PSC_NX/2-0.5)*PSC_VSTREET)); + y = lhs * (-yota*PSC_PIXEL + (PSC_NY/2*PSC_VCELL + (PSC_NY/2-0.5)*PSC_HSTREET)); + +/* Position wrt nominal OTA center after its offset and rotation */ + xp = x; + yp = y; +/* Correct for static OTA offsets? */ + if(DOOFFSET & 0x01) { + xp = psc_otaoff[idx].dx + x*cos(psc_otaoff[idx].rot*MRAD2) - y*sin(psc_otaoff[idx].rot*MRAD2); + yp = psc_otaoff[idx].dy + x*sin(psc_otaoff[idx].rot*MRAD2) + y*cos(psc_otaoff[idx].rot*MRAD2); + } + +/* Correct for rotator induced OTA offsets? */ + if(DOOFFSET & 0x02) { + rot = rotrad / deg2rad; + if(rot < ROTM) { + xp += psc_otarot[idx].dxm + (rot-ROTM)/(ROT0-ROTM) * + (psc_otarot[idx].dx0 - psc_otarot[idx].dxm); + yp += psc_otarot[idx].dym + (rot-ROTM)/(ROT0-ROTM) * + (psc_otarot[idx].dy0 - psc_otarot[idx].dym); + } else { + xp += psc_otarot[idx].dxm + (rot-ROTM)/(ROT1-ROTM) * + (psc_otarot[idx].dx1 - psc_otarot[idx].dxm); + yp += psc_otarot[idx].dym + (rot-ROTM)/(ROT1-ROTM) * + (psc_otarot[idx].dy1 - psc_otarot[idx].dym); + } + } + +/* Center of nominal mechanical PSC_HMECH,PSC_VMECH envelope wrt FP center */ + xc = ((PSC_NX/2-0.5)-ota_xid) * PSC_HMECH; + yc = (ota_yid-(PSC_NY/2-0.5)) * PSC_VMECH + lhs*0.5*(PSC_VMOFF+PSC_VMECH-PSC_VDIE-2*mechgap); + +/* Center of OTA (8x8 set of cells) within mechanical envelope wrt FP center */ + xc -= lhs*(0.5*PSC_HMECH - mechgap - PSC_LBORDER - PSC_NX/2*PSC_HCELL - (PSC_NX/2-0.5)*PSC_VSTREET); + yc -= lhs*(0.5*PSC_VMECH - mechgap - PSC_TBORDER - PSC_NY/2*PSC_VCELL - (PSC_NY/2-0.5)*PSC_HSTREET); + +/* Position within FP */ + x = xp + xc; + y = yp + yc; + +/* Final position from FP offset and rotation */ + *xfp = x; + *yfp = y; +/* Correct for static FP offset? */ + if(DOOFFSET & 0x01) { + *xfp = psc_fpoff.dx + x*cos(psc_fpoff.rot) - y*sin(psc_fpoff.rot); + *yfp = psc_fpoff.dy + x*sin(psc_fpoff.rot) + y*cos(psc_fpoff.rot); + } + + return(0); +} + +/* Convert an OTA pixel position to cell ID and cell coords */ +int psc_pixel_to_cell(double xota, double yota, /* OTA position (pix) */ + int *cell_xid, int *cell_yid, /* Cell ID (0:7) or -1 */ + double *xcell, double *ycell) /* Cell position (pix) */ +{ + double x, y; + x = xota*PSC_PIXEL / (PSC_HCELL+PSC_VSTREET); + *cell_xid = (int)x; + *xcell = (x - (int)x) * (PSC_HCELL+PSC_VSTREET) / PSC_PIXEL; + if(x < 0 || *cell_xid > PSC_NX-1 || + (x-*cell_xid) > 1-PSC_VSTREET/(PSC_HCELL+PSC_VSTREET)) *cell_xid = -1; + *xcell = PSC_HCELL/PSC_PIXEL - *xcell; + + y = yota*PSC_PIXEL / (PSC_VCELL+PSC_HSTREET); + *cell_yid = (int)y; + *ycell = (y - (int)y) * (PSC_VCELL+PSC_HSTREET) / PSC_PIXEL; + if(y < 0 || *cell_yid > PSC_NY-1 || + (y-*cell_yid) > 1-PSC_HSTREET/(PSC_VCELL+PSC_HSTREET)) *cell_yid = -1; + + return(0); +} + +/* Convert a cell position and ID to OTA coords */ +int psc_cell_to_pixel(int cell_xid, int cell_yid, /* Cell ID (0:7) */ + double xcell, double ycell, /* Cell position (pix) */ + double *xota, double *yota) /* OTA position (pix) */ +{ + *xota = (PSC_HCELL/PSC_PIXEL - xcell) + cell_xid*(PSC_HCELL+PSC_VSTREET) / PSC_PIXEL; + *yota = ycell + cell_yid*(PSC_VCELL+PSC_HSTREET) / PSC_PIXEL; + return(0); +} + +/* Evaluate an optical model: arcsec -> microns notionally */ +/* Offset, scale with quadratic and cubic distortion, rotate */ +int psc_psoptics2(double x, double y, double dx, double dy, double dpa, + double pscale, double d2, double d3, double *xfp, double *yfp) +{ + double w2=((x-dx)*(x-dx)+(y-dy)*(y-dy)), w; + w = sqrt(w2); + x = pscale * (1.0 + d2*w + d3*w2) * (x-dx); + y = pscale * (1.0 + d2*w + d3*w2) * (y-dy); +/* Apply any extra TP-FP rotation */ + *xfp = x * cos(dpa) + y * sin(dpa); + *yfp = -x * sin(dpa) + y * cos(dpa); + return(0); +} + +/* Evaluate an inverse optical model: microns -> arcsec notionally */ +/* Offset, scale with *small* quadratic and cubic distortion, rotate */ +int psc_invoptics2(double xfp, double yfp, double dx, double dy, double dpa, + double pscale, double d2, double d3, double *x, double *y) +{ + double xp, yp, w; +/* Undo any extra TP-FP rotation */ + xp = (xfp * cos(dpa) - yfp * sin(dpa)) / pscale; + yp = (xfp * sin(dpa) + yfp * cos(dpa)) / pscale; + w = sqrt(xp*xp+yp*yp); + *x = xp / (1.0 + d2*w + d3*(xp*xp+yp*yp)) + dx; + *y = yp / (1.0 + d2*w + d3*(xp*xp+yp*yp)) + dy; + return(0); +} + + +/* Evaluate an optical model: arcsec -> microns notionally */ +/* Offset, scale with distortion, rotate */ +int psc_psoptics(double x, double y, double dx, double dy, double dpa, + double pscale, double d3, double *xfp, double *yfp) +{ + double w2=((x-dx)*(x-dx)+(y-dy)*(y-dy)); + x = pscale * (1.0 + d3*w2) * (x-dx); + y = pscale * (1.0 + d3*w2) * (y-dy); +/* Apply any extra TP-FP rotation */ + *xfp = x * cos(dpa) + y * sin(dpa); + *yfp = -x * sin(dpa) + y * cos(dpa); + return(0); +} + +/* Evaluate an inverse optical model: microns -> arcsec notionally */ +/* Offset, scale with *small* distortion, rotate */ +int psc_invoptics(double xfp, double yfp, double dx, double dy, double dpa, + double pscale, double d3, double *x, double *y) +{ + double xp, yp; +/* Undo any extra TP-FP rotation */ + xp = (xfp * cos(dpa) - yfp * sin(dpa)) / pscale; + yp = (xfp * sin(dpa) + yfp * cos(dpa)) / pscale; + *x = xp / (1.0 + d3*(xp*xp+yp*yp)) + dx; + *y = yp / (1.0 + d3*(xp*xp+yp*yp)) + dy; + return(0); +} + +/* Project RA,Dec (a,d) to the tangent plane (x,y) at (a0,d0,pa) [radians] */ +/* y is rotated CCW by pa from N; NOTE: x is west when PA = 0 (pos parity) */ +int psc_tproject(double a, double d, double a0, double d0, double pa, + double density, double zenith, double vertical, + double *x, double *y) +{ +#ifndef IRAS + double pdotx, pdoty, pdotz, E, N, dx, dy, z, Up, Horiz, sec2rad=atan(1.0)/(45*3600); + pdotx = cos(a0)*cos(d0)*cos(a)*cos(d) + + sin(a0)*cos(d0)*sin(a)*cos(d) + sin(d0)*sin(d); + pdoty = -sin(a0)*cos(a)*cos(d) + cos(a0)*sin(a)*cos(d); + pdotz = -cos(a0)*sin(d0)*cos(a)*cos(d) - + sin(a0)*sin(d0)*sin(a)*cos(d) + cos(d0)*sin(d); + E = pdoty / pdotx; + N = pdotz / pdotx; + *x = -E*cos(pa) + N*sin(pa); + *y = E*sin(pa) + N*cos(pa); +/* Differential refraction? */ + if(density > 0) { /* NB: set psc_refract_coeff from psc_refract! */ +/* Rotate to vertical */ + Horiz = *x * cos(vertical) + *y * sin(vertical); + Up = -*x * sin(vertical) + *y * cos(vertical); + z = sqrt(Horiz*Horiz+(zenith-Up)*(zenith-Up)); + dx = -sec2rad*psc_refract_coeff * tan(z)*Horiz/z; + dy = sec2rad*psc_refract_coeff * (tan(z)*(zenith-Up)/z - tan(zenith)); + Horiz += dx; + Up += dy; +/* Rotate back */ + *x = Horiz * cos(vertical) - Up * sin(vertical); + *y = Horiz * sin(vertical) + Up * cos(vertical); + } + +#else + double C, F, E, N; + C = cos(d)*cos(a-a0); + F = 1/(sin(d0)*sin(d)+C*cos(d0)); + N = -F*(cos(d0)*sin(d) - C*sin(d0)); + E = -F*cos(d)*sin(a-a0); + *x = E*cos(pa) + N*sin(pa); + *y = -E*sin(pa) + N*cos(pa); +#endif + return(0); +} + +/* Convert tangent plane (x,y) at (a0,d0,pa) to RA,Dec (a,d) [radians] */ +/* y is rotated CCW by pa from N; NOTE: x is west when PA = 0 (pos parity) */ +int psc_tplonglat(double x, double y, double a0, double d0, double pa, + double density, double zenith, double vertical, + double *a, double *d) +{ +#ifndef IRAS + double px, py, pz, E, N, pi=4*atan(1.0), dx, dy, z, sec2rad=atan(1.0)/(45*3600); + double Horiz, Up; +/* Undo differential refraction? */ + if(density > 0) { /* NB: set psc_refract_coeff from psc_refract! */ +/* Rotate to vertical */ + Horiz = x * cos(vertical) + y * sin(vertical); + Up = -x * sin(vertical) + y * cos(vertical); + z = sqrt(Horiz*Horiz+(zenith-Up)*(zenith-Up)); + dx = -sec2rad*psc_refract_coeff * tan(z)*Horiz/z; + dy = sec2rad*psc_refract_coeff * (tan(z)*(zenith-Up)/z - tan(zenith)); + Horiz -= dx; + Up -= dy; +/* Rotate back */ + x = Horiz * cos(vertical) - Up * sin(vertical); + y = Horiz * sin(vertical) + Up * cos(vertical); + } + E = -x*cos(pa) + y*sin(pa); + N = x*sin(pa) + y*cos(pa); + px = cos(a0)*cos(d0) - E*sin(a0) - N*cos(a0)*sin(d0); + py = sin(a0)*cos(d0) + E*cos(a0) - N*sin(a0)*sin(d0); + pz = sin(d0) + N*cos(d0); + *a = atan2(py, px); + if(pz >= 1) *d = pi/2; + else if(pz <= -1) *d = -pi/2; + else *d = asin(pz); +#else + double DD, B, XX, YY, E, N; + E = -x*cos(pa) + y*sin(pa); + N = x*sin(pa) + y*cos(pa); + DD = atan(sqrt(x*x+y*y)); + B = atan2(-x, y); + XX = sin(d0)*sin(DD)*cos(B) + cos(d0)*cos(DD); + YY = sin(DD)*sin(B); + *a = a0 + atan2(YY, XX); + *d = asin(sin(d0)*cos(DD)-cos(d0)*sin(DD)*cos(B)); +#endif + return(0); +} + +#ifdef TEST_MAIN +static +void bomb(char *msg) +{ + fprintf(stderr, "%s", msg); + exit(1); +} + +static +void syntax(char *prog) +{ + printf("Syntax: %s [options] < input_coords > output coords\n", prog); + printf("\nOptions:\n"); + printf(" in=[sky | tp | fp | ota | cell] Select input coord type\n"); + printf(" out=[sky | tp | fp | ota | cell] Select output coord type\n"); + printf("\nRequired extra arguments:\n"); + printf(" sky: requires ra=R dec=D pa=P for boresight [deg]\n"); + printf(" tp-fp: requires dx=X dy=Y [\"] dpa=P [deg] sc=S [um/\"] [d3=D [\"^-2]] for optical model [deg]\n"); + printf("\nOptional extra arguments:\n"); + printf(" offset=[t|f] Apply OTA/FP offsets (default f)\n"); + printf("\nInputs and outputs:\n"); + printf(" sky: RA, Dec [deg]\n"); + printf(" tp: x,y rotated by PA from N [arcsec]\n"); + printf(" fp: xfp,yfp [um] (aka camera coords)\n"); + printf(" ota: OTAxy ID and xota,yota [pix]\n"); + printf(" cell: OTAxy Cellxy IDs and xcell,ycell [pix]\n"); + printf(" (OTAxy can be omitted for OTA-Cell conversions only\n"); + printf("\nSynonyms:\n"); + printf(" ra=R or -ra R Boresight pointing (deg)\n"); + printf(" dec=D or -dec D Boresight pointing (deg)\n"); + printf(" pa=P or -pa P PA of image on sky (deg)\n"); + printf(" dx=X or -dx X FP offset wrt boresight (arcsec)\n"); + printf(" dy=Y or -dy Y FP offset wrt boresight (arcsec)\n"); + printf(" dpa=P or -dpa P FP rotation wrt PA (deg)\n"); + printf(" sc=S or -sc S Plate scale (um/arcsec)\n"); + printf(" d3=D or -d3 D Cubic distortiont (arcsec^-2)\n"); + printf(" offset=t or -dooff Apply OTA offsets? (def=f)\n"); + printf(" offset=f or -nooff Apply OTA offsets? (def=f)\n"); + printf(" in=sky out=tp or -skytp\n"); + printf(" in=tp out=fp or -tpfp\n"); + printf(" in=fp out=ota or -fpota\n"); + printf(" in=cell out=ota or -clota\n"); + printf(" in=ota out=fp or -otafp\n"); + printf("\nTo implement differential refraction in Sky-TP:\n"); + printf(" alt=A or -alt A Altitude (deg)\n"); + printf(" dens=D or -dens D P/T (wrt STP, 0.71 for HK)\n"); + printf(" vert=V or -vert V Rotator angle (deg)\n"); + printf("\nExamples:\n"); + printf(" echo 33 13 277.36 61.79 | pscoords out=sky in=cell ra=105 dec=7 pa=10 dx=0 dy=112 dpa=20 sc=38.856\n"); + printf(" -> 105.100000 7.100022\n"); + printf(" echo 105.1 7.1 | pscoords out=cell in=sky ra=105 dec=7 pa=10 dx=0 dy=112 dpa=20 sc=38.856\n"); + printf(" -> 33 13 277.36 61.79\n"); + printf("(note slight bug in Sky<->TP projections depending on distance from pole)\n"); + printf("\nVersion: %s\n", rcsid); + + exit(0); +} +#endif /* TEST_MAIN */ Index: trunk/tools/eam/pscoords/ps2coords.c =================================================================== --- trunk/tools/eam/pscoords/ps2coords.c (revision 41267) +++ trunk/tools/eam/pscoords/ps2coords.c (revision 41267) @@ -0,0 +1,626 @@ +/* Routine to convert microns in the focal plane to OTA/cell index and posn */ +/* Syntax: pscoords [args] < infile > outfile + * args: + * in=[sky | tp | fp | ota | cell] + * out=[sky | tp | fp | ota | cell] + * + * and other arguments required for various transformations: see syntax() + */ +/* See PSDC-730-003-01 (OTA FITS) and PSDC-730-005-00 (GPC Arch) for details */ +/* 120605 JT, added real refraction calculation and updated offsets */ +/* 080930 Sidik added TEST_MAIN - must be defined if you want main() test program. */ +/* 080904 Rev 1.1, renamed from fpota.c */ +/* 071023 Rev 1.0 John Tonry */ +#include +#include +#include +#include +#include +#include +#include "pscoords.h" + +void *bomb(char *); +void *syntax(char *); + +/* Enumeration of actions */ +#define ACTION_TEST -1 /* Test */ +#define ACTION_CONVERT 0 /* Generic conversion */ +#define ACTION_OTAFP 1 /* OTA -> FP coords */ +#define ACTION_FPOTA 2 /* FP -> OTA coords */ +#define ACTION_SKYTP 3 /* Sky -> tangent plane coords */ +#define ACTION_TPFP 4 /* tangent plane -> FP coords */ +#define ACTION_CLOTA 5 /* cell -> OTA coords */ + +/* Enumeration of coordinates. Note: order matters... */ +#define COORD_NONE 0 /* No coords */ +#define COORD_SKY 1 /* Sky coords: RA, Dec (rad) */ +#define COORD_TP 2 /* Tangent plane coords: (arcsec) */ +#define COORD_FP 3 /* Focal plane (camera) coords: (um) */ +#define COORD_OTA 4 /* OTA coords: (pix) */ +#define COORD_CELL 5 /* Cell coords: (pix) */ + +int main(int argc, char **argv) +{ + double xfp, yfp, xota, yota, xcell, ycell; + int i, ota_xid, ota_yid, cell_xid, cell_yid, verbose=0, action, doof; + int ota_id, cell_id; + int boresight, optical_model, refraction, src=COORD_NONE, dest=COORD_NONE; + double a0, d0, a, d, pa, x, y, alt; + double pi=4*atan(1.0), deg2rad=atan(1.0)/45, sec2rad=atan(1.0)/(45*3600); + double dx, dy, dpa, pscale, d2, d3, vertical, density; + char buf[1024], *src_name=NULL, *dest_name=NULL; + double wave, temp, humid; + char filter; + +/* Parse the arguments */ + if(argc < 2) syntax(argv[0]); + + action = ACTION_CONVERT; + d2 = 0.0; + d3 = 0.0; + doof = 1; /* Use static OTA offsets */ + boresight = optical_model = refraction = 0; + vertical = 0.0; /* PA of vertical on FP (PARANG-POSANG) */ + density = 0.0; /* Air density = P/(1bar) * (273K/T) */ + wave = -1; /* Wavelength for refraction */ + filter = 'r'; /* Filter for refraction */ + temp = PS_temp; /* Temperature for refraction */ + humid = PS_relhumid; /* Humidity for refraction */ + + alt = 90.0; /* Altitude of pointing */ + + for(i=1; i TPx TPy [arcsec]\n"); + } + + if(action == ACTION_TPFP && optical_model < 0x0f) { + bomb("pscoords -tpfp -dx dTPx(\") -dy dTPy(\") -dpa dPA(deg) -sc pscale(um/\") [-d3 D3(\"^-2)] [-d2 D2(\"^-1)]\n < TPx TPy [arcsec] => FPx FPy [um]\n"); + } + +/* Rotation based chip offsets require rotator value! */ + if( (doof&0x02) && !(refraction&0x01) ) { + bomb("pscoords needs rotator (as vert=R) for offset=3 (apply rotator-induced offsets)\n"); + } + + while(fgets(buf, 1024, stdin) != NULL) { + + switch(action) { + case ACTION_CONVERT: +/* Get the input data */ + switch(src) { + case COORD_SKY: + i = sscanf(buf, "%lf %lf", &a, &d); + if(i != 2) bomb("pscoords in=sky ... < RA Dec [deg]\n"); + a *= deg2rad; + d *= deg2rad; + break; + case COORD_TP: + i = sscanf(buf, "%lf %lf", &x, &y); + if(i != 2) bomb("pscoords in=tp ... < TPx TPy [arcsec]\n"); + break; + case COORD_FP: + i = sscanf(buf, "%lf %lf", &xfp, &yfp); + if(i != 2) bomb("pscoords in=fp ... < FPx FPy [um]\n"); + break; + case COORD_OTA: + ota_id = 0; + i = sscanf(buf, "%d %lf %lf", &ota_id, &xota, &yota); + if(i != 3) { + if(dest < COORD_OTA) bomb("pscoords in=ota ... < OTA_ij OTA_x y [pix]\n"); + i = sscanf(buf, "%lf %lf", &xota, &yota); + if(i != 2) bomb("pscoords in=ota out=cell ... < OTA_x y [pix]\n"); + } + ota_xid = ota_id / 10; + ota_yid = ota_id % 10; + break; + case COORD_CELL: + ota_id = 0; + i = sscanf(buf, "%d %d %lf %lf", + &ota_id, &cell_id, &xcell, &ycell); + if(i != 4) { + if(dest < COORD_OTA) bomb("pscoords in=cell ... < OTA_ij Cell_ij Cell_x y [pix]\n"); + i = sscanf(buf, "%d %lf %lf", &cell_id, &xcell, &ycell); + if(i != 3) bomb("pscoords in=cell out=ota ... < Cell_ij Cell_x y [pix]\n"); + } + ota_xid = ota_id / 10; + ota_yid = ota_id % 10; + cell_xid = cell_id / 10; + cell_yid = cell_id % 10; + break; + default: + ;; + break; + } +/* Migrate the coordinates */ + for(i=src; i!=dest; i+=(dest>src?+1:-1) ) { + switch(i) { + case COORD_SKY: + if(dest > i) { + psc_tproject(a, d, a0, d0, pa, density, pi/2-alt, vertical, &x, &y); + x /= sec2rad; + y /= sec2rad; + } + break; + case COORD_TP: + if(dest > i) { + if(d2 == 0.0) { + psc_psoptics(x, y, dx, dy, dpa, pscale, d3, &xfp, &yfp); + } else { + psc_psoptics2(x, y, dx, dy, dpa, pscale, d2, d3, &xfp, &yfp); + } + } else { + psc_tplonglat(x*sec2rad, y*sec2rad, a0, d0, pa, density, pi/2-alt, vertical, &a, &d); + } + break; + case COORD_FP: + if(dest > i) { +// psc_fp_to_pixel(xfp, yfp, &ota_xid, &ota_yid, &xota, &yota); + psc_fp_to_pixel_rot(xfp, yfp, &ota_xid, &ota_yid, &xota, &yota, vertical); + } else { + if(d2 == 0.0) { + psc_invoptics(xfp, yfp, dx, dy, dpa, pscale, d3, &x, &y); + } else { + psc_invoptics2(xfp, yfp, dx, dy, dpa, pscale, d2, d3, &x, &y); + } + } + break; + case COORD_OTA: + if(dest > i) psc_pixel_to_cell(xota, yota, &cell_xid, &cell_yid, &xcell, &ycell); +// else psc_pixel_to_fp(ota_xid, ota_yid, xota, yota, &xfp, &yfp); + else psc_pixel_to_fp_rot(ota_xid, ota_yid, xota, yota, &xfp, &yfp, vertical); + break; + case COORD_CELL: + if(dest > i) ; + else psc_cell_to_pixel(cell_xid,cell_yid, xcell,ycell, &xota,&yota); + break; + default: + bomb("Shouldn't get here...\n"); + break; + } + } +/* Tell us about the result */ + switch(dest) { + case COORD_SKY: + if(a < 0) a += 2*pi; + if(a > 2*pi) a -= 2*pi; + printf("%10.6f %10.6f\n", a/deg2rad, d/deg2rad); + break; + case COORD_TP: + printf("%9.3f %9.3f\n", x, y); + break; + case COORD_FP: + printf("%10.1f %10.1f\n", xfp, yfp); + break; + case COORD_OTA: + printf("%d%d %9.2f %9.2f\n", ota_xid, ota_yid, xota, yota); + break; + case COORD_CELL: + printf("%d%d %d%d %9.2f %9.2f\n", ota_xid, ota_yid, cell_xid, cell_yid, xcell, ycell); + break; + default: + bomb("Shouldn't get here either...\n"); + break; + } + + if(verbose > 0) { + printf("sky: %10.6f %10.6f\n", a/deg2rad, d/deg2rad); + printf("tp: %9.3f %9.3f\n", x, y); + printf("fp: %10.1f %10.1f\n", xfp, yfp); + printf("OTA: %d%d %9.2f %9.2f\n", ota_xid, ota_yid, xota, yota); + printf("cell: %d%d %d%d %9.2f %9.2f\n", ota_xid, ota_yid, cell_xid, cell_yid, xcell, ycell); + } + break; + + case ACTION_TEST: + case ACTION_OTAFP: + i = sscanf(buf, "%d %d %lf %lf", &ota_xid, &ota_yid, &xota, &yota); + if(i != 4) { + fprintf(stderr, "pscoords -otafp OTA_i j OTA_x y Cell_i j Cell_x y [pix] => FP_x y [um]\n"); + exit(1); + } + +// psc_pixel_to_fp(ota_xid, ota_yid, xota, yota, &xfp, &yfp); + psc_pixel_to_fp_rot(ota_xid, ota_yid, xota, yota, &xfp, &yfp, vertical); + psc_pixel_to_cell(xota, yota, &cell_xid, &cell_yid, &xcell, &ycell); + + if(action == 0) { + printf("OTA%d%d%7.1f %7.1f xy%d%d %7.1f %7.1f -> %9.1f %9.1f\n", + ota_xid,ota_yid, xota,yota, + cell_xid,cell_yid, xcell,ycell, xfp,yfp); + + +// psc_fp_to_pixel(xfp, yfp, &ota_xid, &ota_yid, &xota, &yota); + psc_fp_to_pixel_rot(xfp, yfp, &ota_xid, &ota_yid, &xota, &yota, vertical); + psc_pixel_to_cell(xota, yota, &cell_xid, &cell_yid, &xcell, &ycell); + + printf("%9.0f %9.0f -> OTA%d%d %7.1f %7.1f xy%d%d %7.1f %7.1f\n", + xfp,yfp, ota_xid,ota_yid, xota,yota, + cell_xid,cell_yid, xcell,ycell); + } else { + printf("%10.1f %10.1f\n", xfp, yfp); + } + + break; + + case ACTION_CLOTA: + i = sscanf(buf, "%d %d %lf %lf", &cell_xid, &cell_yid, &xcell, &ycell); + if(i != 4) { + fprintf(stderr, "pscoords -clota Cell_i j Cell_x y [pix] => OTA_x y [pix]\n"); + exit(1); + } + psc_cell_to_pixel(cell_xid, cell_yid, xcell, ycell, &xota, &yota); + printf("%9.2f %9.2f\n", xota, yota); + break; + + case ACTION_FPOTA: + i = sscanf(buf, "%lf %lf", &xfp, &yfp); + if(i != 2) { + fprintf(stderr, "pscoords -fpota FP_x y [um] => OTA_i j OTA_x y Cell_i j Cell_x y [pix]\n"); + exit(1); + } + + +// psc_fp_to_pixel(xfp, yfp, &ota_xid, &ota_yid, &xota, &yota); + psc_fp_to_pixel_rot(xfp, yfp, &ota_xid, &ota_yid, &xota, &yota, vertical); + psc_pixel_to_cell(xota, yota, &cell_xid, &cell_yid, &xcell, &ycell); + printf("OTA %d %d %9.2f %9.2f Cell %d %d %9.2f %9.2f\n", + ota_xid,ota_yid, xota,yota, cell_xid,cell_yid, xcell,ycell); + break; + + case ACTION_SKYTP: + i = sscanf(buf, "%lf %lf", &a, &d); + if(i != 2) { + fprintf(stderr, "pscoords -skytp -ra RA0 -dec Dec0 -pa PA\n"); + fprintf(stderr, " < RA Dec [deg] => TPx TPy [arcsec]\n"); + exit(1); + } + a *= deg2rad; + d *= deg2rad; + + psc_tproject(a, d, a0, d0, pa, density, pi/2-alt, vertical, &x, &y); + + printf("%9.3f %9.3f\n", x/deg2rad*3600, y/deg2rad*3600); + break; + + case ACTION_TPFP: + i = sscanf(buf, "%lf %lf", &x, &y); + if(i != 2) { + fprintf(stderr, "pscoords -tpfp -dx dTPx(\") -dy dTPy(\") -dpa dPA(deg) -sc pscale(um/\") [-d3 D3(\"^-2)] [-d2 D2(\"^-1)]\n"); + fprintf(stderr, " < TPx TPy [arcsec] => FPx FPy [um]\n"); + exit(1); + } + + if(d2 == 0.0) { + psc_psoptics(x, y, dx, dy, dpa, pscale, d3, &xfp, &yfp); + } else { + psc_psoptics2(x, y, dx, dy, dpa, pscale, d2, d3, &xfp, &yfp); + } + + printf("%10.1f %10.1f\n", xfp, yfp); + break; + + default: + printf("Unknown request\n"); + break; + } + } + + exit(0); +} + +void *bomb(char *msg) +{ + fprintf(stderr, "%s", msg); + exit(1); +} + +void *syntax(char *prog) +{ + printf("Syntax: %s [options] < input_coords > output coords\n", prog); + printf("\nOptions:\n"); + printf(" in=[sky | tp | fp | ota | cell] Select input coord type\n"); + printf(" out=[sky | tp | fp | ota | cell] Select output coord type\n"); + printf("\nRequired extra arguments:\n"); + printf(" sky: requires ra=R dec=D pa=P for boresight [deg]\n"); + printf(" tp-fp: requires dx=X dy=Y [\"] dpa=P [deg] sc=S [um/\"] [d3=D [\"^-2]] [d2=D [\"^-1]] for optical model [deg]\n"); + printf(" strongly recommended to use dx=X dy=Y [\"] dpa=P [deg] optics=ps\n"); + printf("\nOptional extra arguments:\n"); + printf(" offset=[t|f] Apply OTA/FP offsets (default f)\n"); + printf("\nInputs and outputs:\n"); + printf(" sky: RA, Dec [deg]\n"); + printf(" tp: x,y rotated by PA from N [arcsec]\n"); + printf(" fp: xfp,yfp [um] (aka camera coords)\n"); + printf(" ota: OTAxy ID and xota,yota [pix]\n"); + printf(" cell: OTAxy Cellxy IDs and xcell,ycell [pix]\n"); + printf(" (OTAxy can be omitted for OTA-Cell conversions only\n"); + printf("\nSynonyms:\n"); + printf(" ra=R or -ra R Boresight pointing (deg)\n"); + printf(" dec=D or -dec D Boresight pointing (deg)\n"); + printf(" pa=P or -pa P PA of image on sky (deg)\n"); + printf(" dx=X or -dx X FP offset wrt boresight (arcsec)\n"); + printf(" dy=Y or -dy Y FP offset wrt boresight (arcsec)\n"); + printf(" dpa=P or -dpa P FP rotation wrt PA (deg)\n"); + printf(" sc=S or -sc S Plate scale (um/arcsec)\n"); + printf(" d3=D or -d3 D Cubic distortion (arcsec^-2)\n"); + printf(" d2=D or -d2 D Quadratic distortion (arcsec^-1)\n"); + printf(" offset=t or -dooff Apply OTA offsets? (def=t)\n"); + printf(" offset=f or -nooff Apply OTA offsets? (def=t)\n"); + printf(" optics=ps or -psoptics Use defaults for scale, d2, d3, offsets\n"); + printf(" otaoff=fname Read dx,dy,rot [um,mrad] offset table\n"); + printf(" in=sky out=tp or -skytp\n"); + printf(" in=tp out=fp or -tpfp\n"); + printf(" in=fp out=ota or -fpota\n"); + printf(" in=cell out=ota or -clota\n"); + printf(" in=ota out=fp or -otafp\n"); + printf("\nTo implement differential refraction in Sky-TP:\n"); + printf(" alt=A or -alt A Altitude (deg)\n"); + printf(" vert=V or -vert V Rotator angle (deg)\n"); + printf(" dens=D or -dens D Pressure (bar, 0.71 for HK)\n"); + printf(" (optional additional refraction parameters)\n"); + printf(" filt=F or -filt F Filter (grizyO for eff. wave)\n"); + printf(" wave=W or -wave W Wavelength (nm)\n"); + printf(" temp=T or -temp T Temperature (K)\n"); + printf(" rh=R or -rh R Relative humidity (%%)\n"); + printf("\nExamples:\n"); + printf(" echo 44 02 332.21 197.35 | pscoords out=sky in=cell ra=105 dec=7 pa=20 dx=0 dy=112 dpa=20 optics=ps\n"); + printf(" -> 105.100000 7.100022\n"); + printf(" echo 105.1 7.1 | pscoords out=cell in=sky ra=105 dec=7 pa=20 dx=0 dy=112 dpa=20 optics=ps\n"); + printf(" -> 44 02 332.21 197.35\n"); +// printf("\nVersion: %s\n", rcsid); + + exit(0); +} Index: trunk/tools/eam/pscoords/ps2coords.h =================================================================== --- trunk/tools/eam/pscoords/ps2coords.h (revision 41267) +++ trunk/tools/eam/pscoords/ps2coords.h (revision 41267) @@ -0,0 +1,284 @@ +/* Header file for pscoords.c */ +#define NINT(x) (x<0?(int)((x)-0.5):(int)((x)+0.5)) + +#define OTA_CTE (5.0e-6) /* Coefficient of thermal expansion of Si+Mo */ +#define CFFP_CTE (0.5e-6) /* Coefficient of thermal expansion of CFFP */ +#define DELTA_T (-100.0) /* Operating temp minus metrology temp */ + +#define OTA_SCALE (1.0+OTA_CTE*DELTA_T) /* Room temp um to operating um of Si+Mo */ +#define CFFP_SCALE ((1.0+CFFP_CTE*DELTA_T)*25400) /* Room temp inches to um on CFFP */ + +/* Silicon layout of OTAs (um) */ +#define PSC_PIXEL (10.0*OTA_SCALE) /* Pixel size of an OTA */ +#define PSC_HCELL (5900.0*OTA_SCALE) /* Horizontal cell size */ +#define PSC_VCELL (5980.0*OTA_SCALE) /* Vertical cell size */ +#define PSC_HSTREET (120.0*OTA_SCALE) /* Horizontal street between cells (HS) */ +#define PSC_VSTREET (180.0*OTA_SCALE) /* Vertical street between cells (VS) */ +#define PSC_LBORDER (442.0*OTA_SCALE) /* Left border (to left-most cell) (LB) */ +#define PSC_RBORDER (418.0*OTA_SCALE) /* Right border (to right-most vertical street) (RB) */ +#define PSC_TBORDER (354.5*OTA_SCALE) /* Top border (to top-most cell) (TB) */ +#define PSC_BBORDER (975.5*OTA_SCALE) /* Bottom border (to lowest horizontal street) (BB) */ +#define PSC_HDIE (49500.0*OTA_SCALE) /* Horizontal size of silicon die */ +#define PSC_VDIE (50130.0*OTA_SCALE) /* Vertical size of silicon die */ + +/* Mechanical layout of OTAs (um) */ +#define PSC_HMECH (1.957*CFFP_SCALE) /* Horizontal mech spacing of OTA placement */ +#define PSC_VMECH (2.025*CFFP_SCALE) /* Vertical mech spacing of OTA placement */ +#define PSC_VMOFF (PSC_VDIE-(1.568+0.375)*CFFP_SCALE) /* Vert die offset between sides */ + +/* Number of OTAs populating focal plane */ +#define PSC_NX 8 /* Number of OTAs in FP in x */ +#define PSC_NY 8 /* Number of OTAs in FP in y */ + +/* Cell layout, looking down on back illuminated cell: + * NOTE: center of pixel 0,0 is physical location 0.5,0.5 pixel; + * cell origin is lower left edge of P0,0. + * + * +------------------+ + * | P0,597 P589,597 | + * | ... ... | + * Amplifier<<| P0,0 P589,0 | + * +------------------+ + */ + +/* OTA layout, looking down on back illuminated CCD: + * NOTE: OTA coordinates increase in x from Cell xy0n to Cell xy7n, y from + * Cell xyn0 to Cell xyny, origin is P0,0 in C00, units are pixels. + * "center" of OTA is the middle of the streets between cell 3 and 4 + * therefore center is at (xota, yota) = (2423.0, 2434.0) + * Cell coordinates are flipped wrt OTA coordinates, so the x=589 + * pixel is the left most one in OTA coordinates, i.e in OTA + * coordinates a cell looks like: +---------------+ + * | 589,597 0,597 | + * | ... ... | + * | 589,0 0,0 |> Amp + * +---------------+ + * +--------Top----------+ + * | TB | + * L C07 VS ... C77 VS R Cnn is 590x598 10um pixels + * e HS HS i + * fLB... ... RBg + * t C00 VS ... C70 VS h + * | HS HS t + * | BB | + * +-------Bottom--------+ (amp side) + */ + +/* FP layout, looking down on FP: + * + * NOTE: edge to edge spacing of dice is HMECH,VMECH; + * RHS vertically offset by VMOFF from LHS + * to align pixel i,597 on LHS with pixel i,0 on RHS; + * + * "center" is midpoint: OTA33 Cell00 Pix00 and OTA44 Cell00 Pix00 + * + * Physical coordinate convention runs right to left in x, in + * order to comply with the "sky view" convention. + * + * +--B--+ +--B--+ +--B--+ (VMOFF) + * | | | | | | +--T--+ +--T--+ +--T--+ + * R 17 L R 27 L R 37 L | | | | | | + * | | | | | | L 47 R L 57 R L 67 R + * +--T--+ +--T--+ +--T--+ | | | | | | + * +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ + * | | | | | | | | +--T--+ +--T--+ +--T--+ +--T--+ + * R 06 L R 16 L R 26 L R 36 L | | | | | | | | + * | | | | | | | | L 46 R L 56 R L 66 R L 76 R + * +--T--+ +--T--+ +--T--+ +--T--+ | | | | | | | | + * +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ + * | | | | | | | | +--T--+ +--T--+ +--T--+ +--T--+ + * R 05 L R 15 L R 25 L R 35 L | | | | | | | | + * | | | | | | | | L 45 R L 55 R L 65 R L 75 R + * +--T--+ +--T--+ +--T--+ +--T--+ | | | | | | | | + * +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ + * | | | | | | | | +--T--+ +--T--+ +--T--+ +--T--+ + * R 04 L R 14 L R 24 L R 34 L | | | | | | | | + * | | | | | | | | L 44 R L 54 R L 64 R L 74 R + * +--T--+ +--T--+ +--T--+ +--T--+ | | | | | | | | + * +--B--+ +--B--+ +--B--+ +--B--+X+--B--+ +--B--+ +--B--+ +--B--+ + * | | | | | | | | +--T--+ +--T--+ +--T--+ +--T--+ + * R 03 L R 13 L R 23 L R 33 L | | | | | | | | + * | | | | | | | | L 43 R L 53 R L 63 R L 73 R + * +--T--+ +--T--+ +--T--+ +--T--+ | | | | | | | | + * +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ + * | | | | | | | | +--T--+ +--T--+ +--T--+ +--T--+ + * R 02 L R 12 L R 22 L R 32 L | | | | | | | | + * | | | | | | | | L 42 R L 52 R L 62 R L 72 R + * +--T--+ +--T--+ +--T--+ +--T--+ | | | | | | | | + * +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ + * | | | | | | | | +--T--+ +--T--+ +--T--+ +--T--+ + * R 01 L R 11 L R 21 L R 31 L | | | | | | | | + * | | | | | | | | L 41 R L 51 R L 61 R L 71 R + * +--T--+ +--T--+ +--T--+ +--T--+ | | | | | | | | + * +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ +--B--+ + * | | | | | | +--T--+ +--T--+ +--T--+ + * R 10 L R 20 L R 30 L | | | | | | y + * | | | | | | L 40 R L 50 R L 60 R ^ + * +--T--+ +--T--+ +--T--+ | | | | | | | + * +--B--+ +--B--+ +--B--+ | + * x <----+ + */ + +/* Structure for tweaks of OTA positions wrt mechanical */ +typedef struct { + double dx; /* OTA x origin found at dx wrt mechanical (um) */ + double dy; /* OTA y origin found at dy wrt mechanical (um) */ + double rot; /* OTA rotated CCW by rot wrt mechanical (rad) */ +} PSC_OFFROT_T; + +/* Overall focal plane offset and rotation (um, rad) */ +PSC_OFFROT_T psc_fpoff; + +/* Offsets and rotations of each OTA (um, millirad) */ +PSC_OFFROT_T psc_otaoff[PSC_NX*PSC_NY]; + + +#define ROT0 (-200.0) +#define ROTM (-115.0) +#define ROT1 ( -50.0) + +/* Structure for tweaks of OTA positions wrt rotator 120610 */ +typedef struct { + double dx0; /* OTA x offset at ROT0 (um) */ + double dxm; /* OTA x offset at ROTM (um) */ + double dx1; /* OTA x offset at ROT1 (um) */ + double dy0; /* OTA y offset at ROT0 (um) */ + double dym; /* OTA y offset at ROTM (um) */ + double dy1; /* OTA y offset at ROT1 (um) */ +} PSC_ROTWEAK_T; + +/* Offsets and rotations of each OTA created by rotation */ +PSC_ROTWEAK_T psc_otarot[PSC_NX*PSC_NY]; + +/* This is from Run 3, 090516 */ +#define PS_scale 38.7932 /* Default PS plate scale [um/arcsec] */ +#define PS_d2 9.78e-7 /* Default quadratic PS distortion [arcsec^-1] */ +#define PS_d3 3.16e-11 /* Default cubic PS distortion [arcsec^-2] */ + + +#ifdef RUN_THREE_V1 /* This is from Run 3, 090228 */ +#define PS_scale 38.860 /* Default PS plate scale [um/arcsec] */ +#define PS_d3 1.49e-10 /* Default cubic PS distortion [arcsec^-2] */ +#endif + +#define PS_airdens 0.71 /* Default PS1 air density (Haleakala) */ +#define PS_relhumid 45.0 /* Default PS1 relative humidity (%) */ +#define PS_temp 284.0 /* Default PS1 temperature (K) */ +#define PSC_REFRACT_CONST 60.0 /* Standard refraction ("/tanz) at STP */ + +/**********************/ +/* Normal Application */ +/**********************/ +/* + * To go from a position (a,d) relative to (a0,d0,pa) pointing to cell coords: + * --------------------------------------------------------------------------- + * psc_tproject(a, d, a0, d0, pa, density, pi/2-alt, vertical, &x, &y); + * psc_psoptics2(x/sec2rad, y/sec2rad, dx, dy, dpa, pscale, d2, d3, &xfp, &yfp); + * psc_fp_to_pixel(xfp, yfp, &ota_xid, &ota_yid, &xota, &yota); + * psc_pixel_to_cell(xota, yota, &cell_xid, &cell_yid, &xcell, &ycell); + * + * can also use + * psc_psoptics(x/sec2rad, y/sec2rad, dx, dy, dpa, pscale, d3, &xfp, &yfp); + * + * + * To go from cell coords to a position (a,d) relative to (a0,d0,pa) pointing: + * --------------------------------------------------------------------------- + * psc_cell_to_pixel(cell_xid,cell_yid, xcell,ycell, &xota,&yota); + * psc_pixel_to_fp(ota_xid, ota_yid, xota, yota, &xfp, &yfp); + * psc_invoptics2(xfp, yfp, dx, dy, dpa, pscale, d2, d3, &x, &y); + * psc_tplonglat(x*sec2rad, y*sec2rad, a0, d0, pa, density, pi/2-alt, vertical, &a, &d); + * + * can also use + * psc_invoptics(xfp, yfp, dx, dy, dpa, pscale, d3, &x, &y); + */ + +/**************/ +/* Prototypes */ +/**************/ +/* Enable application of chip offsets? */ +int psc_do_chipoff(int doit); /* Apply chip offsets? (0/1, default=1) */ + +/* Convert a focal plane position to OTA otax,otay, Cell cellx,celly, Pixel */ +int psc_fp_to_pixel(double xfp_um, double yfp_um, /* FP position (um) */ + int *ota_xid, int *ota_yid, /* OTA ID (0:7) */ + double *xota_pix, double *yota_pix); /* OTA position (pix) */ + +/* Convert a focal plane position to OTA otax,otay, Cell cellx,celly, Pixel */ +int psc_fp_to_pixel_rot(double xfp_um, double yfp_um, /* FP position (um) */ + int *ota_xid, int *ota_yid, /* OTA ID (0:7) */ + double *xota_pix, double *yota_pix, /* OTA posn (pix) */ + double rotrat); /* rotator [rad] */ + +/* Convert a pixel position in OTA ota_xid,ota_yid to FP */ +int psc_pixel_to_fp(int ota_xid, int ota_yid, /* OTA ID (0:7) */ + double xota_pix, double yota_pix, /* OTA position (pix) */ + double *xfp_um, double *yfp_um); /* FP position (um) */ + +/* Convert a pixel position in OTA ota_xid,ota_yid to FP, rotator version */ +int psc_pixel_to_fp_rot(int ota_xid, int ota_yid, /* OTA ID (0:7) */ + double xota_pix, double yota_pix, /* OTA posn (pix) */ + double *xfp_um, double *yfp_um, /* FP position (um) */ + double rotrad); /* rotator [rad] */ + +/* Convert an OTA pixel position to cell ID and cell coords */ +int psc_pixel_to_cell(double xota_pix, double yota_pix, /* OTA position (pix) */ + int *cell_xid, int *cell_yid, /* Cell ID (0:7) or -1 */ + double *xcell_pix, double *ycell_pix);/* Cell position (pix) */ + +/* Convert a cell position and ID to OTA coords */ +int psc_cell_to_pixel(int cell_xid, int cell_yid, /* Cell ID (0:7) */ + double xcell_pix, double ycell_pix, /* Cell position (pix) */ + double *xota_pix, double *yota_pix); /* OTA position (pix) */ + +/* Convert tangent plane (x,y) at (a0,d0,pa) to RA,Dec (a,d) [radians] */ +/* y is rotated CCW by pa from N; NOTE: x is west when PA = 0 (pos parity) */ +/* Set density=0 to apply no correction for differential refraction, otherwise P/T in STP */ +int psc_tplonglat(double x_rad, double y_rad, double a0_rad, double d0_rad, double pa_rad, + double density/*0.71 for HK*/, double zenith_rad/*pi/2-ALT*/, + double vertical_rad/*ROT*/, double *a_rad, double *d_rad); + +/* Project RA,Dec (a,d) to the tangent plane (x,y) at (a0,d0,pa) [radians] */ +/* y is rotated CCW by pa from N; NOTE: x is west when PA = 0 (pos parity) */ +/* Set density=0 to apply no correction for differential refraction, otherwise P/T in STP */ +int psc_tproject(double a_rad, double d_rad, double a0_rad, double d0_rad, double pa_rad, + double density/*0.71 for HK*/, double zenith_rad/*pi/2-ALT*/, + double vertical_rad/*ROT*/, double *x_rad, double *y_rad); + +/* Evaluate an inverse optical model: microns -> arcsec notionally */ +/* Offset, scale with *small* distortion, rotate */ +int psc_invoptics(double xfp_um, double yfp_um, double dx_sec, double dy_sec, double dpa_rad, + double pscale/*38.86um/sec*/, double d3/*1.49e-10sec^-2*/, + double *x_sec, double *y_sec); + +/* Evaluate an optical model: arcsec -> microns notionally */ +/* Offset, scale with distortion, rotate */ +int psc_psoptics(double x_sec, double y_sec, double dx_sec, double dy_sec, double dpa_rad, + double pscale/*38.86um/sec*/, double d3/*1.49e-10sec^-2*/, + double *xfp_um, double *yfp_um); + +/* Evaluate an inverse optical model: microns -> arcsec notionally */ +/* Offset, scale with *small* distortion, rotate */ +int psc_invoptics2(double xfp_um, double yfp_um, double dx_sec, double dy_sec, double dpa_rad, + double pscale/*38.793um/sec*/, + double d2/*9.8e-7sec^-2*/, double d3/*3.1e-11sec^-2*/, + double *x_sec, double *y_sec); + +/* Evaluate an optical model: arcsec -> microns notionally */ +/* Offset, scale with distortion, rotate */ +int psc_psoptics2(double x_sec, double y_sec, double dx_sec, double dy_sec, double dpa_rad, + double pscale/*38.793um/sec*/, + double d2/*9.8e-7sec^-2*/, double d3/*3.1e-11sec^-2*/, + double *xfp_um, double *yfp_um); + +/* Apply external set of chip offsets (um, um, mrad; OTA00, OTA10, ... OTA77) */ +int psc_load_otaoff(const char *fname); + +/* Set refraction coefficient for a particular observation */ +int psc_refract(double T /* K */, double RH /* % */, + double P /* mbar */, double wave /* um */); + +/* Return default values for current fit */ +int psc_defaults(double *pscale, double *d2, double *d3, double *airdens); + +/* Load up a new rotator offset table? */ +int psc_load_otarot(const char *fname); Index: trunk/tools/eam/pscoords/pscoords.sh =================================================================== --- trunk/tools/eam/pscoords/pscoords.sh (revision 41267) +++ trunk/tools/eam/pscoords/pscoords.sh (revision 41267) @@ -0,0 +1,190 @@ +# I have a collection of distortion model files (psc.fits). I want to generate summary stats for these files: + +macro merge.psc.data + if ($0 != 2) + echo "USAGE: merge.psc.data (dir)" + break + end + + local outdir + $outdir = $1 + + # summary plots for ROT & SCALE + exec ls $outdir/*.psc.fits | fields ROT1 ROT2 SCALE1 SCALE2 ROTANGLE POSANGLE ALT AZ FILTERID CRVAL1 CRVAL2 DETTEM M2Z MJDOBS > $outdir/exposures.dat + data $outdir/exposures.dat + read ROT1 2 ROT2 3 SCALE1 4 SCALE2 5 ROTANGLE 6 POSANGLE 7 ALT 8 AZ 9 CRVAL1 11 CRVAL2 12 DETTEM 13 M2Z 14 MJDOBS 15 + + $BS = -tickpad 0.2 -xpad 3.5 -labelpadx 3 + + resize -n 0 2000 800; antialias 0.4 + lim MJDOBS 1.00 1.02; clear; label -fn helvetica 18; + box $BS; plot -c blue MJDOBS SCALE1; plot -c red MJDOBS SCALE2 + label -x MJD -y "SCALE1 (blue), SCALE2 (red)" + png -name $outdir/scale.mjd.png + + set dS = 100*(SCALE1 - SCALE2) + vstat -q dS + $MIN_PLOT = $MEDIAN - 0.05 + $MAX_PLOT = $MEDIAN + 0.05 + + resize -n 1 2000 800;; antialias 0.4 + lim POSANGLE $MIN_PLOT $MAX_PLOT; clear; label -fn helvetica 18; + box $BS; + plot -c blue POSANGLE dS where (MJDOBS < 58760.322741) + plot -c red POSANGLE dS where (MJDOBS > 58760.322741) + label -x POSANGLE -y "dS (MJD < 58760 : blue, MJD > 58760 : red)" + png -name $outdir/scale.pos.png + + resize -n 2 2000 800;; antialias 0.4 + lim AZ $MIN_PLOT $MAX_PLOT; clear; label -fn helvetica 18; + box $BS; + plot -c blue AZ dS where (MJDOBS < 58760.322741) + plot -c red AZ dS where (MJDOBS > 58760.322741) + label -x Azimuth -y "dS (MJD < 58760 : blue, MJD > 58760 : red)" + png -name $outdir/scale.az.png + + list name -glob $outdir/*.psc.fits + + data pscoords.corners.dat + read otaID:int 1 + + mcreate zim $name:n otaID[] + set dL_im = NAN + zero(zim) + set dM_im = NAN + zero(zim) + set phi_im = NAN + zero(zim) + set dphi_im = NAN + zero(zim) + + for i 0 $name:n + load.psc.data $name:$i + + mset dL_im dL -y $i + mset dM_im dM -y $i + mset phi_im phi -y $i + mset dphi_im dphi -y $i + end + + set dL_ave = NAN + zero(otaID) + set dL_std = NAN + zero(otaID) + + set dM_ave = NAN + zero(otaID) + set dM_std = NAN + zero(otaID) + + set Phi_ave = NAN + zero(otaID) + set Phi_std = NAN + zero(otaID) + + for i 0 otaID[] + mget dL_im dL_$i -x $i + mget dM_im dM_$i -x $i + mget phi_im phi_$i -x $i + mget dphi_im dphi_$i -x $i + + histogram dL_$i ndL_$i -400 400 1 -range dx + hist.stats dx ndL_$i + dL_ave[$i] = $MEDIAN + dL_std[$i] = $SIGMA + + histogram dM_$i ndM_$i -400 400 1 -range dx + hist.stats dx ndM_$i + dM_ave[$i] = $MEDIAN + dM_std[$i] = $SIGMA + + histogram phi_$i nphi_$i -20 20 0.01 -range dx + hist.stats dx nphi_$i + Phi_ave[$i] = $MEDIAN + Phi_std[$i] = $SIGMA + end + + data pscoords.gpc1.v3.dat + read dL_psc 1 dM_psc 2 Phi_psc 3 ota_psc 5 + sort ota_psc dL_psc dM_psc Phi_psc + + set ddL = dL_ave - dL_psc + set ddM = dM_ave - dM_psc + set dPhi = Phi_ave - Phi_psc + + resize -n 3 2100 800; clear; label -fn helvetica 18; antialias 0.4 + section a 0.00 0.0 0.66 1.0 + lim otaID -250 250; box $BS; plot -c blue -pt cir otaID ddL -dy dL_std; label -x OTA -y "dL (micron) - PSC V3" + section b 0.66 0.0 0.33 1.0 + lim -250 250 -250 250; box $BS; plot -c blue -pt cir dL_psc dL_ave -dy dL_std; label -x "dL PSC V3" -y "dL IPP" + png -name $outdir/dL.png + + resize -n 4 2100 800; clear; label -fn helvetica 18; antialias 0.4 + section a 0.00 0.0 0.66 1.0 + lim otaID -250 250; box $BS; plot -c blue -pt cir otaID ddM -dy dM_std; label -x OTA -y "dM (micron) - PSC V3" + section b 0.66 0.0 0.33 1.0 + lim -250 250 -250 250; box $BS; plot -c blue -pt cir dM_psc dM_ave -dy dM_std; label -x "dM PSC V3" -y "dM IPP" + png -name $outdir/dM.png + + resize -n 5 2100 800; clear; label -fn helvetica 18; antialias 0.4 + section a 0.00 0.0 0.66 1.0 + lim otaID -5 5; box $BS; plot -c blue -pt cir otaID dPhi -dy Phi_std; label -x OTA -y "Rot (millirad) - PSC V3" + section b 0.66 0.0 0.33 1.0 + lim -5 5 -5 5; box $BS; plot -c blue -pt cir Phi_psc Phi_ave -dy Phi_std; label -x "Rot PSC V3" -y "Rot IPP" + png -name $outdir/Phi.png + + data pscoords.corners.dat + read Lo 2 Mo 3 otaIDx:int 1 + join otaID otaIDx + reindex Lo_s = Lo using index2 + reindex Mo_s = Mo using index2 + reindex otaID_s = otaIDx using index2 + + dev -n LM; lim Lo_s Mo_s; clear; box; plot -pt ocir -c blue Lo_s Mo_s; needles Lo_s Mo_s ddL ddM -scale 100 + png -name $outdir/LM.dLdM.png + + write -header $outdir/pscoords.params.dat otaID dL_ave dL_std dM_ave dM_std Phi_ave Phi_std +end + +# NOTE: assumes I have read the list of OTAs from pscoords.corners.dat +macro load.psc.data + if ($0 != 2) + echo "USAGE: load.psc.data (psc.fits)" + break + end + + data $1 + read -fits DATA otaID_mat dLoR dMoR phi_ave phi_std + + join -outer otaID otaID_mat + reindex -keep-unmatched dL = dLoR using index1 + reindex -keep-unmatched dM = dMoR using index1 + reindex -keep-unmatched phi = phi_ave using index1 + reindex -keep-unmatched dphi = phi_std using index1 +end + +$ERR_ONESIGMA_POS = (1 + erf(1/sqrt(2.0)))/2.0 +$ERR_ONESIGMA_MIN = (1 - erf(1/sqrt(2.0)))/2.0 + +## UTIL: +macro hist.stats + if ($0 != 3) + echo "USAGE: stats (dx) (Nx)" + echo " generate histogram with reasonable sampling" + break + end + + $2[0] = 0 + $2[-1] = 0 + + cumulative $2 _Stmp + if (_Stmp[-1] < 5) + $SIGMA = NAN + $MEDIAN = NAN + $NPTS = _Stmp[-1] + return + end + + local sigma_p sigma_m + + set _Ftmp = _Stmp / _Stmp[-1] + + threshold -q $1 _Ftmp $ERR_ONESIGMA_POS + $sigma_p = $threshold + threshold -q $1 _Ftmp $ERR_ONESIGMA_MIN + $sigma_m = $threshold + + $SIGMA = ($sigma_p - $sigma_m) / 2.0 + $MEDIAN = ($sigma_p + $sigma_m) / 2.0 + $NPTS = _Stmp[-1] +end