% /hao/aim3/emery/cetic_climatology/calc_m2g_uam.m
% 11/12 bae Transfer from magnetic to geographic grid, and average over 13.1glat,
%  converting Um and Vm to Un and Vn using the apex declination at the equator

dtr = pi/180.;

% nd = # of days 
nd=61;

dglon = 15.;
lon=[0.:dglon:345.];
nlon=length(lon)
nlonp = nlon+1;

nhlatd=[90. 84.85 80.92 77.82 75.23 72.96 70.87 68.86 66.86 64.84 62.77 60.64 58.44 56.13 53.70 51.14 48.42 45.54 42.48 39.26 35.87 32.34 28.68 24.92 21.07 17.17 13.23 9.27 5.29 0.];
m = 60;
for n=1:29
 m = m - 1;
 latd(m) = -nhlatd(n);
 latd(n) = nhlatd(n);
 if (n>1)
  dlatd(n) = 0.5*(nhlatd(n-1)-nhlatd(n+1));
  dlatd(m) = dlatd(n);
 end
end
 latd(30) = 0.;
 nlatd=length(latd)
 dlatd(30) = latd(29);
 dlatd(1) = nhlatd(1)-nhlatd(2);
 dlatd(nlatd) = dlatd(1);

% datapex=ilon,glon,dinc0,decl0,jlat,glat,alat0,alatm,alat,alatp
% Use glon from 0:15:355 from 0:5:365; decl0, and glat (of mag equ)
   datapex = dlmread('/home/emery/cnofs/out_cnofs_mageq_5deg','',3,0);
   nlon3 = length(datapex)
   j = 0;
   for j3=1:3:nlon3
    j = j + 1;
    nhlat(j) = datapex(j3,6)+13.1;
    shlat(j) = datapex(j3,6)-13.1;
    decl(j) = datapex(j3,4);
   end

   filum=sprintf('/hao/aim3/emery/ccmc_uam/LongRun61_4_UAM/UAM4_07325_08020_375Un')
   filvm=sprintf('/hao/aim3/emery/ccmc_uam/LongRun61_4_UAM/UAM4_07325_08020_375Vn')
  DUM = fopen(filum,'r');
  DVM = fopen(filvm,'r');

   OUTLON = fopen('UAM4_lt_vs_glon_winds','w');
   OUTLT = fopen('UAM4_lt_winds','w');

for iut=1:24
 for j=1:nlon
  numd(iut,j) = 0;
  umd(iut,j) = 0.;
  vmd(iut,j) = 0.;
 end
end

 for md=1:nd
  for nh=1:24
% Skip over 8 headers
   for ns=1:8
    rdhd = fgetl(DUM);
    rdhd = fgetl(DVM);
   end
   for nla=1:nlatd
%  Read lat and values for nlon (25) mlons
    linum = fgetl(DUM);
    inum = sscanf(linum,'%f');
    linvm = fgetl(DVM);
    invm = sscanf(linvm,'%f');
    for j=1:nlon
     if (latd(nla)<nhlat(j) & latd(nla)>shlat(j))
      numd(nh,j) = numd(nh,j) +1;
      umd(nh,j) = umd(nh,j) + inum(j+1);
      vmd(nh,j) = vmd(nh,j) + invm(j+1);
     end
    end
   end
  end	% for nh=1:24
 end	% for md=1:nd
% Average and convert to un and vn geographic coords
  for nh=1:24
   for j=1:nlon
    umd(nh,j) = umd(nh,j)/numd(nh,j);
    vmd(nh,j) = vmd(nh,j)/numd(nh,j);
    und(nh,j) = umd(nh,j)*cos(decl(j)*dtr) + vmd(nh,j)*sin(decl(j)*dtr);
    vnd(nh,j) = vmd(nh,j)*cos(decl(j)*dtr) + umd(nh,j)*sin(decl(j)*dtr);
   end
  end
      
% Convert from UT to LT
for ih=1:24
 for jl=1:nlon
  numk(ih,jl) = 0;
  un(ih,jl) = 0.;
  vn(ih,jl) = 0.;
 end
end
 for ih=1:24
  for jl=1:nlon
   stl = ih + lon(jl)/15.;
   if (stl >= 24.) 
    stl = stl - 24.;
   end
   istl = fix(stl) + 1;
   numk(istl,jl) = numk(istl,jl) + 1;
   un(istl,jl) = un(istl,jl) + und(ih,jl);
   vn(istl,jl) = vn(istl,jl) + vnd(ih,jl);
  end
 end
for ih=1:24
 for jl=1:nlon
  un(ih,jl) = un(ih,jl)/numk(ih,jl);
  vn(ih,jl) = vn(ih,jl)/numk(ih,jl);
 end
end

 fprintf(OUTLON,'LT glon 0-355 of 2 rows of 375 km +/-13.1glat Un then Vn in m/s \n');
lt=0.5:1:23.5;
for istl=1:24
 fprintf(OUTLON,'% 5.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f\n',lt(istl),un(istl,:));
 fprintf(OUTLON,'% 5.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f\n',lt(istl),vn(istl,:));
end

% average over glon
for ih=1:24
 unlt(ih) = 0.;
 vnlt(ih) = 0.;
 for jl=1:nlon
  unlt(ih) = unlt(ih) + un(ih,jl);
  vnlt(ih) = vnlt(ih) + vn(ih,jl);
 end
 unlt(ih) = unlt(ih)/nlon;
 vnlt(ih) = vnlt(ih)/nlon;
end

 fprintf(OUTLT,'LT 375km +/-13.1glat Un Vn in m/s \n');
for ih=1:24
 fprintf(OUTLT,'%5.1f %7.2f %7.2f \n',lt(ih),unlt(ih),vnlt(ih));
end
