3rd commit

.travis.yml

@@ -0,0 +1,24 @@
+language: c
+
+sudo: required
+before_install:
+  - sudo apt-get install gfortran
+
+matrix:
+  include:
+  - os: linux
+    dist: trusty
+    sudo: required
+    env: PLUMED_CC=mpif90 PLUMED_CXX=mpif90  MAKEDOC=yes PLUMED_ALL_TESTS=yes
+  - os: linux
+    dist: trusty
+    sudo: required
+    env: PLUMED_CC=mpif90 PLUMED_CXX=mpif90  PLUMED_CXXFLAGS=-O3 LAPACK=yes CONFIG_FLAGS="--disable-external-lapack"
+
+install:
+  - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then sudo apt-get update -qq ; fi
+  - if test "$PLUMED_CXX" == "mpif90" -a "$TRAVIS_OS_NAME" == "linux" ; then sudo apt-get install -y libopenmpi-dev openmpi-bin ; fi
+  - make
+
+script: 
+  - if test "$PLUMED_CXX" ; then make; fi

README.md

@@ -1,2 +1,32 @@
-# MOCASSIN-3.0
-MOnte CArlo SimulationS of Ionized Nebulae, Version 3.0
+## MOCASSIN-3.0
+[![Build Status](https://travis-ci.org/equib/MOCASSIN-3.0.svg?branch=master)](https://travis-ci.org/equib/MOCASSIN-3.0)
+[![Binder](http://mybinder.org/badge.svg)](http://mybinder.org/repo/equib/mocassin-3.0)
+
+**MOnte CArlo SimulationS of Ionized Nebulae**, Version 3.0
+
+Copyright (C) 2007 Barbara Ercolano 
+
+### Installation
+
+How to compile:
+
+    make mocassin
+    make mocassinWarm
+    make mocassinOutput
+    make mocassinPlot
+
+How to clean:
+
+    make  clean
+     
+How to run:
+
+    mpirun -np N ./mocassin
+
+where N is number of processors for parallel computing.
+
+### References
+
+* Ercolano, B., Clarke, C. J., & Drake, J. J., [ApJ, 699, 1639, 2009](http://adsabs.harvard.edu/abs/2009ApJ...699.1639E)
+
+* Ercolano, B., Drake, J. J., Raymond, J. C., & Clarke, C. C., [ApJ, 688, 398-407, 2008](http://adsabs.harvard.edu/abs/2008ApJ...688..398E)

accessories/dustShellInput.in

@@ -0,0 +1,5 @@
+4.48e15 8.51e17 ! rin, rout
+0.005            ! mdmg
+'mrnsizes.dat'  ! sizefile
+'mrnspecies.dat'! speciesfile
+50 50 50        ! nx ny nz

accessories/makeClothoShell.f90

@@ -0,0 +1,189 @@
+program makeDustShell
+  implicit none 
+
+  real(kind = 8) :: rin, rout, radius, normWeight
+  real(kind = 8) :: dx, dy, dz, den, dx1
+  real(kind = 8) :: mdmg, sums, Pi = 3.1416
+  real(kind = 8), pointer :: xaxis(:), yaxis(:), & 
+       & zaxis(:), rho(:), rad(:), a(:), &
+       & wa(:), abun(:), rhospec(:), da(:)
+
+  integer :: nx, ny, nz, nrho, iskip, irho
+  integer :: i, j, k, nspec, nsizes, ii, ai
+
+  character(len=50) :: sizefile, speciesfile, skip
+
+  open(unit=10, file='dustShellInput.in') 
+
+  read(10, *) rin, rout
+  read(10, *) mdmg
+  read(10, *) sizefile
+  read(10, *) speciesfile
+  read(10, *) nx, ny, nz
+
+  close(10)
+
+  open(unit=10, file='rhoClotho.dat')   
+  read(10,*) nrho
+
+  allocate(rad(nrho))
+  allocate(rho(nrho))
+
+  do i = 1, nrho
+     read(10,*) rad(i), rho(i)
+  end do
+
+  close(10)
+
+  open(unit=10, file=speciesfile)
+  read(10,*) nspec 
+  allocate(abun(nspec))
+  allocate(rhospec(nspec))
+  do i = 1, nspec
+     read(10,*) skip, abun(i), rhospec(i)
+  end do
+  close(10)
+
+  open(unit=10, file=sizefile)
+  read(10,*) nsizes
+  allocate(a(nsizes))
+  allocate(wa(nsizes))
+  allocate(da(nsizes))
+  do i = 1, nsizes
+     read(10,*) iskip, a(i), wa(i)
+  end do
+  close(10)
+
+
+  if (nSizes>1) then
+     da(1) = a(2)-a(1)                    
+     do ai = 2, nSizes-1
+        da(ai) = (a(ai+1)-a(ai-1))/2.              
+     end do
+     da(nSizes) = a(nSizes)-a(nSizes-1)
+  end if
+  normWeight=  0.
+  do ai = 1, nSizes
+     normWeight = normWeight+wa(ai)*da(ai) 
+  end do
+  if (nSizes>1) then
+     do ai = 1, nSizes
+        wa(ai) = (wa(ai)*da(ai))/normWeight
+        if (.not.wa(ai)>=0.) then
+           print*, '! makeDustXSec : Invalid grain weight ', wa(ai), ai
+           stop
+        end if
+     end do
+  end if
+
+  print*, ' index, a [um], da [um], weight '
+  do ai = 1, nSizes
+     print*, ai, real(a(ai)), real(da(ai)), real(wa(ai))
+  end do
+
+  a = a*1.e-4 ! cm
+  sums = 0.
+
+  do i = 1, nspec
+     do j = 1, nsizes
+        sums = sums + rhospec(i)*abun(i)*&
+             & (4./3.)*Pi*wa(j)*a(j)**3
+
+     enddo
+  enddo
+
+  allocate(xaxis(nx+14))
+  allocate(yaxis(ny+14))
+  allocate(zaxis(nz+14))
+
+  xaxis = 0.
+  yaxis = 0.
+  zaxis = 0.
+
+
+!  dx = (rout-rin)/(nx-1)
+!  dy = (rout-rin)/(ny-1)
+!  dz = (rout-rin)/(nz-1)  
+
+  dx = (log10(rout)-log10(rin))/(nx-1)
+  dy = (log10(rout)-log10(rin))/(ny-1)
+  dz = (log10(rout)-log10(rin))/(nz-1)
+
+  dx1 = rin/14.
+
+  do i = 2, 14     
+     xaxis(i) = (xaxis(i-1)+dx1)
+  end do
+  do i = 2, 14
+     yaxis(i) = (yaxis(i-1)+dx1)
+  end do
+  do i = 2, 14  
+     zaxis(i) = (zaxis(i-1)+dx1)
+  end do
+
+  xaxis(2:14) = log10(xaxis(2:14))
+  yaxis(2:14) = log10(yaxis(2:14))
+  zaxis(2:14) = log10(zaxis(2:14))
+
+  do i = 15, nx+14     
+     xaxis(i) = log10(rad(i-14))
+  end do
+  do i = 15, ny+14
+     yaxis(i) = log10(rad(i-14))
+  end do
+  do i = 15, nz+14
+     zaxis(i) = log10(rad(i-14))
+  end do
+
+  xaxis(2:nx+14) = 10.**xaxis(2:nx+14)
+  yaxis(2:nx+14) = 10.**yaxis(2:nx+14)
+  zaxis(2:nx+14) = 10.**zaxis(2:nx+14)
+
+  open(unit=12, file = 'ndust.out') 
+
+  do i = 1, nx+14
+     do j = 1, ny+14
+        do k = 1, nz+14
+
+           radius = 1.e10*sqrt( (xaxis(i)/1.e10)*(xaxis(i)/1.e10) + &
+                & (yaxis(j)/1.e10)*(yaxis(j)/1.e10) + &
+                & (zaxis(k)/1.e10)*(zaxis(k)/1.e10) )
+
+           irho = 0
+
+           do ii = 1, nrho
+
+              if(rad(ii) > radius) then
+                 irho = ii-1
+                 exit
+              end if
+           end do
+           if (irho == 0) then
+              den = 0.
+           elseif (irho ==nrho) then
+              den = rho(nrho)
+           elseif (irho > 0 .and. irho < nrho) then
+              den = rho(irho+1) + (rho(irho)-rho(irho+1))*&
+                   & (rad(irho+1)-radius)/(rad(irho+1)-rad(irho))              
+           else
+              print*, 'insanity in rho location'
+              stop
+           end if           
+
+           den = den*0.7/1.67d-24
+
+           den = den*3.32d-10
+
+
+           write(12, *) real(xaxis(i)), real(yaxis(j)), real(zaxis(k)), real(den)
+
+        end do
+
+     end do
+  end do
+
+  close(12)
+
+!  print*, mdmg/sums, 3.32d-10*0.7/1.67d-24, (3.32d-10*0.7/1.67d-24)/(mdmg/sums)
+
+end program makeDustShell

accessories/mrnsizes.dat

@@ -0,0 +1,32 @@
+          28  sizes 
+           1  5.0000008E-03  1.1313702E+08
+           2  5.7795621E-03  6.8134608E+07
+           3  6.6806665E-03  4.1032760E+07
+           4  7.7222600E-03  2.4711240E+07
+           5  8.9262547E-03  1.4881871E+07
+           6  1.0317967E-02   8962321.    
+           7  1.1926661E-02   5397386.    
+           8  1.3786174E-02   3250475.    
+           9  1.5935604E-02   1957538.    
+          10  1.8420158E-02   1178890.    
+          11  2.1292087E-02   709963.8    
+          12  2.4611777E-02   427562.5    
+          13  2.8449053E-02   257491.4    
+          14  3.2884602E-02   155069.3    
+          15  3.8011711E-02   93387.56    
+          16  4.3938190E-02   56240.86    
+          17  5.0788693E-02   33870.01    
+          18  5.8707271E-02   20397.57    
+          19  6.7860432E-02   12284.06    
+          20  7.8440711E-02   7397.835    
+          21  9.0670548E-02   4455.209    
+          22  0.1048072       2683.064    
+          23  0.1211479       1615.826    
+          24  0.1400363       973.1000    
+          25  0.1618697       586.0309    
+          26  0.1871071       352.9262    
+          27  0.2162794       212.5430    
+          28  0.2500000       128.0001    
+
+
+ slope (p):    3.500000    

accessories/mrnspecies.dat

@@ -0,0 +1,4 @@
+2
+'dustData/sil-dlee.nk' 0.5288 3.6
+'dustData/grph2-dl.nk' 0.47 2.2
+'dustData/grph1-dl.nk' 0.1178 2.2

accessories/rhoClotho.dat

@@ -0,0 +1,51 @@
+50
+  4.4800001E+15  1.3687948E-15
+  4.9863309E+15  9.1459748E-16
+  5.5498880E+15  6.8187640E-16
+  6.1771379E+15  5.3306070E-16
+  6.8752800E+15  4.2759727E-16
+  7.6523261E+15  3.4837793E-16
+  8.5171945E+15  2.8665477E-16
+  9.4798110E+15  2.3737357E-16
+  1.0551222E+16  1.9735550E-16
+  1.1743725E+16  1.6447244E-16
+  1.3071004E+16  1.3722924E-16
+  1.4548293E+16  1.1453396E-16
+  1.6192545E+16  9.5562965E-17
+  1.8022632E+16  7.9677365E-17
+  2.0059557E+16  6.6369333E-17
+  2.2326694E+16  5.5226542E-17
+  2.4850062E+16  4.5908193E-17
+  2.7658625E+16  3.8128579E-17
+  3.0784612E+16  3.1645707E-17
+  3.4263898E+16  2.6253306E-17
+  3.8136415E+16  2.1775053E-17
+  4.2446603E+16  1.8060331E-17
+  4.7243936E+16  1.4980997E-17
+  5.2583461E+16  1.2428760E-17
+  5.8526458E+16  1.0312863E-17
+  6.5141138E+16  8.5578617E-18
+  7.2503412E+16  7.1014266E-18
+  8.0697771E+16  5.8922267E-18
+  8.9818254E+16  4.8880162E-18
+  9.9969539E+16  4.0539758E-18
+  1.1126813E+17  3.3613417E-18
+  1.2384369E+17  2.7862878E-18
+  1.3784054E+17  2.3086652E-18
+  1.5341932E+17  1.9122937E-18
+  1.7075883E+17  1.5836033E-18
+  1.9005804E+17  1.3112292E-18
+  2.1153846E+17  1.0856665E-18
+  2.3544662E+17  8.9897648E-19
+  2.6205686E+17  7.4453751E-19
+  2.9167463E+17  6.1683360E-19
+  3.2463978E+17  5.1127543E-19
+  3.6133065E+17  4.2404928E-19
+  4.0216837E+17  3.5198944E-19
+  4.4762156E+17  2.9247111E-19
+  4.9821191E+17  2.4332022E-19
+  5.5451996E+17  2.0273792E-19
+  6.1719202E+17  1.6923756E-19
+  6.8694724E+17  1.4159198E-19
+  7.6458623E+17  1.1878992E-19
+  8.5100001E+17  9.9999997E-20

benchmarks/dust/1D/p0tau1/input/input.in

@@ -0,0 +1,19 @@
+autoPackets 0.20 2. 10000000
+output
+symmetricXYZ
+contShape  blackbody
+nebComposition noGas
+maxIterateMC  30 95.
+nPhotons 100000
+nx 16
+ny 16
+nz 16
+LStar 38.26
+Rin 0.
+Rout 2.1800001e17
+TStellar 2500.
+Ndust file 'input/p0tau1.ndust'
+dustFile 'input/p0tau1_grainspecies.dat' 'input/p0tau1_grainsizes.dat'
+writeGrid 80. 
+convLimit 0.05
+