# axc_ul/views.py

# from astropy_healpix import HEALPix does not have an option to
#                                     search for pixels non-inclusively
import healpy as hp
import astropy.units as u
from astropy.coordinates import SkyCoord

import scipy.special as sp

from rest_framework.views import APIView
from rest_framework.response import Response
from rest_framework import status

from django.shortcuts import get_object_or_404

from axc_ul.models import Pixel

from .serializers import PixelSerializer

class PixelDetailView(APIView):
    """
    API endpoint that returns the pixel data (counts, exposure, rate)
    for a given survey number and hpid.
    """
    def get(self, request, hpid):
        # Get the survey using the survey_number field.
        # survey = get_object_or_404(Survey, number=survey_number)
        
        # Get the pixel corresponding to the survey and hpid.
        pixel = get_object_or_404(Pixel, hpid=hpid)
        
        # Serialize the pixel data to JSON.
        serializer = PixelSerializer(pixel)
        return Response(serializer.data, status=status.HTTP_200_OK)


class UpperLimitView(APIView):
    """
    Calculate and return upper limits (UL, CRUL, FLUL) based on provided RA, DEC, and CL parameters.

    Args:
        request: HTTP GET request containing 'ra', 'dec', and 'cl' query parameters.

    Returns:
        Response object with JSON data containing calculated UL, CRUL, and FLUL values.

    Raises:
        ValueError if provided parameters are invalid.
    """
    def get(self, request):
        
        try:
            ra = float(request.query_params.get('ra'))
            dec = float(request.query_params.get('dec'))
            confidence_level = float(request.query_params.get('cl'))
        except (TypeError, ValueError):
            return Response(
                {"error": "Invalud parameters, provide RA, DEC, and CL"},
                status = status.HTTP_400_BAD_REQUEST
            )
        
        # hp = HEALPix(
        #     nside = 4096,
        #     order = 'ring',
        #     frame = 'icrs'
        # )
        
        src_coord = SkyCoord(
            ra, dec, unit = 'deg', frame = 'icrs'
        )
        gal = src_coord.galactic
        src_vec = hp.ang2vec(gal.l.deg, gal.b.deg, lonlat = True)

        aperture_radius = 50 # radius of the aperture in arc seconds
        # HPD ~48 arcseconds
        # 90% ~100 arcseconds
        annulus_inner = 100 # 2 * aperture_radius
        annulus_outer = 200 # 4 * aperture_radius
        
        source_pixel_list = hp.query_disc(
            nside = 4096,
            vec = src_vec,
            inclusive = False,
            nest = False,
            radius = (aperture_radius * u.arcsecond).to(u.radian).value
        )
        
        inner_pixel_list = hp.query_disc(
            nside = 4096,
            vec = src_vec,
            inclusive = False,
            nest = False,
            radius = (annulus_inner * u.arcsecond).to(u.radian).value
        )
        
        outer_pixel_list = hp.query_disc(
            nside = 4096,
            vec = src_vec,
            inclusive = False,
            nest = False,
            radius = (annulus_outer * u.arcsecond).to(u.radian).value
        )
        
        
        annulus_pixel_list = [
            item for item in outer_pixel_list if item not in inner_pixel_list
        ]
                
                
        source_pixels = Pixel.objects.filter(hpid__in = source_pixel_list) 
        
        annulus_pixels = Pixel.objects.filter(hpid__in = annulus_pixel_list)
        
        
        N = sum(obj.counts for obj in source_pixels)
        
        Nnpix = len(source_pixels)
        
        Bcounts = sum(obj.counts for obj in annulus_pixels)
        
        Bnpix = len(annulus_pixels)
        
        B = Bcounts / Bnpix * Nnpix
        
        tsum = sum(obj.exposure for obj in source_pixels)
        
        t = tsum / Nnpix
        
        
        EEF = .9      # eclosed energy fraction, .5 for hpd, .9 for w90
        ECF = 4e-11   # energy conversion factor
        
        UL = (
            sp.gammaincinv(
                N + 1, 
                confidence_level * sp.gammaincc(N + 1, B) + sp.gammainc(N + 1, B)
            ) - B
        )                   # count upper limit
                            
         
        CRUL = UL / t / EEF # count rate upper limit
        
        FLUL = CRUL * ECF   # flux upper limit
        
        return Response({
            'UpperLimit' : UL,
            'CountRateUpperLimit' : CRUL,
            'FluxUpperLimit' : FLUL,
            'N' : N,
            'Nnpix' : Nnpix,
            'Bcounts' : Bcounts,
            'Bnpix' : Bnpix,
            'B' : B,
            'tsum' : tsum,
            't' : t
        }, status=status.HTTP_200_OK)