#!/usr/local/bin/perl
#
# testcombine.pl
#
# Program to create some dithered/rotated frames and combine them
# to see how the combination works.
#
#
# Copyright (C) 2003 Paul A. Price
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#                                                                             
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#                                                                             
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
#
#
# CVS: $Id: mkfakeimages.pl,v 1.2 2004-05-20 00:15:57 price Exp $
#
# History:
# Added orderly star placement.  Fixed so that the star positions are
# accurate once mapped.  --- PAP, 20/10/03.
# Allow different seeing between different images --- PAP 19/5/04.
#


# Parameters.
$DEG_IN_RAD = 57.2957795130823;   # Degrees per radian.

$GAIN = 1;			# Number of electrons per ADU
$READNOISE = 10;		# Detector read noise (electrons);
$NUMSTARS = 500;		# Number of stars on images.
$INSIZE = 1024;		    # Size of input images (pixels on a side).
$OUTSIZE = 1024;	    # Size of output image (pixels on a side).
$SEEING = 0; # Seeing in input images (input pixels).  Zero means use array.
$NUMCRS = 0;			# Number of CRs on input images.
$STARFILE = "stars";	      # File containing list of stars (x,y,f).
$IMAGE = "test";		# Base file name for images and stuff.
$BACKGROUND = 10000;		# Sky background for input images
$OUTSCALE = sqrt(1);  # Pixel scale of output image compared to input.
$STARANGLE = 1.5; # Angle at which a regular star field is tilted (degrees)
$STARFLUXMIN = 500;		# Minimum flux of stars (ADU)
$STARFLUXMAX = 20000;		# Maximum flux of stars (ADU)
$MAXOFFSET = 5;	 # Maximum offset for each image (output image scale).
@ROTATIONS=(0,3);		# Rotations of the images
@SEEING = (2.5, 4.5);	       # Seeing in input images (input pixels)


$CREATE = "./fakeimage -d -x $INSIZE -y $INSIZE -g $GAIN -r $READNOISE -b $BACKGROUND";   # Command to create an input image
$CREATELIST = "-m ";   # Switch to take stars from a list.
$CREATESEEING = "-f ";   # Switch for seeing

# Parse command-line arguments
for ($i=0;$i<=$#ARGV;$i++)
{
    if ($ARGV[$i] eq "-numstars") { $NUMSTARS = $ARGV[++$i]; }
    elsif ($ARGV[$i] eq "-insize") { $INSIZE = $ARGV[++$i]; }
    elsif ($ARGV[$i] eq "-outsize") { $OUTSIZE = $ARGV[++$i]; }
    elsif ($ARGV[$i] eq "-seeing") { $SEEING = $ARGV[++$i]; }
    elsif ($ARGV[$i] eq "-numcrs") { $NUMCRS = $ARGV[++$i]; }
    elsif ($ARGV[$i] eq "-switches") { $COMBINESWITCHES = $ARGV[++$i]; }
    else
    {
	die "Unknown command-line switch: $ARGV[$i]\n";
    }
}



# Generate list of stars.
@x=();   # List of star x positions
@y=();   # List of star y positions
@f=();   # List of star fluxes
open STARS, "> $IMAGE.stars";
$num=0;
for ($i=0;$i<sqrt($NUMSTARS);$i++)
{
    # Random distribution
#    push @x, rand($OUTSIZE);
#    push @y, rand($OUTSIZE);
#    push @f, rand(32*1024-$BACKGROUND);

    # Ordered distribution
    for ($j=0;$j<sqrt($NUMSTARS);$j++)
    {
	push @x, ($i+1)*$OUTSIZE/(sqrt($NUMSTARS)+1)*cos($STARANGLE/$DEG_IN_RAD) + ($j+1)*$OUTSIZE/(sqrt($NUMSTARS)+1)*sin($STARANGLE/$DEG_IN_RAD);
	push @y, - ($i+1)*$OUTSIZE/(sqrt($NUMSTARS)+1)*sin($STARANGLE/$DEG_IN_RAD) + ($j+1)*$OUTSIZE/(sqrt($NUMSTARS)+1)*cos($STARANGLE/$DEG_IN_RAD);
	if ($STARFLUXMIN == $STARFLUXMAX) { push @f, $STARFLUXMAX; }
	else { push @f, rand($STARFLUXMAX-$STARFLUXMIN) + $STARFLUXMIN; }

	print STARS "$x[$num] $y[$num] $f[$num]\n";
	$num++;
    }
}


# Iterate over the input images
for ($i=0;$i<=$#ROTATIONS;$i++)
{
    # Get transformation.
    $B00 = cos($ROTATIONS[$i]/$DEG_IN_RAD) / $OUTSCALE;
    $B01 = sin($ROTATIONS[$i]/$DEG_IN_RAD) / $OUTSCALE;
    $B10 = - $B01;
    $B11 = $B00;
    $A0 = $INSIZE/2;
    $A1 = $INSIZE/2;
    if ($MAXOFFSET == 0)
    {
	$offsetx = 0;
	$offsety = 0;
    }
    else
    {
	$offsetx = rand(2*$MAXOFFSET) - $MAXOFFSET;
	$offsety = rand(2*$MAXOFFSET) - $MAXOFFSET;
    }

    # Output star positions to file
    open STARS, "> ${IMAGE}_$i.stars";
    for ($j=0;$j<$num;$j++)
    {
	$x = $B00*($x[$j] - 1.0 - $OUTSIZE/2 - $offsetx) + $B10*($y[$j] - 1.0 - $OUTSIZE/2 - $offsety) + $A0;
	$y = $B01*($x[$j] - 1.0 - $OUTSIZE/2 - $offsetx) + $B11*($y[$j] - 1.0 - $OUTSIZE/2 - $offsety) + $A1;


	print STARS "$x $y $f[$j]\n";
    }
    close STARS;


    # Create image
    if ($SEEING == 0) { $seeing = $SEEING[$i]; }
    else { $seeing = $SEEING; }
    print "$CREATE $CREATESWITCHES $CREATELIST ${IMAGE}_$i.stars $CREATESEEING $seeing ${IMAGE}_$i.fits\n";
    system "$CREATE $CREATESWITCHES $CREATELIST ${IMAGE}_$i.stars $CREATESEEING $seeing ${IMAGE}_$i.fits";

    # Create map
    $B00 = $OUTSCALE*cos($ROTATIONS[$i]/$DEG_IN_RAD);
    $B01 = - $OUTSCALE*sin($ROTATIONS[$i]/$DEG_IN_RAD);
    $B10 = - $B01;
    $B11 = $B00;

    $A0=$OUTSIZE/2 - ($B00 + $B10)*$INSIZE/2 + $offsetx;
    $A1=$OUTSIZE/2 - ($B01 + $B11)*$INSIZE/2 + $offsety;

    # Output map
    open MAP, "> ${IMAGE}_$i.fits.map";
    print MAP "$A0 $A1\n$B00 $B01\n$B10 $B11\n";
    close MAP;
}

