#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Oct 29 18:37:03 2024

@author: sam
"""

import numpy as np



def dot_product(v1,v2):
    dp = 0
    for i in range(len(v1)):
        dp += v1[i]*v2[i]
        # adds to c, the product of the ith elements of each vector
    return dp

def cross_product(v1,v2):
    return np.array([v1[1]*v2[2]-v1[2]*v2[1],v1[2]*v2[0]-v1[0]*v2[2],v1[0]*v2[1]-v1[1]*v2[0]])

def partial_derivative(f,u,v): # using central differences
    h = 1e-9
    r_u = (f(u+h,v) - f(u-h,v))/(2*h) # ∂r/∂u
    r_v = (f(u,v+h) - f(u,v-h))/(2*h) # ∂r/∂v
    return r_u,r_v

def integral(function,surface,u_range,v_range):
# integral function that performs surface and flux integrals
# takes scalar/vector fields as function
# parametrically or explicitly defined surface to integrate over
# u and v parameters range in the form of two lists [u1,u2], [v1,v2]
    
    num_divisions = 100
    # number of subdivisions of the range of u and v
    
    du = (u_range[1] - u_range[0])/num_divisions
    # change in u calculated as the (max value - min value)/number of divisions
    dv = (v_range[1] - v_range[0])/num_divisions
    # change in v calculated as the (max value - min value)/number of divisions

    integral = 0
    # initialise integral sum

    for i in range(num_divisions):
    # loops through segments in u range
    
        for j in range(num_divisions):
        # loops through segments in v range
            
            u = u_range[0] + (i + 0.5) * du
            # for finding the mid point of each du
            v = v_range[0] + (j + 0.5) * dv
            # for finding the mid point of each dv
            
            r = surface(u,v)
            # obtain surface function evaluated at the mid points of each dS
            
            r_u,r_v = partial_derivative(surface,u,v)
            # obtain partial derivatives ∂r/∂u, ∂r/∂v
            
            dS = du * dv
            # small change in surface area
            
            if isinstance(r,np.ndarray):
            # this checks the case that the surface function (r) is defined as an array
            # which implies that it is parametrically defined in the testfile
                
                F = function(r[0],r[1],r[2])
                # substitute function with [r[0],r[1],r[2]] from surface parameterisation
                
                if isinstance(F,np.ndarray):
                # this checks the case that the function (F) is defined as an array
                # which implies that it is a vector field and we are performing a flux integral
                    
                    integral += np.linalg.norm(dot_product(F,cross_product(r_u,r_v))) * dS
                    # evaluates flux integral at each dS and adds it to the sum 'integral'
                    
                else:
                # else the function (F) is not defined as an array
                # therefore it is a scalar field and we are performing a surface integral
                
                    integral += F * np.linalg.norm(cross_product(r_u,r_v)) * dS
                    # linalg.norm calculates the magnitude of the vector, in this case, the magnitude of the cross product
                    
            else:
            # else the surface (r) is not defined as an array
            # therefore it is a scalar function which is explicitly defined in the testfile
                
                F = function(u,v,r)
                # substitute function with [u,v,r] where u = x, v = y, r = g(x,y)
            
                if isinstance(F,np.ndarray):
                # this checks the case that the function (F) is defined as an array
                # which implies that it is a vector field and we are performing a flux integral
                    
                    n = np.array([-r_u,-r_v,1])
                    # obtains normal vectors to the explicit surface using the partial_derivative function
                    # evaluated using central differences
                    nhat = n/np.linalg.norm(n)
                    # obtains unit normal vectors
                    
                    integral += np.linalg.norm(dot_product(F,nhat)) * np.linalg.norm(n) * dS
                    # dot product of F with the unit normal since surface is not parametrically defined
                    
                else:
                # else the function (F) is not defined as an array
                # therefore it is a scalar field and we are performing a surface integral
                
                    integral += F * np.sqrt(1 + r_u**2 + r_v**2) * dS
                    # surface integral formula for explicit surfaces
                    
    return integral