profile_read_oc.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4b06b283",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "% function [head, phy, raw] = profile_read(fname, dbin, dtime, ach)\n",
    "%\n",
    "% Reads a Raymetrics/Licel data file applying an analog displacement,\n",
    "% and a dead time correction. Data can be output in both physical or\n",
    "% raw units.\n",
    "%\n",
    "% Input:\n",
    "% \n",
    "%   fname: full path to data file\n",
    "% \n",
    "%    dbin: delay, in number of bins, between analog and PC channels\n",
    "%\n",
    "%   dtime: dead time (in sec) for a correction like like S/(1-S*dtime)\n",
    "%\n",
    "%     ach: read only channel number 'ach', instead of all channels\n",
    "%\n",
    "% Output:\n",
    "%\n",
    "%    head: Is a matlab structure which contains all information in the\n",
    "%          file header. It has two extra fields: the matlab\n",
    "%          julian-day-number of the initial and final dates.\n",
    "%\n",
    "%     phy: Is a matrix with vertical bins as rows, and channels as\n",
    "%          columns. Values are in physical units.\n",
    "%\n",
    "%     raw: Is a matrix with vertical bins as rows, and channels as\n",
    "%          columns. Values are in raw units.\n",
    "%\n",
    "function [head, phy, raw] = profile_read(fname, dbin, dtime, ach, maxz)\n",
    "\n",
    "% if dbin not given, displace by zero #####(~ significa NOT)\n",
    "if ~exist('dbin','var') dbin=0; end\n",
    "if isempty(dbin) dbin=0; end\n",
    "% if dtime not given, no dead time correction\n",
    "if ~exist('dtime','var') dtime=0; end\n",
    "if isempty(dtime) dtime=0; end\n",
    "% if ach not requested, return all channels\n",
    "if ~exist('ach','var') ach=0; end \n",
    "if isempty(ach) ach=0; end \n",
    "% if maxz not requested, return all levels\n",
    "if ~exist('maxz','var') maxz=0; end\n",
    "if isempty(maxz) maxz=0; end\n",
    "\n",
    "%% OPEN FILE\n",
    "fp=fopen(fname,'r');\n",
    "if (fp < 0) %% ##Isso é caso não consiga abrir\n",
    "  error(['fail to open file: ' fname]);\n",
    "end\n",
    "\n",
    "%% READ HEADER LINE #1\n",
    "head.file=fscanf(fp,'%s',[1,1]);\n",
    "\n",
    "%% LINE #2\n",
    "head.site=fscanf(fp,'%s',[1,1]);\n",
    "\n",
    "%% some site might have spaces.\n",
    "%% keep reading until we find something like ??/??/????\n",
    "pos=ftell(fp);\n",
    "[tmp, n]=fscanf(fp,'%s',[1,1]);\n",
    "while (length(tmp)~=10 || tmp(3)~='/' || tmp(6)~='/')\n",
    "  head.site=[head.site ' ' tmp];\n",
    "  pos=ftell(fp);\n",
    "  [tmp,n]=fscanf(fp,'%s',[1,1]);\n",
    "end\n",
    "\n",
    "\n",
    "fseek(fp,pos,'bof');\n",
    "\n",
    "head.datei=fscanf(fp,'%2d/%2d/%4d',[1,3]); %% DD MM YY\n",
    "head.houri=fscanf(fp,'%2d:%2d:%2d',[1,3]); %% hh mn ss\n",
    "head.datef=fscanf(fp,'%2d/%2d/%4d',[1,3]); %% DD MM YY\n",
    "head.hourf=fscanf(fp,'%2d:%2d:%2d',[1,3]); %% hh mn ss\n",
    "head.jdi=datenum([head.datei(3:-1:1) head.houri]);\n",
    "\n",
    "head.jdf=datenum([head.datef(3:-1:1) head.hourf]);\n",
    "\n",
    "\n",
    "head.alt=fscanf(fp,'%d',[1,1]);\n",
    "\n",
    "lon=fscanf(fp,'%s',[1,1]);  lon(lon==',')='.'; head.lon=str2num(lon); \n",
    "\n",
    "lat=fscanf(fp,'%s',[1,1]); lat(lat==',')='.'; head.lat=str2num(lat);\n",
    "\n",
    "% Some old Licel does not include 00, T0 and P0 in this line. So we\n",
    "% need to read the rest of the line, and from that try to read what we\n",
    "% want. \n",
    "tmp=fgets(fp);\n",
    "[A,n]=sscanf(tmp,'%s %s %s %s');\n",
    "\n",
    "if n>=1\n",
    "  zen=A(1); zen(zen==',')='.'; head.zen=str2num(zen);   \n",
    "end\n",
    "if n>=2\n",
    "  idum=A(2); idum(idum==',')='.'; head.idum=str2num(idum);\n",
    "end\n",
    "if n>=3\n",
    "  T0=A(3); T0(T0==',')='.'; head.T0=str2num(T0);\n",
    "end\n",
    "if n>=4\n",
    "  P0=A(4); P0(P0==',')='.'; head.P0=str2num(P0);        \n",
    "end\n",
    "\n",
    "%% LINE #3\n",
    "head.nshoots=fscanf(fp,'%d',[1,1]);\n",
    "head.nhz=fscanf(fp,'%d',[1,1]);\n",
    "head.nshoots2=fscanf(fp,'%d',[1,1]);\n",
    "head.nhz2=fscanf(fp,'%d',[1,1]);\n",
    "head.nch=fscanf(fp,'%d',[1,1]);\n",
    "\n",
    "\n",
    "\n",
    "%% READ CHANNEL LINES\n",
    "for i = 1:head.nch\n",
    "  head.ch(i).active =fscanf(fp,'%d',[1,1]);\n",
    "  head.ch(i).photons=fscanf(fp,'%d',[1,1]);             \n",
    "  head.ch(i).elastic=fscanf(fp,'%d',[1,1]);             \n",
    "  head.ch(i).ndata  =fscanf(fp,'%d 1',[1,1]);           \n",
    "  head.ch(i).pmtv   =fscanf(fp,'%d',[1,1]);             \n",
    "  head.ch(i).binw   =fscanf(fp,'%f',[1,1]);             \n",
    "  head.ch(i).wlen   =fscanf(fp,'%5d.',[1,1]);            \n",
    "  head.ch(i).pol    =fscanf(fp,'%1c 0 0 00 000 ',[1,1]);\n",
    "  head.ch(i).bits   =fscanf(fp,'%d',[1,1]);             \n",
    "  head.ch(i).nshoots=fscanf(fp,'%d',[1,1]);             \n",
    "  head.ch(i).discr  =fscanf(fp,'%f',[1,1]);             \n",
    "  head.ch(i).tr     ={fscanf(fp,'%s',[1,1])};\n",
    "end\n",
    "\n",
    "%% FIND FIRST END OF LINE\n",
    "a=' ';\n",
    "while (a~=10)\n",
    "  a=fread(fp,1,'schar');\n",
    "\n",
    "end\n",
    "\n",
    "\n",
    "for ch = 1:head.nch\n",
    "  nz=head.ch(ch).ndata;\n",
    "  trash=fread(fp,2,'schar');\n",
    "  disp('TRASH AQUI')\n",
    "  disp(trash)\n",
    "  %disp('AQUI')\n",
    "  %disp(nz)                                % Teste\n",
    "  %% READ RAW CHANNELS\n",
    "  tmpraw=fread(fp,nz,'int32');\n",
    "  #disp(tmpraw(end))                       % Teste\n",
    "  #disp('Printa')\n",
    "  %disp('Foi')                             % Teste\n",
    "  if (ch==ach | ach==0)\n",
    "    % conversion factor from raw to physical units\n",
    "    if (head.ch(ch).photons==0)\n",
    "      dScale = head.ch(ch).nshoots*2^head.ch(ch).bits/(head.ch(ch).discr*1.e3);\n",
    "      %disp(head.ch(ch).discr)             % Teste\n",
    "      %disp('Contando!')\n",
    "      %disp(dScale)                        % Teste\n",
    "    else \n",
    "      dScale = head.ch(ch).nshoots/20.;\n",
    "      %disp('Contando!')\n",
    "      %disp(dScale)                        % Teste\n",
    "    end\n",
    "    %disp(class(tmpraw))                   % Teste\n",
    "    tmpphy=tmpraw/dScale;\n",
    "    \n",
    "    tmpphy2 = tmpphy(1+dbin:nz)(end);                      % Teste\n",
    "                               % Teste; Vendo quais sao os numeros pra comparar\n",
    "    % for channel 1 and 3, displace by dbin bins\n",
    "    % for channel 2, 4 and 5, correct dead-time\n",
    "    if (ch==1 | ch==3)\n",
    "      % displace by dbin's\n",
    "      tmpphy(1:nz-dbin) = tmpphy(1+dbin:nz);\n",
    "      \n",
    "      % repeat the last dbin values to keep size of vectors\n",
    "      tmpphy(nz-dbin+1:nz) = tmpphy(nz-dbin+1:nz);\n",
    "    else\n",
    "      % correct for dead-time\n",
    "      tmpphy(1:nz) = tmpphy(1:nz)./(1-tmpphy(1:nz)*dtime);\n",
    "    end\n",
    "    \n",
    "    % copy to final destination\n",
    "    if (maxz==0)\n",
    "      maxz=nz;\n",
    "    else\n",
    "      maxz=min(nz, maxz);\n",
    "      head.ch(ch).ndata=maxz;\n",
    "    end\n",
    "    if (ach==0) \n",
    "      phy(1:maxz, ch) = tmpphy(1:maxz);\n",
    "      raw(1:maxz, ch) = tmpraw(1:maxz);\n",
    "    else\n",
    "      phy(1:maxz, 1) = tmpphy(1:maxz);\n",
    "      raw(1:maxz, 1) = tmpraw(1:maxz);\n",
    "    end\n",
    "  end\n",
    "end\n",
    "%disp(trash)\n",
    "%disp(tmpphy2)\n",
    "% erase uncessary header\n",
    "if ~(ach==0)\n",
    "  head.ch(1).active  = head.ch(ach).active ;\n",
    "  head.ch(1).photons = head.ch(ach).photons;\n",
    "  head.ch(1).elastic = head.ch(ach).elastic;\n",
    "  head.ch(1).ndata   = head.ch(ach).ndata  ;\n",
    "  head.ch(1).pmtv    = head.ch(ach).pmtv   ;\n",
    "  head.ch(1).binw    = head.ch(ach).binw   ;\n",
    "  head.ch(1).wlen    = head.ch(ach).wlen   ;\n",
    "  head.ch(1).pol     = head.ch(ach).pol    ;\n",
    "  head.ch(1).bits    = head.ch(ach).bits   ;\n",
    "  head.ch(1).nshoots = head.ch(ach).nshoots;\n",
    "  head.ch(1).discr   = head.ch(ach).discr  ;\n",
    "  head.ch(1).tr      = head.ch(ach).tr     ;\n",
    "\n",
    "  for i = head.nch:-1:2\n",
    "    head.ch(i).active  = [];             \n",
    "    head.ch(i).photons = [];             \n",
    "    head.ch(i).elastic = [];             \n",
    "    head.ch(i).ndata   = [];           \n",
    "    head.ch(i).pmtv    = [];             \n",
    "    head.ch(i).binw    = [];             \n",
    "    head.ch(i).wlen    = [];            \n",
    "    head.ch(i).pol     = [];\n",
    "    head.ch(i).bits    = [];             \n",
    "    head.ch(i).nshoots = [];             \n",
    "    head.ch(i).discr   = [];             \n",
    "    head.ch(i).tr      = [];\n",
    "  end\n",
    "  head.nch=1;\n",
    "end\n",
    "\n",
    "%% CLOSE FILE\n",
    "fclose(fp);\n",
    "\n",
    "%%\n",
    "\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "839d97cd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TRASH AQUI\n",
      "   13\n",
      "   10\n",
      "warning: Matlab-style short-circuit operation performed for operator |\n",
      "warning: called from\n",
      "    profile_read at line 139 column 5\n",
      "\n",
      "TRASH AQUI\n",
      "   13\n",
      "   10\n",
      "TRASH AQUI\n",
      "   13\n",
      "   10\n",
      "TRASH AQUI\n",
      "   13\n",
      "   10\n",
      "TRASH AQUI\n",
      "   13\n",
      "   10\n",
      "ans =\n",
      "\n",
      "  scalar structure containing the fields:\n",
      "\n",
      "    file = RM1183010.264\n",
      "    site = Embrapa_\n",
      "    datei =\n",
      "\n",
      "         30      8   2011\n",
      "\n",
      "    houri =\n",
      "\n",
      "       10   25   46\n",
      "\n",
      "    datef =\n",
      "\n",
      "         30      8   2011\n",
      "\n",
      "    hourf =\n",
      "\n",
      "       10   26   47\n",
      "\n",
      "    jdi = 7.3475e+05\n",
      "    jdf = 7.3475e+05\n",
      "    alt = 100\n",
      "    lon = -60\n",
      "    lat = -2.9000\n",
      "    zen = [](0x0)\n",
      "    idum = 9\n",
      "    T0 = 0\n",
      "    P0 = 0\n",
      "    nshoots = 600\n",
      "    nhz = 10\n",
      "    nshoots2 = 600\n",
      "    nhz2 = 10\n",
      "    nch = 5\n",
      "    ch =\n",
      "\n",
      "      1x5 struct array containing the fields:\n",
      "\n",
      "        active\n",
      "        photons\n",
      "        elastic\n",
      "        ndata\n",
      "        pmtv\n",
      "        binw\n",
      "        wlen\n",
      "        pol\n",
      "        bits\n",
      "        nshoots\n",
      "        discr\n",
      "        tr\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "profile_read('/Volumes/GoogleDrive/My Drive/Estudos/IC/Henrique/Programs/RMlicelPython/30/RM1183010.264')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "10a499f6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans = 1\n"
     ]
    }
   ],
   "source": [
    "3*2^6/2 == 3*(2^6)/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "3643c2f7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans = 1000\n"
     ]
    }
   ],
   "source": [
    "1.e3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d291377",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Octave",
   "language": "octave",
   "name": "octave"
  },
  "language_info": {
   "file_extension": ".m",
   "help_links": [
    {
     "text": "GNU Octave",
     "url": "https://www.gnu.org/software/octave/support.html"
    },
    {
     "text": "Octave Kernel",
     "url": "https://github.com/Calysto/octave_kernel"
    },
    {
     "text": "MetaKernel Magics",
     "url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
    }
   ],
   "mimetype": "text/x-octave",
   "name": "octave",
   "version": "6.4.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}