###### tags: `code`
# Solar dotprod issue

>11/2/2021

## dotprod and diag

## Corripio implementation is good - we just missed the normalisation by zenith angle (eq9 in Tscale 2014)


require(insol)
require(ncdf4)
nc=nc_open("/home/joel/sim/tmapp_points_test/forcing/SURF_201509.nc") 
ssrd =ncvar_get(nc, 'ssrd')
ssrd_wfj = ssrd[2,2,]/3600

startDate="2015-09-01 00:00"
endDate="2015-09-01"
strptime(startDate,format="%Y-%m-%d %H:%M")

start = as.POSIXct(strptime(startDate,format="%Y-%m-%d %H:%M"))

meteodate <- seq(from = start, length.out = 240, by = 3600*3)
metjd = JD(meteodate)
sunv = sunvector(metjd,46.830,9.809,-3)
cos_inc_N = sunv%*%as.vector(normalvector(45,0))
cos_inc_S = sunv%*%as.vector(normalvector(45,180))
theta0 = sunv%*%as.vector(normalvector(0,0))

# here we nomalise subgrid sosine angle bz the grid cosine angle (zenith angle where slope = 0)
IsimN = ssrd_wfj * (cos_inc_N/theta0)
IsimS = ssrd_wfj * (cos_inc_S/theta0)
IsimF = ssrd_wfj * (cos_inc_F/theta0)


plot(meteodate,IsimS, type='l',col='red', ylab='Wm-2', xlab="2015")  
lines(meteodate,IsimN,col='blue')
lines(meteodate,IsimF,col='black')
lines(meteodate,ssrd_wfj ,col='green')
legend("topright", c("south", "north", "flat", "era5"), col =c("red", "blue", "black", "green"), lty=1, bg='white')

> 9/2/2021

Issue 1: I think I multipy by an unnecessary dimension making this multipoint case not scale well

Issue 2: It doesnt work! Solar is reduced strongly on both a N/S 45deg slope.

so heres a minimal example, i have two contrasting north/south 45deg slope points as shown in listpoints:

,id,asp,ele,slp,svf,lon,lat,surfRough,tz,name
0,1,180,2543,45,0.980,9.809,46.830,0.002,0,south
1,2,0,2543,45,0.960,9.847,46.812,0.002,0,north

### Files (by email)
- listpoints:
/home/joel/src/tmapp-evo/tests/testdata/listpoints_sunvtest.txt
- normal vector of slope:
/home/joel/src/tmapp-evo/tests/testdata/nv_sunvtest.npy
- time varying sunvector:
/home/joel/src/tmapp-evo/tests/testdata/sunv_sunvtest.npy

### Single point case that had before when looped over points:

dotprod=np.dot(sunv ,np.transpose(nv))

### Multipoint example we use now:

Dimensions:

In [24]: nv.shape
Out[24]: (2, 3)
points * xyx
In [25]: sunv.shape
Out[25]: (3, 248, 2)   
xyz * time * points

Code:

import numpy as np
nv = np.load("nv_sunvtest.npy")
sunv = np.load("sunv_sunvtest.npy")
dotprod=np.tensordot(sunv  ,np.transpose(nv) , axes=([0],[0]))

this runs but doesnt scale as do unneccessay multiplication on unnecessay dimension - i think this is the issue at least! I drop that dimension (middle one) here:

dprod=dotprod[:,0,:] # drop unneccessary dimension
dprod[dprod<0] = 0 #negative indicates selfshading

this obviously doesnt work as mean dprod is low for both N/S:

In [56]: dprod.mean(0)                                                                                                                                                             
Out[56]: array([0.30212406, 0.13537485])   

other things i looked at but couldnt make work:

a =np.einsum('ijk,jk->ik', sunvr, nv)  

### Docs
https://numpy.org/doc/stable/reference/generated/numpy.dot.html
https://numpy.org/doc/stable/reference/generated/numpy.tensordot.html
https://numpy.org/doc/stable/reference/generated/numpy.einsum.html
einsum tutorial:
https://stackoverflow.com/questions/26089893/understanding-numpys-einsum/33641428#33641428