0
Help us understand the problem. What are the problem?

More than 1 year has passed since last update.

posted at

updated at

TypeScriptで蘇る、太陽電波放射輸送計算プログラム

本記事はCBcloud Advent Calendar 2020の14日目の記事です。

Vue.jsの実用的なテクニックに関しては我らが@dadayamaさんが素敵な記事を書いてくれていますので、そちらをご参照ください。

OpenAPI Generatorを使ったスキーマ駆動開発でAPI仕様と型安全な実装を同時に得るための方法
複雑な画面開発に立ち向かう時に捗るVue Composition APIのTips

はじめに

時は遡ること8年前の2012年、私は大学で太陽物理学の研究室に配属されました。
私は太陽物理の世界に熱中したと同時に、プログラミングの世界にも踏み入れることにもなりました。

そこで出会った言語がIDLです。
IDLはもともと水星探査機の画像解析を目的として作られた言語で、長い歴史を持ちながらも太陽物理学界では現役で使われています。

当時このIDLを利用して画像解析をしていた私ですが、プログラムに関しては全くのド素人でした。
「IDLを使いこなすと画像解析が出来る」以外の知識はほとんどなく、それ以外は少年期にホームページビルダーを触ったことがある程度でした。

大学院に進学しながらも研究者の道は目指さず、Webエンジニアとしてのキャリアを歩み始めて5年になりますが、
ふと当時のプログラムをTypeScriptで書いて、Web上で動かせるようにしたらどうだろう?多分どうにもならないけどやってみます。

蘇らせるもの

当時の私の専門は太陽電波天文学?で、太陽大気中で電波がどのように伝搬するか、
また観測値から太陽の物理量(特に磁場強度)をどのように得るかに関心を割いていました。

そこで今回は「太陽大気中を電波が伝搬した際に、各波長においてどの程度の放射強度で測定できるか」を計算するプログラムを復刻してみます。

当時IDLで書いたプログラムを見てみよう

大学時代に使用していたUSBメモリにプログラムが丸々残っていました。
当時はプラグインが何も入っていない素のvimで、コードフォーマッタも何もない状態で書いていたため、綺麗なコードではありません。
コメントの書き方も独特です。

いきなり光速が4行目あたりに定義されているのもいい味が出ていますね(しかも結局どの処理でも使われていない)。

rad_trans.pro
pro rad_trans,h,t,n,f,param,b=b,l=l,co=co

;============================================================
c=299792458 ;Speed of light[m/s]
wlnth=c*1e2/f ;[cm]

;============================================================
;================ set parameter =============================
sz=(size(h))[1]
if not keyword_set(l) then l=(h[1]-h[0]) & l*=1e5;[cm] thickness of each layer
if not keyword_set(co) then co=1.2 ;heavy ion correction (Gary and Hurford 2004)

;===========================================================
;================= Delta Opacity ===========================
gfact=fltarr(sz)
if min(t) lt 2.0*1e5 then begin
    gfact[where(t lt 2.0*1e5)]=(18.2+alog((t[where(t lt 2.0*1e5)]^1.5)/f))
endif 
if max(t) gt 2.0*1e5 then begin
    gfact[where(t gt 2.0*1e5)]=(24.5+alog(t[where(t gt 2.0*1e5)]/f))
endif
xi=gfact*9.78*1e-3*co
k=((n^2)*xi)/((t^(1.5))*(f)^2) ;abrorption coefficient
dtau=k*L

;===========================================================
;============== transfer equation ==========================
tb=fltarr(sz)
dtb=fltarr(sz)
tau=fltarr(sz)

for i=0,sz-1 do begin
    tau[i]=total(dtau[i:sz-1])
endfor

for i=0,sz-1 do begin
    dtb[i] =t[i]*(1-exp(-dtau[i]))
    if i eq 0 then begin
        tb[i] =dtb[i]
    endif else begin
        tb[i] =dtb[i]+tb[i-1]*exp(-dtau[i])
    endelse
endfor

;===========================================================
;=================== make stracture ========================

param={ h:h,        $
    t:t,        $
    n:n,        $
    l:l,        $
    f:f,        $
    dtb:dtb,    $
    tb:tb,      $
    dtau:dtau,  $
    tau:tau,    $
    iobs:tb[sz-1],  $
    sz:sz       }

;===========================================================
;================== with magnetic field ====================
if keyword_set(b) then begin
    fB=b*2.8*1e6;[Hz]
    kr=((n^2)*xi)/((t^(1.5))*(f-fB)^2)
    kl=((n^2)*xi)/((t^(1.5))*(f+fB)^2)
    dtaur=kr*L
    dtaul=kl*L
    tbr=fltarr(sz)
    tbl=fltarr(sz)
    dtbr=fltarr(sz)
    dtbl=fltarr(sz)
    taur=fltarr(sz)
    taul=fltarr(sz)

    for i=0,sz-1 do begin
        taur[i]=total(dtaur[i:sz-1])
        taul[i]=total(dtaul[i:sz-1])
    endfor

    for i=0,sz-1 do begin
        dtbr[i] =t[i]*(1-exp(-dtaur[i]))
        dtbl[i] =t[i]*(1-exp(-dtaul[i]))
        if i eq 0 then begin
            tbr[i]=dtbr[i]
            tbl[i]=dtbl[i]
        endif else begin
            tbr[i]=dtbr[i]+tbr[i-1]*exp(-dtaur[i])
            tbl[i]=dtbl[i]+tbl[i-1]*exp(-dtaul[i])
        endelse
    endfor

    itot=(tbr+tbl)*0.5
    vtot=(tbr-tbl)*0.5

    param={ h:h,        $
        t:t,        $
        n:n,        $
        l:l,        $
        b:b,        $
        f:f,        $
        dtb:dtb,    $
        dtbr:dtbr,  $
        dtbl:dtbl,  $
        tb:tb,      $
        tbr:tbr,    $
        tbl:tbl,    $
        dtau:dtau,  $
        dtaur:dtaur,    $
        dtaul:dtaul,    $
        tau:tau,    $
        taur:taur,  $
        taul:taul,  $
        iobs:itot[sz-1],$
        vobs:vtot[sz-1],$
        sz:sz       }
endif

end

ご丁寧にPDFのマニュアルまでありました。
スクリーンショット 2020-12-12 000503.png

様々な記憶が蘇ってきました。
それでは約5年前に自分が遺した資料をもとに紐解いて行きましょう。

TypeScriptで書く

具体的にTypeScriptで書き下していく前に、太陽電波の放射強度の計算式を改めて理解する必要があります。

その計算式に関しては私が2016年にAstro Physical Journalに投稿した論文を利用します。
本論文の第2章METHODS OF MAGNETIC FIELD MEASUREMENTSにて詳しい説明がありますので、物好きの方はそちらをご覧ください。
本記事にも和訳した説明を書きながら理解しようと思ったのですが、数式の記述が面倒だったので諦めました。

また今回は簡単のため、磁場の存在は無視します。

光学的厚さ

まずは太陽大気の光学的厚さの式を作っていきましょう。
光学的厚さとは、ある波長で大気を観測した際に、どの程度透明に見えるかを示す指標で、0が透明、1を超えると不透明(光学的に厚い)となります。

IDLのコードでは

;===========================================================
;================= Delta Opacity ===========================
gfact=fltarr(sz)
if min(t) lt 2.0*1e5 then begin
    gfact[where(t lt 2.0*1e5)]=(18.2+alog((t[where(t lt 2.0*1e5)]^1.5)/f))
endif 
if max(t) gt 2.0*1e5 then begin
    gfact[where(t gt 2.0*1e5)]=(24.5+alog(t[where(t gt 2.0*1e5)]/f))
endif
xi=gfact*9.78*1e-3*co
k=((n^2)*xi)/((t^(1.5))*(f)^2) ;abrorption coefficient
dtau=k*L

と記述されていました。
これをTypeScriptで書くと

/**
 * GauntFactor(gff)を返す
 * @param temperature 電子温度 (K)
 * @param frequency 観測周波数 (Hz)
 * @returns GauntFactor(gff)
 */
function calcGauntFactor(temperature: number, frequency: number) {
  if (temperature > 2e5) {
    return 18.2 + Math.log(temperature * (3 / 2)) - Math.log(frequency);
  } else {
    return 24.5 + Math.log(temperature) - Math.log(frequency);
  }
}

/**
 * プラズマパラメータ(ξ)の計算
 * @param temperature 電子温度 (K)
 * @param frequency 観測周波数 (Hz)
 * @param heavyIonCorrection 重イオンの補正係数
 * @returns プラズマパラメータ(ξ)
 */
function calcPlasmaParameter(
  temperature: number,
  frequency: number,
  heavyIonCorrection = 1.2
) {
  return (
    9.78 * 1e-3 * calcGauntFactor(temperature, frequency) * heavyIonCorrection
  );
}

/**
 * エミッションメジャー(ΔEM)を返す
 * 電子は視線方向に対して一様に分布していると仮定する
 * @param electronDensity 電子密度 (cm^-3)
 * @param length 大気の厚さ (cm)
 * @returns エミッションメジャー(ΔEM)
 */
function calcEmissionMeasure(electronDensity: number, length = 1e7) {
  return Math.pow(electronDensity, 2) * length;
}

/**
 * 光学的厚さ(Δτ)を返す
 * @param temperature 電子温度 (K)
 * @param frequency 観測周波数 (Hz)
 * @param electronDensity 電子密度 (cm^-3)
 * @returns 光学的厚さ(Δτ)
 */
function calcOpticalThickness(
  temperature: number,
  frequency: number,
  electronDensity: number
) {
  const plasmaParameter = calcPlasmaParameter(temperature, frequency);
  const emissionMeasure = calcEmissionMeasure(electronDensity);
  return (
    (1.2 * plasmaParameter * emissionMeasure) /
    (Math.pow(temperature, 1.5) * Math.pow(frequency, 2))
  );
}

こんな感じになりました。
各変数の計算をfunctionに切り出したため、多少可読性は上がったでしょうか。

TypeScriptで実装しているので、各引数がnumber型で有ることは保証できていますが、jsdocコメントで記載してあるような単位の保証は出来ていません。この点については後ほど修正していきます。

放射強度

続いて、放射強度の計算です
各層の放射強度は、その層の奥から放射される明るさと、その層自体の明るさの足し算になります。
IDLで書いたコードでは、fltarr(sz)で浮動小数点の配列を作り、そこに各層の計算結果を放り込んでいるみたいです。

;===========================================================
;============== transfer equation ==========================
tb=fltarr(sz)
dtb=fltarr(sz)
tau=fltarr(sz)

for i=0,sz-1 do begin
    tau[i]=total(dtau[i:sz-1])
endfor

for i=0,sz-1 do begin
    dtb[i] =t[i]*(1-exp(-dtau[i]))
    if i eq 0 then begin
        tb[i] =dtb[i]
    endif else begin
        tb[i] =dtb[i]+tb[i-1]*exp(-dtau[i])
    endelse
endfor

とりあえずその層の放射強度だけを計算するようにします。

/**
 * 放射強度(Tb)を返す
 * @param temperature 電子温度 (K)
 * @param frequency 観測周波数 (Hz)
 * @param electronDensity 電子密度 (cm^-3)
 * @param backgroundBrightnessTemperature 背景からの放射強度 (K)
 * @retunrs 放射強度(Tb)
 */
function calcBrightnessTemperature(
  temperature: number,
  frequency: number,
  electronDensity: number,
  backgroundBrightnessTemperature?: number
) {
  const opticalThickness = calcOpticalThickness(
    temperature,
    frequency,
    electronDensity
  );
  const brightnessTemperature =
    temperature * (1 - Math.exp(-1 * opticalThickness));

  if (backgroundBrightnessTemperature == null) {
    return brightnessTemperature;
  } else {
    return (
      brightnessTemperature +
      backgroundBrightnessTemperature * Math.exp(-1 * opticalThickness)
    );
  }
}

これで、視線方向の大気をn層に分割したとき、n番目の太陽大気の放射強度を計算する関数が出来ました。

ここで、上記のcalcBrightnessTemperatureの引数を見てもらうと、各層の電子温度と電子密度が必要なことがわかります。これらのパラメータはどのように算出すれば良いのでしょうか。

今回は、私の論文でも利用したSelhorst2005の太陽大気モデルを利用します。ありがたいことに、高度40000kmまでの各層の温度・電子密度を表で提供してくれているので、GoogleSpreadsheetを利用して値を線形補間しつつ、jsonに変換して利用することにします。

csvからjsonに変換する際には適当なWebサービス(https://csvjson.com/) の力をお借りしました。

変換したjsonはこんな形です

[
  {
    "height": 0,
    "temperature": 6520,
    "density": 76800000000000
  },
  {
    "height": 100,
    "temperature": 5410,
    "density": 9890000000000
  },
  {
    "height": 200,
    "temperature": 4990,
    "density": 3930000000000
  },
  {
    "height": 300,
    "temperature": 4770,
    "density": 1710000000000
  },
]

このjsonを利用し、観測で期待される放射強度を計算します。
IDLのコードでは愚直にfor文で回していましたが、reduceを利用することで完結に記述できます。

/**
 * 放射強度の観測値を返す
 * @param frequency 観測周波数 (Hz)
 * @param atmosphereModel 太陽大気モデル
 * @returns 放射強度の観測値
 */
function calcTotalBrightnesstemperature(
  frequency: number,
  atmosphereModel: Atmosphere
) {
  const brightnessTemperature = atmosphereModel.reduce((previous, current) => {
    return calcBrightnessTemperature(
      current.temperature,
      frequency,
      current.density,
      previous
    );
  }, 0);

  return brightnessTemperature;
}

実行してみる

それでは叩いてみます。
特に意味はないですが、好きなのでVue.jsに載せました。

<template>
  <div id="app">
    <div>観測周波数(GHz)<input type="number" v-model="freqnency" /></div>
    <div>輝度温度(K): {{ brightnessTemperature }}</div>
  </div>
</template>

<script lang="ts">
import { Component, Vue } from "vue-property-decorator";
import { calcTotalBrightnesstemperature } from "./radTrans";
import { atmosphereModel } from "./selhorst2005";

@Component
export default class App extends Vue {
  freqnency = 17;

  get brightnessTemperature() {
    return Math.floor(
      calcTotalBrightnesstemperature(this.freqnency * 1e9, atmosphereModel)
    );
  }
}
</script>

vue.gif

それっぽい結果が出ました!

ここで、野辺山電波ヘリオグラフで利用されている17GHzと34GHzの周波数において太陽静穏領域を観測した際の典型的な輝度温度は、それぞれ約10000K、約9000Kで有ることが知られており、これらの値にもよく一致する結果を得ることが出来ました。

image.png

image.png

math.jsで単位の整合性を担保する

ここまである程度正しい結果が得られたのですが、序盤の方に

TypeScriptで実装しているので、各引数がnumber型で有ることは保証できていますが、jsdocコメントで記載してあるような単位の保証は出来ていません。この点については後ほど修正していきます。

と記載したとおり、全ての引数がnumber型で扱われており、単位の担保が出来ていません。
この状態では、今回の実装に手を入れた際に意図せず単位の整合性が取れなくなってしまう可能性があります。

そこで、今回はJavaScriptの計算ライブラリであるmath.jsを利用してみます。
まずは今回利用しているプロジェクトにmath.jsを導入します。

yarn add mathjs @types/mathjs

今回はTypeScriptで利用するため、@types/mathjsも忘れず導入します。
導入できたら、まずは単位の変換を試してみます。

import { unit } from 'mathjs'

// kmの単位情報を持ったmath.Unit型で長さを保持する
const length = mathjs.unit(1, 'km')
// mに変換して出力する
console.log(length .toNumber('m')) // -> 1000

このように、単位付きの数値を扱うことが出来ます。
一方で、以下のように変換不可能な単位を設定するとエラーになります。

import { unit } from 'mathjs'

// kmの単位情報を持ったmath.Unit型で長さを保持する
const length = unit(1, 'km')
// gに変換して出力する
console.log(length.toNumber('g')) // Error: Units do not match ('g' != '1 km')

これの仕組みを利用し、先程実装した放射強度の計算メソッドたちの引数を全てmath.Unit型で受け取り、不正が単位が入力された場合は実行時エラーで検知できるようにしてみます(本当はTypeScriptのコンパイル時に静的に検知したいのですが、いい方法が思いつきませんでした)。

Unitを利用して書き直すと、各functionは以下のようになりました。

type AtmosphereLayer = {
  height: Unit;
  temperature: Unit;
  density: Unit;
};

/**
 * GauntFactor(gff)を返す
 * @param temperature 電子温度 (K)
 * @param frequency 観測周波数 (Hz)
 * @returns GauntFactor(gff)
 */
function calcGauntFactor(temperature: Unit, frequency: Unit) {
  const t = temperature.toNumber("K");
  const f = frequency.toNumber("Hz");
  if (t > 2e5) {
    return 18.2 + Math.log(t * (3 / 2)) - Math.log(f);
  } else {
    return 24.5 + Math.log(t) - Math.log(f);
  }
}

/**
 * プラズマパラメータ(ξ)の計算
 * @param temperature 電子温度 (K)
 * @param frequency 観測周波数 (Hz)
 * @param heavyIonCorrection 重イオンの補正係数
 * @returns プラズマパラメータ(ξ)
 */
function calcPlasmaParameter(
  temperature: Unit,
  frequency: Unit,
  heavyIonCorrection = 1.2
) {
  return (
    9.78 * 1e-3 * calcGauntFactor(temperature, frequency) * heavyIonCorrection
  );
}

/**
 * エミッションメジャー(ΔEM)を返す
 * 電子は視線方向に対して一様に分布していると仮定する
 * @param electronDensity 電子密度 (cm^-3)
 * @param length 大気の厚さ (cm)
 * @returns エミッションメジャー(ΔEM)
 */
function calcEmissionMeasure(electronDensity: Unit, length?: Unit) {
  const l = length?.toNumber("cm") || unit(10, "km").toNumber("cm");
  const ed = electronDensity.toNumber("cm^3");
  return Math.pow(ed, 2) * l;
}

/**
 * 光学的厚さ(Δτ)を返す
 * @param temperature 電子温度 (K)
 * @param frequency 観測周波数 (Hz)
 * @param electronDensity 電子密度 (cm^-3)
 * @returns 光学的厚さ(Δτ)
 */
function calcOpticalThickness(
  temperature: Unit,
  frequency: Unit,
  electronDensity: Unit
) {
  const plasmaParameter = calcPlasmaParameter(temperature, frequency);
  const emissionMeasure = calcEmissionMeasure(electronDensity);
  const t = temperature.toNumber("K");
  const f = frequency.toNumber("Hz");
  return (
    (1.2 * plasmaParameter * emissionMeasure) /
    (Math.pow(t, 1.5) * Math.pow(f, 2))
  );
}

/**
 * 放射強度(Tb)を返す
 * @param temperature 電子温度 (K)
 * @param frequency 観測周波数 (Hz)
 * @param electronDensity 電子密度 (cm^-3)
 * @param backgroundBrightnessTemperature 背景からの放射強度 (K)
 * @retunrs 放射強度(Tb)
 */
function calcBrightnessTemperature(
  temperature: Unit,
  frequency: Unit,
  electronDensity: Unit,
  backgroundBrightnessTemperature?: Unit
) {
  const opticalThickness = calcOpticalThickness(
    temperature,
    frequency,
    electronDensity
  );
  const t = temperature.toNumber("K");
  const bbt = backgroundBrightnessTemperature?.toNumber("K") || 0;
  const brightnessTemperature = t * (1 - Math.exp(-1 * opticalThickness));

  if (backgroundBrightnessTemperature == null) {
    return brightnessTemperature;
  } else {
    return brightnessTemperature + bbt * Math.exp(-1 * opticalThickness);
  }
}

/**
 * 放射強度の観測値を返す
 * @param frequency 観測周波数 (Hz)
 * @param atmosphereModel 太陽大気モデル
 * @returns 放射強度の観測値
 */
function calcTotalBrightnesstemperature(
  frequency: Unit,
  atmosphereModel: Atmosphere
) {
  const brightnessTemperature = atmosphereModel.reduce((previous, current) => {
    return calcBrightnessTemperature(
      current.temperature,
      frequency,
      current.density,
      unit(previous, "K")
    );
  }, 0);

  return brightnessTemperature;
}

export { calcTotalBrightnesstemperature };

また、大気モデルのjsonに関してもUnit形式に変換した上でexportするようにします。

const atmosphereModelOrigin = [
  {
    "height": 0,
    "temperature": 6520,
    "density": 76800000000000
  },
  {
    "height": 100,
    "temperature": 5410,
    "density": 9890000000000
  },
  {
    "height": 200,
    "temperature": 4990,
    "density": 3930000000000
  },
  {
    "height": 300,
    "temperature": 4770,
    "density": 1710000000000
  },
]

const atmosphereModel = atmosphereModelOrigin.map(val => ({
  height: unit(val.height, "km"),
  temperature: unit(val.temperature, "K"),
  density: unit(val.density, "cm^3")
}));

export { atmosphereModel };

また、上記関数群の呼び出し口のvue.jsのコードは

App.vue
<template>
  <div id="app">
    <div>観測周波数(GHz)<input type="number" v-model="freqnency" /></div>
    <div>輝度温度(K): {{ brightnessTemperature }}</div>
  </div>
</template>

<script lang="ts">
import { Component, Vue } from "vue-property-decorator";
import { calcTotalBrightnesstemperature } from "./radTrans";
import { atmosphereModel } from "./selhorst2005";
import { unit } from "mathjs";

@Component
export default class App extends Vue {
  freqnency = 17;

  get brightnessTemperature() {
    const freqnency = unit(this.freqnency, "GHz");
    return Math.floor(
      calcTotalBrightnesstemperature(freqnency, atmosphereModel)
    );
  }
}
</script>

となります。
このように、TypeScriptを活用して引数の型をmath.jsのUnit型に設定し、各関数内で単位の変換を行うことで、単位に不整合が起きにくいコードにすることが出来ました。

ただしここまでやっても、単位の不整合は実行時エラーでしか検知できません。
そこで、最低限の動作を保証するためのユニットテストJestで書いてみます。

ユニットテスト

JestはJavaScript(TypeScript)のテストを書くためのフレームワークで、Vue.jsなどにも対応しています。

まずはjestを導入します。

yarn add --dev jest @types/jest ts-jest

jest.config.jsは以下のようにし、JestをTypeScriptで実行できるようにします。

jest.config.js
module.exports = {
  roots: ["<rootDir>/src"],
  testMatch: [
    "**/?(*.)+(spec|test).+(ts|tsx|js)"
  ],
  transform: {
    "^.+\\.(ts|tsx)$": "ts-jest"
  }
};

その上で、test/radTrans.spec.tsファイルを作成し、以下のようにテストを書きました。

radTrans.spec.ts
import { unit } from "mathjs";
import { Atmosphere, calcTotalBrightnesstemperature } from "../radTrans";

describe("calcTotalBrightnesstemperature", () => {
  const atmosphereModel = [
    {
      height: 1,
      temperature: 2,
      density: 3
    },
    {
      height: 4,
      temperature: 5,
      density: 6
    }
  ].map(v => ({
    height: unit(v.height, "km"),
    temperature: unit(v.temperature, "K"),
    density: unit(v.density, "cm^3")
  })) as Atmosphere;

  it("正しい単位で引数を指定すると0K以上の値が返る", () => {
    const frequency = unit(10, "GHz");
    const result = calcTotalBrightnesstemperature(frequency, atmosphereModel);
    expect(result).toBeGreaterThan(0);
  });

  it("不正な単位で入力するとエラーになる", () => {
    const frequency = unit(10, "km");
    expect(() => {
      calcTotalBrightnesstemperature(frequency, atmosphereModel);
    }).toThrow();
  });
});

お試し実装なのでテストケースのバリエーションは少ないですが、このような形でユニットテストを記述していけば、リリース前の実行時エラーの検出はある程度出来るかと思います。

また、Vue.js側のコンポーネントテストを記述したい場合はVueTestUtilsを使うと良いですね。

やってみて

昔にIDLで書いたコードをTypeScriptで書き直してみましたがいかがでしたか?

この5年間で私自身のスキルが上がったことももちろんありますが、TypeScriptで記述することで型安全なコードの記述や、単体テストも簡単に書くことが出来ました(もしかしたらIDLでも出来たのかもしれませんが)。

また、math.jsを利用することで単位の整合性チェックも簡単に実装することが出来ました。この辺りはやはりTypeScript(JavaScript)を取り巻くエコシステムが成熟していることの強さですかね。

今までフロントエンドの実装をしてきた中では日付の整合性ぐらいしか単位の意識はしてこなかった気がしますが、math.jsを使うことで今回のような実装が広がることが分かったのは今回の収穫かと思います。

古いコードを別の言語で書き直すと新たな知見が得られるのは楽しい経験でしたので、今後もやっていきたいと思いました。

書いたコード

今回書いたコードは以下のリポジトリにおいてあります。
https://github.com/yamish123/solar-radiation-transfer

それでは。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Sign upLogin
0
Help us understand the problem. What are the problem?