根据经纬度 时间 计算太阳高度角和方位角

% function [hs As] = sunpoint（picName）
%%通过经纬度，当地观测时间计算太阳高度角和方位角；
clear;
clc;
% readDir = ‘d:\我的文档\桌面\1‘;
% readPath = [readDir ‘\*.jpg‘];
% readList = dir（readPath）;
% picName = readList.name;
% NF = str2num（picName（1:4）） 
% Y = str2num（picName（5:6））
% R = str2num（picName（7:8））
% S = str2num（picName（9:10））
% F = str2num（picName（11:12））;
NF  = 2011; Y =3;R = 25;S =7;F=2;%%NFYRS F 分别是观测时的年月日时值和分值
w = 30+35/60+6/3600;  %w为当地纬度（成都w = 30+35/60+6/3600）
D =103;E=59+23/60;%D、 E 分别是观测处经度的度值和分度值（成都D =103;E=59+23/60）
%%（1）求时差Eq
%%时差Eq = 0.0028 - 1.9857sina+ 9.9059sin2a-7.0924cosa - 0.6882cos2a;
%%式中: a为日角其求解公式为:a= 2*pi* t/ 365.2422;
%%t 由3部分组成 即 t = N +N1 - N 0;
N0 = 79.6764 + 0.2422*（ NF - 1985） - floor（（ NF - 1985）/4）; %式中: NF为年份 floor（ X） 为 MAT LAB语言中求出不
%大于 X 的最大整数的标准函数。
%%N 为积日 所谓积日即日期在年内的顺序号 例如 平年12月 31日的积日为 365闰年则为366。
%%当给定年月日分别为 NF  Y R 时 相应的积日N 可通过以下MAT LAB伪代码来实现:
A = NF/4;%%NF为年份，Y为月，R为日
B= A - floor（A） ; 
C= 32.8;
if （Y<=2）
    C= 30.6; 
end
if （ B== 0） & （Y> 2） 
    C= 31.8;