% /hao/aim3/emery/cetic_climatology/calc_uam4_cnofs.m
% 11/12 bae Calculate from 61-day hourly avs vs 24 glons +/-13.1glat from UAM4
%  bins to compare with C/NOFS data:  1 day LT vs glon 400-550km; 450,530,610,690,770
%  375,407,443,482,525,(573,625,683,746,816)km n(i) or 407-525, 407-482, 482-573, 573-683, 683, 746-816
%   

% nd = # of days 
nd=61;

dglon = 15.;
glon=[0.:dglon:355.];
nlon=length(glon)
lt=[0.5:1.:23.5];

for ih=1:24
  tnlt(ih) = 0.;
  tilt(ih) = 0.;
  oplt(ih) = 0.;
  for k=1:6
  numk(ih,k) = 0;
  nelt(ih,k) = 0.;
  telt(ih,k) = 0.;
  end
 for jl=1:24
   tn2(ih,jl) = 0.;
   ti2(ih,jl) = 0.;
   op2(ih,jl) = 0.;
  for k=1:6
   ne2(ih,jl,k) = 0.;
   te2(ih,jl,k) = 0.;
  end
   tn(ih,jl) = 0.;
   ti(ih,jl) = 0.;
   op(ih,jl) = 0.;
  for k=1:6
   ne(ih,jl,k) = 0.;
   te(ih,jl,k) = 0.;
  end
  for k=1:10
   num3(ih,jl,k) = 0;
   tn3(ih,jl,k) = 0.;
   ti3(ih,jl,k) = 0.;
   op3(ih,jl,k) = 0.;
   ne3(ih,jl,k) = 0.;
   te3(ih,jl,k) = 0.;
  end
 end
end

   filtn=sprintf('/hao/aim3/emery/ccmc_uam/LongRun61_4_UAM_ave/avtn')
   filti=sprintf('/hao/aim3/emery/ccmc_uam/LongRun61_4_UAM_ave/avti')
   filte=sprintf('/hao/aim3/emery/ccmc_uam/LongRun61_4_UAM_ave/avte')
   filne=sprintf('/hao/aim3/emery/ccmc_uam/LongRun61_4_UAM_ave/avne')
   filop=sprintf('/hao/aim3/emery/ccmc_uam/LongRun61_4_UAM_ave/avop')
   DTN = fopen(filtn,'r');
   DTI = fopen(filti,'r');
   DTE = fopen(filte,'r');
   DNE = fopen(filne,'r');
   DOP = fopen(filop,'r');

   OUTLON = fopen('UAM4_lt_vs_glon','w');
   OUTLT = fopen('UAM4_lt','w');

%  375,407,443,482,525,(573,625,683,746,816)km n(i) or 407-525, 407-482, 482-573, 573-683, 683, 746-816
% Read over 5 alts for neutrals, 10 for ions
for k=1:10
% Skip over 2 headers first alt, 3 headers the next alts
   for ns=1:2
     rdhd = fgetl(DTI);
     rdhd = fgetl(DTE);
     rdhd = fgetl(DNE);
     rdhd = fgetl(DOP);
    if (k<=5)
     rdhd = fgetl(DTN);
    end
   end
  if (k>1)
     rdhd = fgetl(DTI);
     rdhd = fgetl(DTE);
     rdhd = fgetl(DNE);
     rdhd = fgetl(DOP);
     if (k<=5)
     rdhd = fgetl(DTN);
     end
  end
% Read 61 days to get single day values in val2 arrays
 for md=1:nd
  for ih=1:24
     rvti = fgetl(DTI);
     rvte = fgetl(DTE);
     rvne = fgetl(DNE);
     rvop = fgetl(DOP);
     datti = sscanf(rvti,'%f');
     datte = sscanf(rvte,'%f');
     datne = sscanf(rvne,'%f');
     datop = sscanf(rvop,'%f');
     if (k<=5)
     rvtn = fgetl(DTN);
     dattn = sscanf(rvtn,'%f');
     end
   for jl=1:nlon
     num3(ih,jl,k) = num3(ih,jl,k) + 1;
     ti3(ih,jl,k) = ti3(ih,jl,k) + datti(3+jl);
     te3(ih,jl,k) = te3(ih,jl,k) + datte(3+jl);
     ne3(ih,jl,k) = ne3(ih,jl,k) + datne(3+jl);
     op3(ih,jl,k) = op3(ih,jl,k) + datop(3+jl);
     if (k<=5)
     tn3(ih,jl,k) = tn3(ih,jl,k) + dattn(3+jl);
     end
   end
  end
 end	% for md=1:nd
  for ih=1:24
   for jl=1:nlon
     ti3(ih,jl,k) = ti3(ih,jl,k)/nd;
     te3(ih,jl,k) = te3(ih,jl,k)/nd;
     ne3(ih,jl,k) = ne3(ih,jl,k)/nd;
     op3(ih,jl,k) = op3(ih,jl,k)/nd;
     if (k<=5)
     tn3(ih,jl,k) = tn3(ih,jl,k)/nd;
     end
   end
  end
end	% for k=1:10

%  375,407,443,482,525,(573,625,683,746,816)km n(i) or 407-525, 407-482, 482-573, 573-683, 683, 746-816
% Now average in alt
 for ih=1:24
  for jl=1:nlon
   m = 1;
   for k=2:5
    tn2(ih,jl) = tn2(ih,jl) + tn3(ih,jl,k);
    ti2(ih,jl) = ti2(ih,jl) + ti3(ih,jl,k);
    op2(ih,jl) = op2(ih,jl) + op3(ih,jl,k);
    te2(ih,jl,m) = te2(ih,jl,m) + te3(ih,jl,k);
    ne2(ih,jl,m) = ne2(ih,jl,m) + ne3(ih,jl,k);
   end
    tn2(ih,jl) = tn2(ih,jl)/4.;
    ti2(ih,jl) = ti2(ih,jl)/4.;
    op2(ih,jl) = op2(ih,jl)/4.;
    te2(ih,jl,m) = te2(ih,jl,m)/4.;
    ne2(ih,jl,m) = ne2(ih,jl,m)/4.;
   m = 2;
   for k=2:4
    te2(ih,jl,m) = te2(ih,jl,m) + te3(ih,jl,k);
    ne2(ih,jl,m) = ne2(ih,jl,m) + ne3(ih,jl,k);
   end
    te2(ih,jl,m) = te2(ih,jl,m)/3.;
    ne2(ih,jl,m) = ne2(ih,jl,m)/3.;
   m = 3;
   for k=4:6
    te2(ih,jl,m) = te2(ih,jl,m) + te3(ih,jl,k);
    ne2(ih,jl,m) = ne2(ih,jl,m) + ne3(ih,jl,k);
   end
    te2(ih,jl,m) = te2(ih,jl,m)/3.;
    ne2(ih,jl,m) = ne2(ih,jl,m)/3.;
   m = 4;
   for k=6:8
    te2(ih,jl,m) = te2(ih,jl,m) + te3(ih,jl,k);
    ne2(ih,jl,m) = ne2(ih,jl,m) + ne3(ih,jl,k);
   end
    te2(ih,jl,m) = te2(ih,jl,m)/3.;
    ne2(ih,jl,m) = ne2(ih,jl,m)/3.;
   m = 5;
    te2(ih,jl,m) = te3(ih,jl,8);
    ne2(ih,jl,m) = ne3(ih,jl,8);
   m = 6;
   for k=9:10
    te2(ih,jl,m) = te2(ih,jl,m) + te3(ih,jl,k);
    ne2(ih,jl,m) = ne2(ih,jl,m) + ne3(ih,jl,k);
   end
    te2(ih,jl,m) = te2(ih,jl,m)/2.;
    ne2(ih,jl,m) = ne2(ih,jl,m)/2.;
   end	% for jl=1:nlon
  end	% for ih=1:24

% Convert from UT to LT
 for ih=1:24
  for jl=1:nlon
   stl = ih + glon(jl)/15.;
   if (stl >= 24.) 
    stl = stl - 24.;
   end
   istl = fix(stl) + 1;
   m = 1;
   numk(istl,m) = numk(istl,m) + 1;
   tn(istl,jl) = tn(istl,jl) + tn2(ih,jl);
   ti(istl,jl) = ti(istl,jl) + ti2(ih,jl);
   te(istl,jl,m) = te(istl,jl,m) + te2(ih,jl,m);
   op(istl,jl) = op(istl,jl) + op2(ih,jl);
   ne(istl,jl,m) = ne(istl,jl,m) + ne2(ih,jl,m);
   for m=2:6
    numk(istl,m) = numk(istl,m) + 1;
    te(istl,jl,m) = te(istl,jl,m) + te2(ih,jl,m);
    ne(istl,jl,m) = ne(istl,jl,m) + ne2(ih,jl,m);
   end
   end	% for jl=1:nlon
  end	% for ih=1:24
for istl=1:24
 m = 1;
   tn(istl,jl) = tn(istl,jl)/numk(istl,m);
   ti(istl,jl) = ti(istl,jl)/numk(istl,m);
   op(istl,jl) = op(istl,jl)/numk(istl,m);
   te(istl,jl,m) = te(istl,jl,m)/numk(istl,m);
   ne(istl,jl,m) = ne(istl,jl,m)/numk(istl,m);
 for m=2:6
   te(istl,jl,m) = te(istl,jl,m)/numk(istl,m);
   ne(istl,jl,m) = ne(istl,jl,m)/numk(istl,m);
 end
end

 fprintf(OUTLON,'LT glon 0-355 of 15 rows of: 400-550 Tn Ti Op Te Ne; 450km Te Ne; 530km Te Ne; 610km Te Ne; 690km Te Ne; 770km Te Ne \n');

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),tn(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),ti(istl,:));
 fprintf(OUTLON,'%5.1f %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e \n',lt(istl),op(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),te(istl,:,1));
 fprintf(OUTLON,'%5.1f %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e \n',lt(istl),ne(istl,:,1));
 for m=2:6
 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),te(istl,:,m));
 fprintf(OUTLON,'%5.1f %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e \n',lt(istl),ne(istl,:,m));
 end
end

% average over glon
for ih=1:24
 for jl=1:nlon
  tnlt(ih) = tnlt(ih) + tn(ih,jl);
  tilt(ih) = tilt(ih) + ti(ih,jl);
  oplt(ih) = oplt(ih) + op(ih,jl);
  for k=1:6
  nelt(ih,k) = nelt(ih,k)+ne(ih,jl,k);
  telt(ih,k) = telt(ih,k)+te(ih,jl,k);
  end
 end
  tnlt(ih) = tnlt(ih)/nlon;
  tilt(ih) = tilt(ih)/nlon;
  oplt(ih) = oplt(ih)/nlon;
  for k=1:6
  nelt(ih,k) = nelt(ih,k)/nlon;
  telt(ih,k) = telt(ih,k)/nlon;
  end
end	% for ih=1:24

 fprintf(OUTLT,'LT 400-550 Tn Ti Op Te Ne; 450km Te Ne; 530km Te Ne; 610km Te Ne; 690km Te Ne; 770km Te Ne \n');
for ih=1:24
 fprintf(OUTLT,'%5.1f %7.2f %7.2f %e %7.2f %e %7.2f %e %7.2f %e %7.2f %e %7.2f %e %7.2f %e \n',lt(ih),tnlt(ih),tilt(ih),oplt(ih),telt(ih,1),nelt(ih,1),telt(ih,2),nelt(ih,2),telt(ih,3),nelt(ih,3),telt(ih,4),nelt(ih,4),telt(ih,5),nelt(ih,5),telt(ih,6),nelt(ih,6));
end






