## mkxrsreg.pl by Ken Ebisawa (ebisawa@subaru.gsfc.nasa.gov)
## version 1.0 2004-05-16
if($#ARGV!=2){
    print "\n usage: mkxrsreg.pl Euler1 Euler2 Euler3\n\n";
    print " where Euler1,2,3 are the satellite Euler angles.\n\n";
    print " This script creates a ds9 format XRS pixel region file \n";
    print " (output in STDOUT) for given satellite Euler angles.\n\n";
    print " In ASTRO-E2, Z-axis is the pointing direction, and we adopt\n";
    print " Z-Y-Z rotation for the Euler angle definition.\n";
    print " Consequently, (Euler1, 90.0-Euler2) is the equatorial coordinates\n";
    print " of the pointing direction, and (90-Euler3) gives the roll angle that is \n";
    print " measured from North to DETY axis in counter-clockwise.\n";
    exit(0);
}
$Euler[0]=$ARGV[0];
$Euler[1]=$ARGV[1];
$Euler[2]=$ARGV[2];

$rad2deg = 45.0/atan2(1.0,1.0);
&set_matrix(@Euler,$matrix1,$matrix2,$matrix3);

#Definition of these parameters are needed, so that
#these values are taken from subroutines.
my @outvect1=(0.0, 0.0,0.0);
my @outvect2=(0.0, 0.0,0.0);
my @outvect3=(0.0, 0.0,0.0);

$FOCALLEN=4500.0; # XRS focal length 
for($pixel=0;$pixel<=31;$pixel++){
    if($pixel!=2&&$pixel!=3){
	printf("fk5;polygon(");
	for($corner=0;$corner<=3;$corner++){
	    
# Get the pixel corner positions from the teldef file
# Unit is mm, and center of XRS is (0,0), in
# "look-down" coordinates.
	    $focX =  &get_PixelXcorners($pixel,$corner);
	    $focY =  &get_PixelYcorners($pixel,$corner);
	    
# Determine the direction vector in FOC coordinates
	    $phi =atan2(sqrt($focX**2+$focY**2),$FOCALLEN);
	    $theta= atan2($focY,$focX);
	    $Foc[0]=-sin($phi)*cos($theta);
	    $Foc[1]=-sin($phi)*sin($theta);
	    $Foc[2]= cos($phi);
# Convert the direction vector in the equatorial coordinates
# with Euler rotation
	    &multi_matrix(@Foc,       $matrix1, @outvect1);
	    &multi_matrix(@outvect1,  $matrix2, @outvect2);
	    &multi_matrix(@outvect2,  $matrix3, @outvect3);
	    
# Project the tangential vector in alpha and delta
	    if($outvect3[0]==0.0E0){
		if($outvect3[1]==0.0E0){
		    $alpha = 0.0E0;
		}
		elsif($outvect3[1] > 0.0){
		    $alpha= 90.0E0;
		}
		else{
		    $alpha= -90.0E0;
		}
	    }
	    else{
		$alpha = atan2($outvect3[1],$outvect3[0])*$rad2deg;
	    }
	    if ($alpha<0.0){
		$alpha = $alpha + 360.0E0;
	    }
	    if(abs($outvect3[2]) > 1.0) {
		print "Error!\n";
		print "argument of asin below -1 or above 1\n";
		print "argument of asin',$outvect3[2]<\n";
		exit(0);
	    }
	    $delta = &asin($outvect3[2])*$rad2deg;
# Write the ds9 region file for each pixel
	    printf("%10.6f,%10.6f",$alpha,$delta);
	    if($corner<3){printf(",");}
	}
    printf(")\n");
    }
}

sub get_PixelXcorners{
    local($pixel, $corner);
    $pixel=$_[0];$corner=$_[1];
# This is taken from xrs_teldef_2003-02-15.fit
### for pixel X ###
    if($pixel==20){if($corner==0){return   -1.2720;}
		   if($corner==1){return   -0.6480;}
		   if($corner==2){return   -0.6480;}
		   if($corner==3){return   -1.2720;}
  	       }
    if($pixel==22){if($corner==0){return   -0.6320;}
		   if($corner==1){return   -0.0080;}
		   if($corner==2){return   -0.0080;}
		   if($corner==3){return   -0.6320;}
	       }
    if($pixel==31){if($corner==0){return    0.0080;}
		   if($corner==1){return    0.6320;}
		   if($corner==2){return    0.6320;}
		   if($corner==3){return    0.0080;}
	       }
    if($pixel==29){if($corner==0){return    0.6480;}
		   if($corner==1){return    1.2720;}
		   if($corner==2){return    1.2720;}
		   if($corner==3){return    0.6480;}
	       }
    if($pixel==18){if($corner==0){return   -1.9120;}
		   if($corner==1){return   -1.2880;}
		   if($corner==2){return   -1.2880;}
		   if($corner==3){return   -1.9120;}
	       }
    if($pixel==19){if($corner==0){return   -1.2720;}
		   if($corner==1){return   -0.6480;}
		   if($corner==2){return   -0.6480;}
		   if($corner==3){return   -1.2720;}
	       }
    if($pixel==21){if($corner==0){return   -0.6320;}
		   if($corner==1){return   -0.0080;}
		   if($corner==2){return   -0.0080;}
		   if($corner==3){return   -0.6320;}
	       }
    if($pixel==30){if($corner==0){return    0.0080;}
		   if($corner==1){return    0.6320;}
		   if($corner==2){return    0.6320;}
		   if($corner==3){return    0.0080;}
	       }
    if($pixel==28){if($corner==0){return    0.6480;}
		   if($corner==1){return    1.2720;}
		   if($corner==2){return    1.2720;}
		   if($corner==3){return    0.6480;}
	       }
    if($pixel==27){if($corner==0){return    1.2880;}
		   if($corner==1){return    1.9120;}
		   if($corner==2){return    1.9120;}
		   if($corner==3){return    1.2880;}
	       }
    if($pixel==16){if($corner==0){return   -1.9120;}
		   if($corner==1){return   -1.2880;}
		   if($corner==2){return   -1.2880;}
		   if($corner==3){return   -1.9120;}
	       }
    if($pixel==17){if($corner==0){return   -1.2720;}
		   if($corner==1){return   -0.6480;}
		   if($corner==2){return   -0.6480;}
		   if($corner==3){return   -1.2720;}
	       }
    if($pixel==23){if($corner==0){return   -0.6320;}
		   if($corner==1){return   -0.0080;}
		   if($corner==2){return   -0.0080;}
		   if($corner==3){return   -0.6320;}
	       }
    if($pixel==24){if($corner==0){return    0.0080;}
		   if($corner==1){return    0.6320;}
		   if($corner==2){return    0.6320;}
		   if($corner==3){return    0.0080;}
	       }
    if($pixel==26){if($corner==0){return    0.6480;}
		   if($corner==1){return    1.2720;}
		   if($corner==2){return    1.2720;}
		   if($corner==3){return    0.6480;}
	       }
    if($pixel==25){if($corner==0){return    1.2880;}
		   if($corner==1){return    1.9120;}
		   if($corner==2){return    1.9120;}
		   if($corner==3){return    1.2880;}
	       }
    if($pixel== 9){if($corner==0){return   -1.9120;}
		   if($corner==1){return   -1.2880;}
		   if($corner==2){return   -1.2880;}
		   if($corner==3){return   -1.9120;}
	       }
    if($pixel==10){if($corner==0){return   -1.2720;}
		   if($corner==1){return   -0.6480;}
		   if($corner==2){return   -0.6480;}
		   if($corner==3){return   -1.2720;}
	       }
    if($pixel== 8){if($corner==0){return   -0.6320;}
		   if($corner==1){return   -0.0080;}
		   if($corner==2){return   -0.0080;}
		   if($corner==3){return   -0.6320;}
	       }
    if($pixel== 7){if($corner==0){return    0.0080;}
		   if($corner==1){return    0.6320;}
		   if($corner==2){return    0.6320;}
		   if($corner==3){return    0.0080;}
	       }
    if($pixel== 1){if($corner==0){return    0.6480;}
		   if($corner==1){return    1.2720;}
		   if($corner==2){return    1.2720;}
		   if($corner==3){return    0.6480;}
              }
    if($pixel== 0){if($corner==0){return    1.2880;}
		   if($corner==1){return    1.9120;}
		   if($corner==2){return    1.9120;}
		   if($corner==3){return    1.2880;}
	       }
    if($pixel==11){if($corner==0){return   -1.9120;}
		   if($corner==1){return   -1.2880;}
		   if($corner==2){return   -1.2880;}
		   if($corner==3){return   -1.9120;}
	       }
    if($pixel==12){if($corner==0){return   -1.2720;}
		   if($corner==1){return   -0.6480;}
		   if($corner==2){return   -0.6480;}
		   if($corner==3){return   -1.2720;}
	       }
    if($pixel==14){if($corner==0){return   -0.6320;}
		   if($corner==1){return   -0.0080;}
		   if($corner==2){return   -0.0080;}
		   if($corner==3){return   -0.6320;}
	       }
    if($pixel== 5){if($corner==0){return    0.0080;}
		   if($corner==1){return    0.6320;}
		   if($corner==2){return    0.6320;}
		   if($corner==3){return    0.0080;}
	       }
    if($pixel== 3){if($corner==0){return    0.6480;}
		   if($corner==1){return    1.2720;}
		   if($corner==2){return    1.2720;}
		   if($corner==3){return    0.6480;}
	       }
    if($pixel==13){if($corner==0){return   -1.2720;}
		   if($corner==1){return   -0.6480;}
		   if($corner==2){return   -0.6480;}
		   if($corner==3){return   -1.2720;}
	       }
    if($pixel==15){if($corner==0){return   -0.6320;}
		   if($corner==1){return   -0.0080;}
		   if($corner==2){return   -0.0080;}
		   if($corner==3){return   -0.6320;}
	       }
    if($pixel== 6){if($corner==0){return    0.0080;}
		   if($corner==1){return    0.6320;}
		   if($corner==2){return    0.6320;}
		   if($corner==3){return    0.0080;}
	       }
    if($pixel== 4){if($corner==0){return    0.6480;}
		   if($corner==1){return    1.2720;}
		   if($corner==2){return    1.2720;}
		   if($corner==3){return    0.6480;}
	       }
    if($pixel== 2){if($corner==0){return    4.0440;}
		   if($corner==1){return    4.6680;}
		   if($corner==2){return    4.6680;}
		   if($corner==3){return    4.0440;}
	       }
}
sub get_PixelYcorners{
    local($pixel, $corner);
    $pixel=$_[0];$corner=$_[1];
# This is taken from xrs_teldef_2003-02-15.fit
### for pixel Y ###
    if($pixel==20){if($corner==0){return    1.2880;}
		   if($corner==1){return    1.2880;}
		   if($corner==2){return    1.9120;}
		   if($corner==3){return    1.9120;}
	       }
    if($pixel==22){if($corner==0){return    1.2880;}
		   if($corner==1){return    1.2880;}
		   if($corner==2){return    1.9120;}
		   if($corner==3){return    1.9120;}
	       }
    if($pixel==31){if($corner==0){return    1.2880;}
		   if($corner==1){return    1.2880;}
		   if($corner==2){return    1.9120;}
		   if($corner==3){return    1.9120;}
	       }
    if($pixel==29){if($corner==0){return    1.2880;}
		   if($corner==1){return    1.2880;}
		   if($corner==2){return    1.9120;}
		   if($corner==3){return    1.9120;}
	       }
    if($pixel==18){if($corner==0){return    0.6480;}
		   if($corner==1){return    0.6480;}
		   if($corner==2){return    1.2720;}
		   if($corner==3){return    1.2720;}
	       }
    if($pixel==19){if($corner==0){return    0.6480;}
		   if($corner==1){return    0.6480;}
		   if($corner==2){return    1.2720;}
		   if($corner==3){return    1.2720;}
	       }
    if($pixel==21){if($corner==0){return    0.6480;}
		   if($corner==1){return    0.6480;}
		   if($corner==2){return    1.2720;}
		   if($corner==3){return    1.2720;}
	       }
    if($pixel==30){if($corner==0){return    0.6480;}
		   if($corner==1){return    0.6480;}
		   if($corner==2){return    1.2720;}
		   if($corner==3){return    1.2720;}
	       }
    if($pixel==28){if($corner==0){return    0.6480;}
		   if($corner==1){return    0.6480;}
		   if($corner==2){return    1.2720;}
		   if($corner==3){return    1.2720;}
	       }
    if($pixel==27){if($corner==0){return    0.6480;}
		   if($corner==1){return    0.6480;}
		   if($corner==2){return    1.2720;}
		   if($corner==3){return    1.2720;}
	       }
    if($pixel==16){if($corner==0){return    0.0080;}
		   if($corner==1){return    0.0080;}
		   if($corner==2){return    0.6320;}
		   if($corner==3){return    0.6320;}
	       }
    if($pixel==17){if($corner==0){return    0.0080;}
		   if($corner==1){return    0.0080;}
		   if($corner==2){return    0.6320;}
		   if($corner==3){return    0.6320;}
	       }
    if($pixel==23){if($corner==0){return    0.0080;}
		   if($corner==1){return    0.0080;}
		   if($corner==2){return    0.6320;}
		   if($corner==3){return    0.6320;}
	       }
    if($pixel==24){if($corner==0){return    0.0080;}
		   if($corner==1){return    0.0080;}
		   if($corner==2){return    0.6320;}
		   if($corner==3){return    0.6320;}
	       }
    if($pixel==26){if($corner==0){return    0.0080;}
		   if($corner==1){return    0.0080;}
		   if($corner==2){return    0.6320;}
		   if($corner==3){return    0.6320;}
	       }
    if($pixel==25){if($corner==0){return    0.0080;}
		   if($corner==1){return    0.0080;}
		   if($corner==2){return    0.6320;}
		   if($corner==3){return    0.6320;}
	       }
    if($pixel== 9){if($corner==0){return   -0.6320;}
		   if($corner==1){return   -0.6320;}
		   if($corner==2){return   -0.0080;}
		   if($corner==3){return   -0.0080;}
	       }
    if($pixel==10){if($corner==0){return   -0.6320;}
		   if($corner==1){return   -0.6320;}
		   if($corner==2){return   -0.0080;}
		   if($corner==3){return   -0.0080;}
	       }
    if($pixel== 8){if($corner==0){return   -0.6320;}
		   if($corner==1){return   -0.6320;}
		   if($corner==2){return   -0.0080;}
		   if($corner==3){return   -0.0080;}
	       }
    if($pixel== 7){if($corner==0){return   -0.6320;}
		   if($corner==1){return   -0.6320;}
		   if($corner==2){return   -0.0080;}
		   if($corner==3){return   -0.0080;}
	       }
    if($pixel== 1){if($corner==0){return   -0.6320;}
		   if($corner==1){return   -0.6320;}
		   if($corner==2){return   -0.0080;}
		   if($corner==3){return   -0.0080;}
	       }
    if($pixel== 0){if($corner==0){return   -0.6320;}
		   if($corner==1){return   -0.6320;}
		   if($corner==2){return   -0.0080;}
		   if($corner==3){return   -0.0080;}
	       }
    if($pixel==11){if($corner==0){return   -1.2720;}
		   if($corner==1){return   -1.2720;}
		   if($corner==2){return   -0.6480;}
		   if($corner==3){return   -0.6480;}
	       }
    if($pixel==12){if($corner==0){return   -1.2720;}
		   if($corner==1){return   -1.2720;}
		   if($corner==2){return   -0.6480;}
		   if($corner==3){return   -0.6480;}
	       }
    if($pixel==14){if($corner==0){return   -1.2720;}
		   if($corner==1){return   -1.2720;}
		   if($corner==2){return   -0.6480;}
		   if($corner==3){return   -0.6480;}
	       }
    if($pixel== 5){if($corner==0){return   -1.2720;}
		   if($corner==1){return   -1.2720;}
		   if($corner==2){return   -0.6480;}
		   if($corner==3){return   -0.6480;}
	       }
    if($pixel== 3){if($corner==0){return   -1.2720;}
		   if($corner==1){return   -1.2720;}
		   if($corner==2){return   -0.6480;}
		   if($corner==3){return   -0.6480;}
	       }
    if($pixel==13){if($corner==0){return   -1.9120;}
		   if($corner==1){return   -1.9120;}
		   if($corner==2){return   -1.2880;}
		   if($corner==3){return   -1.2880;}
	       }
    if($pixel==15){if($corner==0){return   -1.9120;}
		   if($corner==1){return   -1.9120;}
		   if($corner==2){return   -1.2880;}
		   if($corner==3){return   -1.2880;}
	       }
    if($pixel== 6){if($corner==0){return   -1.9120;}
		   if($corner==1){return   -1.9120;}
		   if($corner==2){return   -1.2880;}
		   if($corner==3){return   -1.2880;}
	       }
    if($pixel== 4){if($corner==0){return   -1.9120;}
		   if($corner==1){return   -1.9120;}
		   if($corner==2){return   -1.2880;}
		   if($corner==3){return   -1.2880;}
              }
    if($pixel== 2){if($corner==0){return   -4.1870;}
		   if($corner==1){return   -4.1870;}
		   if($corner==2){return   -3.5630;}
		   if($corner==3){return   -3.5630;}
	       }
}

sub set_matrix{
# Set rotation matrix for Euler angles
    my($alpha,$delta,$theta);
    my($matrix1,$matrix2,$matrix3);
    my($dec2rag);
    $dec2rag = atan2(1.0,1.0)/45.0;
    $alpha = $_[0];
    $delta = $_[1];
    $theta = $_[2];
    $alpha = $dec2rag*$alpha;
    $delta = $dec2rag*$delta;
    $theta = $dec2rag*$theta;
    $matrix1->[0][0]= cos($theta);
    $matrix1->[0][1]= -sin($theta);
    $matrix1->[0][2]= 0.0E0;
    $matrix1->[1][0]= sin($theta);
    $matrix1->[1][1]= cos($theta);
    $matrix1->[1][2]= 0.0E0;
    $matrix1->[2][0]= 0.0E0;
    $matrix1->[2][1]= 0.0E0;
    $matrix1->[2][2]= 1.0E0;

    $matrix2->[0][0]= cos($delta);
    $matrix2->[0][1]= 0.0E0;
    $matrix2->[0][2]= sin($delta);
    $matrix2->[1][0]= 0.0E0;
    $matrix2->[1][1]= 1.0E0;
    $matrix2->[1][2]= 0.0E0;
    $matrix2->[2][0]= -sin($delta);
    $matrix2->[2][1]= 0.0E0;
    $matrix2->[2][2]= cos($delta);

    $matrix3->[0][0]= cos($alpha);
    $matrix3->[0][1]= -sin($alpha);
    $matrix3->[0][2]= 0.0E0;
    $matrix3->[1][0]= sin($alpha);
    $matrix3->[1][1]= cos($alpha);
    $matrix3->[1][2]= 0.0E0;
    $matrix3->[2][0]= 0.0E0;
    $matrix3->[2][1]= 0.0E0;
    $matrix3->[2][2]= 1.0E0;

    $_[3]=$matrix1; 
    $_[4]=$matrix2; 
    $_[5]=$matrix3; 
}
sub multi_matrix{
    my(@invect,$matrix,@outvect,$i,$j);
    $invect[0] =$_[0];
    $invect[1] =$_[1];
    $invect[2] =$_[2];
    $matrix    =$_[3];
    for($i=0;$i<=2;$i++){
	$outvect[$i]=0.0;
	for($j=0;$j<=2;$j++){
	    $outvect[$i] = $outvect[$i]+$matrix->[$i][$j]*$invect[$j];
	}
    }
    $_[4]=$outvect[0];
    $_[5]=$outvect[1];
    $_[6]=$outvect[2];
}

sub asin { atan2($_[0], sqrt(1 - $_[0] * $_[0])) }