text:
//===--- Sudoku/Solver_find.h                                           ---===//
//
// Helper functions, to find available options
//===----------------------------------------------------------------------===//
// list_where_*
// appearance_*: collecting options by appearance count in a given set.
//===----------------------------------------------------------------------===//
#pragma once

#include "Board_Section_traits.h"
#include "Location.h"
#include "Options.h"
#include "Size.h"
#include "Value.h"
#include "traits.h"

#include <gsl/gsl>

#include <array>
#include <vector>

#include <algorithm> // find, find_if
#include <iterator>
#include <type_traits> // is_base_of

#include "Board.fwd.h" // Forward declarations

#include <cassert>

namespace Sudoku
{
//===----------------------------------------------------------------------===//
template<int N, typename SectionT>
[[nodiscard]] auto
	list_where_option(SectionT, Value, ptrdiff_t rep_count = elem_size<N>);

template<int N, typename ItrT>
[[nodiscard]] auto
	list_where_option(ItrT begin, ItrT end, Value, ptrdiff_t rep_count = 0);

template<int N, typename SectionT, typename Options = Options<elem_size<N>>>
[[nodiscard]] auto
	list_where_option(const SectionT, const Options sample) noexcept(true);

template<int N, typename SectionT, typename Options = Options<elem_size<N>>>
[[nodiscard]] auto list_where_equal(SectionT, Options sample) noexcept(true);

template<int N, typename Options = Options<elem_size<N>>>
[[nodiscard]] auto
	list_where_subset(const Board<Options, N>&, const Options) noexcept(true);

template<int N, typename SectionT, typename Options = Options<elem_size<N>>>
[[nodiscard]] auto
	list_where_subset(const SectionT, const Options) noexcept(true);

template<int N, typename SectionT, typename Options = Options<elem_size<N>>>
[[nodiscard]] auto
	list_where_any_option(const SectionT, const Options sample) noexcept(true);

template<int N, typename Options = Options<elem_size<N>>, typename SectionT>
[[nodiscard]] Options appearance_once(SectionT section) noexcept;

template<int N, typename Options = Options<elem_size<N>>, typename InItr_>
[[nodiscard]] Options
	appearance_once(const InItr_ begin, const InItr_ end) noexcept;

template<int N, typename SectionT>
[[nodiscard]] auto appearance_sets(const SectionT section);

template<int N, typename InItr_>
[[nodiscard]] auto appearance_sets(const InItr_ begin, const InItr_ end);

//===----------------------------------------------------------------------===//

// List locations in [section] where [value] is an option
template<int N, typename SectionT>
inline auto list_where_option(
	const SectionT section,
	Value value,
	const ptrdiff_t rep_count /*= elem_size<N>*/)
{
	{
		static_assert(Board_Section::traits::is_Section_v<SectionT>);
		assert(rep_count > 0 && rep_count <= elem_size<N>);
	}
	const auto begin = section.cbegin();
	const auto end   = section.cend();

	return list_where_option<N>(begin, end, value, rep_count);
}

// List locations where [value] is an option
template<int N, typename ItrT>
auto list_where_option(
	const ItrT begin,
	const ItrT end,
	const Value value,
	ptrdiff_t rep_count /*= 0*/)
{
	using Options = Options<elem_size<N>>;
	{
		static_assert(traits::is_input<ItrT>);
		static_assert(
			traits::iterator_to<ItrT, const Options> ||
			traits::iterator_to<ItrT, Options>);
		assert(is_valid<N>(value));

		if (rep_count == 0)
			rep_count = std::distance(begin, end);
		else
			assert(rep_count > 0 && rep_count <= std::distance(begin, end));
	}

	std::vector<Location<N>> locations{};
	locations.reserve(gsl::narrow_cast<size_t>(rep_count));

	const auto check_option = [value](Options O) {
		return is_option(O, value);
	};

#if true
	// slightly faster, it saves one run of find_if on each execution,
	// but rep_count can be too small
	auto itr = begin;
	for (gsl::index i{0}; i < rep_count; ++i)
	{
		itr = std::find_if(itr, end, check_option);
		if (itr == end)
		{ // rep_count too large
			break;
		}
		locations.push_back(itr.location());
		++itr;
	}
	assert( // rep_count not too low
		itr == end || std::find_if(itr, end, check_option) == end);
#else
	for (auto last{std::find_if(begin, end, check_option)}; last != end;
		 last = std::find_if(last, end, check_option))
	{
		locations.push_back(last.location());
		++last;
	}
	assert( // rep_count not too low
		locations.size() <= gsl::narrow_cast<size_t>(rep_count));
#endif
	return locations;
}

// List locations in [section] containing [sample]
template<int N, typename SectionT, typename Options>
auto list_where_option( // NOLINT(bugprone-exception-escape)
	const SectionT section,
	const Options sample) noexcept(true)
{ // vector creation and growth could potentially throw, out of memory.
	{
		static_assert(Board_Section::traits::is_Section_v<SectionT>);
	}
	std::vector<Location<N>> locations{};
	locations.reserve(size_t{elem_size<N>});

	const auto check_option = [sample](Options O) noexcept {
		return sample == shared(O, sample);
	};
	const auto end = section.cend();

	for (auto last{std::find_if(section.cbegin(), end, check_option)};
		 last != end;
		 last = std::find_if(last, end, check_option))
	{
		if (!is_answer_fast(*last))
		{
			locations.push_back(last.location());
		}
		++last;
	}
	return locations;
}

// List locations in [section] equal to [sample]
template<int N, typename SectionT, typename Options>
auto list_where_equal( // NOLINT(bugprone-exception-escape)
	const SectionT section,
	const Options sample) noexcept(true)
{ // vector creation and growth could potentially throw, out of memory.
	{
		static_assert(Board_Section::traits::is_Section_v<SectionT>);
		static_assert(std::is_same_v<typename SectionT::value_type, Options>);
	}
	std::vector<Location<N>> locations{};
	locations.reserve(size_t{elem_size<N>});

	const auto end = section.cend();

	for (auto last{std::find(section.cbegin(), end, sample)}; last != end;
		 last = std::find(last, end, sample))
	{
		locations.push_back(last.location());
		++last;
	}
	return locations;
}

// All locations where available options are a subset of the sample.
template<int N, typename Options>
auto list_where_subset( // NOLINT(bugprone-exception-escape)
	const Board<Options, N>& board,
	const Options sample) noexcept(true)
{ // vector creation and growth could potentially throw, out of memory.
	using Location = Location<N>;
	std::vector<Location> list{};
	list.reserve(sample.count()); // minimum

	for (int i{}; i < full_size<N>; ++i)
	{
		const auto& other = board[Location(i)];
		if (!is_answer_fast(other) && // exclude processed answers
			(other == sample || shared(sample, other) == other))
		{
			list.push_back(Location(i));
		}
	}
	return list;
}

// All locations where available options are a subset of the sample.
template<int N, typename SectionT, typename Options>
auto list_where_subset( // NOLINT(bugprone-exception-escape)
	const SectionT section,
	const Options sample) noexcept(true)
{ // vector creation and growth could potentially throw, out of memory.
	{
		static_assert(Board_Section::traits::is_Section_v<SectionT>);
	}
	using Location = Location<N>;
	std::vector<Location> list{};
	list.reserve(sample.count());

	const auto end = section.cend();
	for (auto itr = section.cbegin(); itr != end; ++itr)
	{
		const auto& other = *itr;
		if (!is_answer_fast(other) && // exclude processed answers
			(other == sample || shared(sample, other) == other))
		{
			list.push_back(itr.location());
		}
	}
	return list;
}

// All locations containing an option present in [sample]
template<int N, typename SectionT, typename Options>
auto list_where_any_option( // NOLINT(bugprone-exception-escape)
	const SectionT section,
	const Options sample) noexcept(true)
{ // vector creation and growth could potentially throw, out of memory.
	{
		static_assert(Board_Section::traits::is_Section_v<SectionT>);
	}
	using Location = Location<N>;
	std::vector<Location> locations{};
	locations.reserve(sample.count_all());

	const auto check_option = [sample](Options O) noexcept {
		return is_answer_fast(O) ? false
								 : !(shared(O, sample).count_all() == 0u);
	};
	const auto end = section.cend();

	for (auto itr{std::find_if(section.cbegin(), end, check_option)};
		 itr != end;
		 itr = std::find_if(itr, end, check_option))
	{
		locations.push_back(itr.location());
		++itr;
	}
	return locations;
}

//===----------------------------------------------------------------------===//

//	returning options collected by appearance count in input-dataset
template<int N, typename InItr_>
auto appearance_sets(const InItr_ begin, const InItr_ end)
{
	// clang-format off
/*
Example illustration
Step 1) Collect
input			worker[n]
				|	[0]			|	[1]			|	[2]			|	[3]
1	100 000 001	|	100	000	001	|	000	000	000	|				|
2	110 100 010	|	110	100	011	|	100	000	000	|	000	000	000	|
3	010 000 011	|	110	100	011	|	110	000	011	|	000	000	000	|	000	000	000
4	000 000 011	|	110	100	011	|	110	000	011	|	000	000	011	|	000	000	000
5	101 100	111	|	111	100	111	|	110	100	011	|	100	000	011	|	000	000	011
6 A	000	001	000	|	111	101	111	|	110	101	011	|	100	001	011	|	000	001	011
7	101	100	110	|	111	101	111	|	111	101	111	|	100	101	011	|	100	001	011
8	100	110	010	|	111	111	111	|	111	101	111	|	100	101	011	|	100	101	011
9	010	000	101	|				|	111	101	111	|	110	101	111	|	100	101	011
				worker[n] contains options appearing more than n times
Step 2) flip
				|	000	000	000	|	000	010	000	|	001	010	000	|	011	010	100
532	419	365		|	==0x		|	==1x		|	<=2x		|	<=3x
				worker[n] contains options appearing n times or less
Step 3) XOR [n-1]
				|	== empty!	|	000	010	000	|	001	000	000	|	010	000	100
				|	==0x		|	==1x		|	==2x		|	==3x
				worker[n] contains options appearing exactly n times
// Q: What about the answer bit?
// A: Always true (==not answered)
//? Working with more or less than [elem_size] input elements?
// Less than 9: possible worker[0] is not empty
// Less than 3: no use for worker[3]
// More than 9: shouldn't be an issue
*/
	// clang-format on

	using Options = Options<elem_size<N>>;
	{
		static_assert(
			traits::iterator_to<InItr_, Options> ||
			traits::iterator_to<InItr_, const Options>);
	}
	// To limit processing time, counting up to N
	constexpr auto max = size_t{N}; // default: (9x9 board) up-to 3 times
	std::array<Options, max + 1> worker{};
	worker.fill(Options(Value{0}));

	// Collect options by appearance count
	// worker[n] contains options appearing more than n times (or answer)
	for (auto elem_itr = begin; elem_itr != end; ++elem_itr)
	{
		if (is_answer(*elem_itr))
		{
			// add answer to all
			for (auto& set : worker)
			{
				set += *elem_itr;
			} // OR
		}
		else
		{
			for (size_t i{max}; i > 0; --i)
			{
				worker.at(i) += (worker.at(i - 1) & *elem_itr); // OR( AND )
			}
			worker[0] += *elem_itr; // OR
		}
	}
	// flip -> worker[n] contains options appearing n times or less
	for (auto& option_set : worker)
	{
		option_set.flip();
	}
	// TODO test: can trigger on partial sections; improve?
	// fails if not all options exist
	assert(worker[0].count_all() == 0);

	// XOR -> worker[n] options appearing n times
	for (size_t i{max}; i > 1; --i)
	{
		worker.at(i).XOR(worker.at(i - 1));
		worker.at(i) += Options(Value{0}); // set not-answered
	}
	return worker;
}

//	returning options collected by appearance count in section
template<int N, typename SectionT>
inline auto appearance_sets(const SectionT section)
{ // referral function, creates iterators
	return appearance_sets<N>(section.cbegin(), section.cend());
}

// return a mask for values with a single appearance
template<int N, typename Options, typename InItr_>
Options appearance_once(const InItr_ begin, const InItr_ end) noexcept
{
	{
		static_assert(traits::is_input<InItr_>);
		static_assert(
			traits::iterator_to<InItr_, Options> ||
			traits::iterator_to<InItr_, const Options>);
	}
	Options sum(Value{0});    // helper all used
	Options worker(Value{0}); // multiple uses OR answer
	for (auto itr = begin; itr != end; ++itr)
	{
		if (is_answer(*itr))
		{
			worker = *itr + worker;
		}
		else
		{
			worker += (sum & *itr);
			sum += *itr;
		}
	}
	return worker.flip(); // multiple uses -> single-use
}

// range-based-for
template<int N, typename Options, typename SectionT>
Options appearance_once(SectionT section) noexcept
{
	{
		static_assert(Board_Section::traits::is_Section_v<SectionT>);
		static_assert(std::is_same_v<typename SectionT::value_type, Options>);
	}
	Options sum(Value{0});    // helper all used
	Options worker(Value{0}); // multiple uses OR answer
	for (const Options& item : section)
	{
		if (is_answer(item))
		{
			worker = item + worker; // OR
		}
		else
		{
			worker += (sum & item); // OR ( AND )
			sum += item;
		}
	}
	return worker.flip(); // multiple uses -> single-use
}

} // namespace Sudoku

max_stars_repo_name: FeodorFitsner/fwkSudoku

max_issues_repo_name: FeodorFitsner/fwkSudoku

max_forks_repo_name: FeodorFitsner/fwkSudoku


text:
/*
FRACTAL - A program growing fractals to benchmark parallelization and drawing
libraries.

Copyright 2009-2021, Javier Burguete Tolosa.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this
  list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``AS IS'' AND ANY EXPRESS
OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
EVENT SHALL Javier Burguete Tolosa OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

/**
 * \file simulator.c
 * \brief Source file to define the windows data and functions.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2009-2021, Javier Burguete Tolosa.
 */
#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <glib.h>
#include <libintl.h>
#include <png.h>
#include <ft2build.h>
#include FT_FREETYPE_H
#include <gtk/gtk.h>
#include <GL/glew.h>
#if HAVE_FREEGLUT
#include <GL/freeglut.h>
#elif HAVE_SDL
#include <SDL.h>
#elif HAVE_GLFW
#include <GLFW/glfw3.h>
extern GLFWwindow *window;
#endif

#include "config.h"
#include "fractal.h"
#include "image.h"
#include "text.h"
#include "graphic.h"
#include "draw.h"
#include "simulator.h"

#if HAVE_SDL
SDL_Event exit_event[1];
#endif

static float phid = -45.;       ///< Horizontal perspective angle (in degrees).
static float thetad = 80.;      ///< Vertical perspective angle (in degrees).

/**
 * Function the set the perspective of a point.
 */
static void
perspective (int x,             ///< Point x-coordinate.
             int y,             ///< Point y-coordinate.
             int z,             ///< Point z-coordinate.
             float *X,          ///< Perspective X-coordinate.
             float *Y)          ///< Perspective Y-coordinate.
{
  float cp, sp, st, ct;
  sincosf (graphic->phi, &sp, &cp);
  sincosf (graphic->theta, &st, &ct);
  *X = x * cp - y * sp;
  *Y = z * st + (y * cp + x * sp) * ct;
}

/**
 * Function to set the view perspective.
 */
static void
set_perspective ()
{
  float k1, k2;
  phid = gtk_range_get_value (GTK_RANGE (dialog_simulator->hscale));
  thetad = gtk_range_get_value (GTK_RANGE (dialog_simulator->vscale));
  graphic->phi = phid * M_PI / 180.;
  graphic->theta = thetad * M_PI / 180.;
  if (fractal_3D)
    {
      perspective (0, 0, 0, &k1, &k2);
      graphic->xmin = k1;
      perspective (length, width, 0, &k1, &k2);
      graphic->xmax = k1;
      perspective (length, 0, 0, &k1, &k2);
      graphic->ymin = k2;
      perspective (0, width, height, &k1, &k2);
      graphic->ymax = k2;
    }
  draw ();
}

/**
 * Function to show an error message
 */
void
show_error (const char *message)        ///< Error message string.
{
  GtkMessageDialog *dlg;
  dlg = (GtkMessageDialog *)
    gtk_message_dialog_new (dialog_simulator->window,
                            GTK_DIALOG_DESTROY_WITH_PARENT,
                            GTK_MESSAGE_ERROR,
                            GTK_BUTTONS_CLOSE, "%s", message);
#if !GTK4
  gtk_widget_show_all (GTK_WIDGET (dlg));
  g_signal_connect_swapped (dlg, "response", G_CALLBACK (gtk_widget_destroy),
                            GTK_WIDGET (dlg));
#else
  gtk_widget_show (GTK_WIDGET (dlg));
  g_signal_connect_swapped (dlg, "response", G_CALLBACK (gtk_window_destroy),
                            GTK_WINDOW (dlg));
#endif
}

/**
 * Function to update the enabled fractal characteristics.
 */
static void
dialog_options_update ()
{
  int i;
  DialogOptions *dlg = dialog_options;
  i = gtk_check_button_get_active (dlg->button_3D);
  gtk_widget_set_sensitive (GTK_WIDGET (dlg->label_length), i);
  gtk_widget_set_sensitive (GTK_WIDGET (dlg->entry_length), i);
  for (i = 0; i < N_RANDOM_SEED_TYPES; ++i)
    if (gtk_check_button_get_active (dlg->array_seeds[i]))
      break;
  i = (i == RANDOM_SEED_TYPE_FIXED) ? 1 : 0;
  gtk_widget_set_sensitive (GTK_WIDGET (dlg->label_seed), i);
  gtk_widget_set_sensitive (GTK_WIDGET (dlg->entry_seed), i);
}

/**
 * Function to close a dialog to set the fractal options.
 */
static void
dialog_options_close (DialogOptions * dlg,      ///< DialogOptions struct.
                      int response_id)  ///< Response identifier.
{
  unsigned int i;
  if (response_id == GTK_RESPONSE_OK)
    {
      fractal_diagonal = gtk_check_button_get_active (dlg->button_diagonal);
      fractal_3D = gtk_check_button_get_active (dlg->button_3D);
      for (i = 0; i < N_FRACTAL_TYPES; ++i)
        if (gtk_check_button_get_active (dlg->array_fractals[i]))
          fractal_type = i;
      length = gtk_spin_button_get_value_as_int (dlg->entry_length);
      width = gtk_spin_button_get_value_as_int (dlg->entry_width);
      height = gtk_spin_button_get_value_as_int (dlg->entry_height);
#if HAVE_GTKGLAREA
      gtk_widget_set_size_request (GTK_WIDGET (dialog_simulator->gl_area),
                                   width, height);
#endif
      random_seed = gtk_spin_button_get_value_as_int (dlg->entry_seed);
      nthreads = gtk_spin_button_get_value_as_int (dlg->entry_nthreads);
      animating = gtk_check_button_get_active (dlg->button_animate);
      for (i = 0; i < N_RANDOM_TYPES; ++i)
        if (gtk_check_button_get_active (dlg->array_algorithms[i]))
          random_algorithm = i;
      for (i = 0; i < N_RANDOM_SEED_TYPES; ++i)
        if (gtk_check_button_get_active (dlg->array_seeds[i]))
          random_seed_type = i;
      medium_start ();
      set_perspective ();
      breaking = 1;
    }
  else if (response_id == GTK_RESPONSE_CANCEL);
  else
    return;
#if !GTK4
  gtk_widget_destroy (GTK_WIDGET (dlg->dialog));
#else
  gtk_window_destroy (GTK_WINDOW (dlg->dialog));
#endif
}

/**
 * Function to create a dialog to set the fractal options.
 */
static void
dialog_options_create ()
{
  unsigned int i;
  const char *array_fractals[N_FRACTAL_TYPES] =
    { _("_Tree"), _("_Forest"), _("_Neuron") };
  const char *array_algorithms[N_RANDOM_TYPES] = {
    "_mt19937",
    "_ranlxs0",
    "r_anlxs1",
    "ra_nlxs2",
    "ran_lxd1",
    "ranlx_d2",
    "ranlu_x",
    "ranlux_389",
    "_cmrg",
    "mr_g",
    "_taus2",
    "g_fsr4"
  };
  const char *array_seeds[N_RANDOM_SEED_TYPES] =
    { _("_Default"), _("_Clock based"), _("_Fixed") };
  DialogOptions *dlg = dialog_options;
#if !GTK4
  GtkContainer *content;
#else
  GtkBox *content;
#endif

  dlg->button_diagonal = (GtkCheckButton *) gtk_check_button_new_with_mnemonic
    (_("_Diagonal movement"));
  dlg->button_3D =
    (GtkCheckButton *) gtk_check_button_new_with_mnemonic ("3_D");
  g_signal_connect (dlg->button_3D, "clicked", dialog_options_update, NULL);
  dlg->label_length = (GtkLabel *) gtk_label_new (_("Length"));
  dlg->label_width = (GtkLabel *) gtk_label_new (_("Width"));
  dlg->label_height = (GtkLabel *) gtk_label_new (_("Height"));
  dlg->label_seed = (GtkLabel *) gtk_label_new (_("Random seed"));
  dlg->label_nthreads = (GtkLabel *) gtk_label_new (_("Threads number"));
  dlg->entry_length =
    (GtkSpinButton *) gtk_spin_button_new_with_range (32., 2400., 1.);
  dlg->entry_width =
    (GtkSpinButton *) gtk_spin_button_new_with_range (32., 2400., 1.);
  dlg->entry_height =
    (GtkSpinButton *) gtk_spin_button_new_with_range (20., 2400., 1.);
  dlg->entry_seed =
    (GtkSpinButton *) gtk_spin_button_new_with_range (0., 4294967295., 1.);
  dlg->entry_nthreads =
    (GtkSpinButton *) gtk_spin_button_new_with_range (0., 64., 1.);

  dlg->grid_fractal = (GtkGrid *) gtk_grid_new ();
  dlg->array_fractals[0] = NULL;
  for (i = 0; i < N_FRACTAL_TYPES; ++i)
    {
#if !GTK4
      dlg->array_fractals[i] =
        (GtkRadioButton *) gtk_radio_button_new_with_mnemonic_from_widget
        (dlg->array_fractals[0], array_fractals[i]);
#else
      dlg->array_fractals[i] = (GtkCheckButton *)
        gtk_check_button_new_with_mnemonic (array_fractals[i]);
      gtk_check_button_set_group (dlg->array_fractals[i],
                                  dlg->array_fractals[0]);
#endif
      gtk_grid_attach (dlg->grid_fractal, GTK_WIDGET (dlg->array_fractals[i]),
                       0, i, 1, 1);
    }
  gtk_check_button_set_active (dlg->array_fractals[fractal_type], 1);
  dlg->frame_fractal = (GtkFrame *) gtk_frame_new (_("Fractal type"));
  gtk_frame_set_child (dlg->frame_fractal, GTK_WIDGET (dlg->grid_fractal));

  dlg->button_animate = (GtkCheckButton *) gtk_check_button_new_with_mnemonic
    (_("_Animate"));
  gtk_check_button_set_active (dlg->button_animate, animating);

  dlg->grid_algorithm = (GtkGrid *) gtk_grid_new ();
  dlg->array_algorithms[0] = NULL;
  for (i = 0; i < N_RANDOM_TYPES; ++i)
    {
#if !GTK4
      dlg->array_algorithms[i] =
        (GtkRadioButton *) gtk_radio_button_new_with_mnemonic_from_widget
        (dlg->array_algorithms[0], array_algorithms[i]);
#else
      dlg->array_algorithms[i] = (GtkCheckButton *)
        gtk_check_button_new_with_mnemonic (array_algorithms[i]);
      gtk_check_button_set_group (dlg->array_algorithms[i],
                                  dlg->array_algorithms[0]);
#endif
      gtk_grid_attach (dlg->grid_algorithm,
                       GTK_WIDGET (dlg->array_algorithms[i]), 0, i, 1, 1);
    }
  gtk_check_button_set_active (dlg->array_algorithms[random_algorithm], 1);
  dlg->frame_algorithm = (GtkFrame *) gtk_frame_new (_("Random algorithm"));
  gtk_frame_set_child (dlg->frame_algorithm, GTK_WIDGET (dlg->grid_algorithm));

  dlg->grid_seed = (GtkGrid *) gtk_grid_new ();
  dlg->array_seeds[0] = NULL;
  for (i = 0; i < N_RANDOM_SEED_TYPES; ++i)
    {
#if !GTK4
      dlg->array_seeds[i] =
        (GtkRadioButton *) gtk_radio_button_new_with_mnemonic_from_widget
        (dlg->array_seeds[0], array_seeds[i]);
#else
      dlg->array_seeds[i] = (GtkCheckButton *)
        gtk_check_button_new_with_mnemonic (array_seeds[i]);
      gtk_check_button_set_group (dlg->array_seeds[i], dlg->array_seeds[0]);
#endif
      gtk_grid_attach (dlg->grid_seed, GTK_WIDGET (dlg->array_seeds[i]),
                       0, i, 1, 1);
      g_signal_connect (dlg->array_seeds[i], "clicked", dialog_options_update,
                        NULL);
    }
  gtk_check_button_set_active (dlg->array_seeds[random_seed_type], 1);
  dlg->frame_seed = (GtkFrame *) gtk_frame_new (_("Random seed type"));
  gtk_frame_set_child (dlg->frame_seed, GTK_WIDGET (dlg->grid_seed));

  dlg->grid = (GtkGrid *) gtk_grid_new ();
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->button_diagonal), 0, 0, 2, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->button_3D), 0, 1, 2, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->label_length), 0, 2, 1, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->entry_length), 1, 2, 1, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->label_width), 0, 3, 1, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->entry_width), 1, 3, 1, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->label_height), 0, 4, 1, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->entry_height), 1, 4, 1, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->frame_fractal), 0, 5, 2, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->button_animate), 0, 6, 2, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->frame_algorithm), 2, 0, 1, 8);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->frame_seed), 0, 7, 2, 2);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->label_seed), 0, 9, 1, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->entry_seed), 1, 9, 1, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->label_nthreads), 0, 10, 1, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->entry_nthreads), 1, 10, 1, 1);

  dlg->logo = (GtkImage *) gtk_image_new_from_file ("logo2.png");

  dlg->bar = (GtkHeaderBar *) gtk_header_bar_new ();
#if !GTK4
  gtk_header_bar_set_title (dlg->bar, _("Options"));
  gtk_header_bar_set_subtitle (dlg->bar, _("Set the fractal options"));
  gtk_header_bar_set_show_close_button (dlg->bar, 1);
#endif
  gtk_header_bar_pack_start (dlg->bar, GTK_WIDGET (dlg->logo));

  dlg->dialog
    = (GtkDialog *) gtk_dialog_new_with_buttons (NULL,
                                                 dialog_simulator->window,
                                                 GTK_DIALOG_MODAL |
                                                 GTK_DIALOG_DESTROY_WITH_PARENT,
                                                 _("_OK"),
                                                 GTK_RESPONSE_OK,
                                                 _("_Cancel"),
                                                 GTK_RESPONSE_CANCEL, NULL);
  gtk_window_set_titlebar (GTK_WINDOW (dlg->dialog), GTK_WIDGET (dlg->bar));
#if !GTK4
  content = (GtkContainer *) gtk_dialog_get_content_area (dlg->dialog);
  gtk_container_add (content, GTK_WIDGET (dlg->grid));
  gtk_widget_show_all (GTK_WIDGET (dlg->dialog));
#else
  gtk_window_set_title (GTK_WINDOW (dlg->dialog), _("Options"));
  content = (GtkBox *) gtk_dialog_get_content_area (dlg->dialog);
  gtk_box_append (content, GTK_WIDGET (dlg->grid));
  gtk_widget_show (GTK_WIDGET (dlg->dialog));
#endif

  gtk_check_button_set_active (dlg->button_diagonal, fractal_diagonal);
  gtk_check_button_set_active (dlg->button_3D, fractal_3D);
  gtk_spin_button_set_value (dlg->entry_length, length);
  gtk_spin_button_set_value (dlg->entry_width, width);
  gtk_spin_button_set_value (dlg->entry_height, height);
  gtk_spin_button_set_value (dlg->entry_seed, random_seed);
  gtk_spin_button_set_value (dlg->entry_nthreads, nthreads);
  dialog_options_update ();

  g_signal_connect_swapped (dlg->dialog, "response",
                            G_CALLBACK (dialog_options_close), dlg);
}

/**
 * Function to show a help dialog.
 */
static void
dialog_simulator_help ()
{
  gchar *authors[] = {
    "Javier Burguete Tolosa (jburguete@eead.csic.es)",
    NULL
  };
  gtk_show_about_dialog (dialog_simulator->window,
                         "program_name", "Fractal",
                         "comments",
                         _("A program using growing fractals to benchmark "
                           "parallelization and drawing libraries"),
                         "authors", authors,
                         "translator-credits",
                         _("Javier Burguete Tolosa (jburguete@eead.csic.es)"),
                         "version", "3.4.15",
                         "copyright",
                         "Copyright 2009-2021 Javier Burguete Tolosa",
                         "license-type", GTK_LICENSE_BSD,
                         "logo", dialog_simulator->logo,
                         "website-label", _("Website"),
                         "website", "https://github.com/jburguete/fractal",
                         NULL);
}

/**
 * Function to update the enabled window properties.
 */
void
dialog_simulator_update ()
{
  gtk_widget_set_sensitive
    (GTK_WIDGET (dialog_simulator->button_start), !simulating);
  gtk_widget_set_sensitive
    (GTK_WIDGET (dialog_simulator->button_stop), simulating);
  set_perspective ();
}

/**
 * Function to show the fractal progress.
 */
void
dialog_simulator_progress ()
{
  GMainContext *context;
  register unsigned int k;
  register float x;
  if (fractal_3D)
    {
      switch (fractal_type)
        {
        case 0:
        case 1:
          x = max_d / (float) (height - 1);
          break;
        default:
          k = length;
          if (width < k)
            k = width;
          if (height < k)
            k = height;
          k = k / 2;
          if (k)
            --k;
          x = fmin (1., max_d / (float) k);
        }
    }
  else
    {
      switch (fractal_type)
        {
        case 0:
        case 1:
          x = max_d / (float) (height - 1);
          break;
        default:
          k = width;
          if (height < k)
            k = height;
          k = k / 2 - 1;
          x = fmin (1., max_d / (float) k);
        }
    }
  gtk_progress_bar_set_fraction (dialog_simulator->progress, x);
  gtk_spin_button_set_value
    (dialog_simulator->entry_time, difftime (time (NULL), t0));
  context = g_main_context_default ();
  while (g_main_context_pending (context))
    g_main_context_iteration (context, 0);
}

/**
 * Function to close the dialog to save the graphical view.
 */
static void
dialog_simulator_graphic_close (GtkFileChooserDialog * dlg,
                                ///< GtkFileChooserDialog struct.
                                int response_id)       ///< Response identifier.
{
  char *filename;
  if (response_id == GTK_RESPONSE_ACCEPT)
    filename = gtk_file_chooser_get_current_name (GTK_FILE_CHOOSER (dlg));
  else
    filename = NULL;
#if !GTK4
  gtk_widget_destroy (GTK_WIDGET (dlg));
#else
  gtk_window_destroy (GTK_WINDOW (dlg));
#endif
  if (filename)
    {
      graphic_save (filename);
      g_free (filename);
    }
}

/**
 * Function to save the graphical view on a PNG file.
 */
static void
dialog_simulator_graphic_save ()
{
  GtkFileChooserDialog *dlg;
  dlg = (GtkFileChooserDialog *)
    gtk_file_chooser_dialog_new (_("Save graphical"),
                                 dialog_simulator->window,
                                 GTK_FILE_CHOOSER_ACTION_SAVE,
                                 _("_Cancel"), GTK_RESPONSE_CANCEL,
                                 _("_Open"), GTK_RESPONSE_ACCEPT, NULL);
#if !GTK4
  gtk_file_chooser_set_do_overwrite_confirmation (GTK_FILE_CHOOSER (dlg), 1);
  gtk_widget_show_all (GTK_WIDGET (dlg));
#else
  gtk_widget_show (GTK_WIDGET (dlg));
#endif
  g_signal_connect (dlg, "response",
                    G_CALLBACK (dialog_simulator_graphic_close), NULL);
}

#if HAVE_GTKGLAREA

static void
dialog_simulator_draw_init (GtkGLArea * gl_area)
{
#if DEBUG
  printf ("dialog_simulator_draw_init: start\n");
#endif
  gtk_gl_area_make_current (gl_area);
  graphic_init (graphic, "logo.png");
#if DEBUG
  printf ("dialog_simulator_draw_init: end\n");
#endif
}

#elif HAVE_GLFW

/**
 * Function to close the GLFW window.
 */
static void
window_close ()
{
  glfwSetWindowShouldClose (window, GL_TRUE);
}

#endif

/**
 * Function to create the main window.
 */
void
dialog_simulator_create ()
{
  static char *tip_exit, *tip_options, *tip_start, *tip_stop, *tip_save,
    *tip_help;
  DialogSimulator *dlg;

#if DEBUG
  printf ("dialog_simulator_create: start\n");
#endif

  dlg = dialog_simulator;
#if HAVE_GTKGLAREA
  dlg->loop = g_main_loop_new (NULL, 0);
#elif HAVE_SDL
  exit_event->type = SDL_QUIT;
#endif

  tip_options = _("Fractal options");
  tip_start = _("Start fractal growing");
  tip_stop = _("Stop fractal growing");
  tip_save = _("Save graphical");
  tip_help = _("Help");
  tip_exit = _("Exit");

#if !GTK4
  dlg->logo = gtk_image_get_pixbuf
    (GTK_IMAGE (gtk_image_new_from_file ("logo.png")));
  dlg->logo_min = gtk_image_get_pixbuf
    (GTK_IMAGE (gtk_image_new_from_file ("logo2.png")));
#else
  dlg->logo = gtk_image_get_paintable
    (GTK_IMAGE (gtk_image_new_from_file ("logo.png")));
  dlg->logo_min = gtk_image_get_paintable
    (GTK_IMAGE (gtk_image_new_from_file ("logo2.png")));
#endif

  dlg->box = (GtkBox *) gtk_box_new (GTK_ORIENTATION_HORIZONTAL, 0);

  dlg->button_options = (GtkButton *)
#if !GTK4
    gtk_button_new_from_icon_name ("preferences-system", GTK_ICON_SIZE_BUTTON);
#else
    gtk_button_new_from_icon_name ("preferences-system");
#endif
  gtk_widget_set_tooltip_text (GTK_WIDGET (dlg->button_options), tip_options);
  gtk_box_append (dlg->box, GTK_WIDGET (dlg->button_options));
  g_signal_connect_swapped
    (dlg->button_options, "clicked", dialog_options_create, dlg->logo_min);

  dlg->button_start = (GtkButton *)
#if !GTK4
    gtk_button_new_from_icon_name ("system-run", GTK_ICON_SIZE_BUTTON);
#else
    gtk_button_new_from_icon_name ("system-run");
#endif
  gtk_widget_set_tooltip_text (GTK_WIDGET (dlg->button_start), tip_start);
  gtk_box_append (dlg->box, GTK_WIDGET (dlg->button_start));
  g_signal_connect (dlg->button_start, "clicked", fractal, NULL);

  dlg->button_stop = (GtkButton *)
#if !GTK4
    gtk_button_new_from_icon_name ("process-stop", GTK_ICON_SIZE_BUTTON);
#else
    gtk_button_new_from_icon_name ("process-stop");
#endif
  gtk_widget_set_tooltip_text (GTK_WIDGET (dlg->button_stop), tip_stop);
  gtk_box_append (dlg->box, GTK_WIDGET (dlg->button_stop));
  g_signal_connect (dlg->button_stop, "clicked", fractal_stop, NULL);

  dlg->button_save = (GtkButton *)
#if !GTK4
    gtk_button_new_from_icon_name ("document-save", GTK_ICON_SIZE_BUTTON);
#else
    gtk_button_new_from_icon_name ("document-save");
#endif
  gtk_widget_set_tooltip_text (GTK_WIDGET (dlg->button_save), tip_save);
  gtk_box_append (dlg->box, GTK_WIDGET (dlg->button_save));
  g_signal_connect (dlg->button_save, "clicked", dialog_simulator_graphic_save,
                    NULL);

  dlg->button_help = (GtkButton *)
#if !GTK4
    gtk_button_new_from_icon_name ("help-about", GTK_ICON_SIZE_BUTTON);
#else
    gtk_button_new_from_icon_name ("help-about");
#endif
  gtk_widget_set_tooltip_text (GTK_WIDGET (dlg->button_help), tip_help);
  gtk_box_append (dlg->box, GTK_WIDGET (dlg->button_help));
  g_signal_connect (dlg->button_help, "clicked", dialog_simulator_help, NULL);

  dlg->button_exit = (GtkButton *)
#if !GTK4
    gtk_button_new_from_icon_name ("application-exit", GTK_ICON_SIZE_BUTTON);
#else
    gtk_button_new_from_icon_name ("application-exit");
#endif
  gtk_widget_set_tooltip_text (GTK_WIDGET (dlg->button_exit), tip_exit);
  gtk_box_append (dlg->box, GTK_WIDGET (dlg->button_exit));
#if HAVE_GTKGLAREA
  g_signal_connect_swapped (dlg->button_exit, "clicked",
                            (GCallback) g_main_loop_quit, dlg->loop);
#elif HAVE_FREEGLUT
  g_signal_connect (dlg->button_exit, "clicked", glutLeaveMainLoop, NULL);
#elif HAVE_SDL
  g_signal_connect_swapped (dlg->button_exit, "clicked",
                            (void (*)) SDL_PushEvent, exit_event);
#elif HAVE_GLFW
  g_signal_connect (dlg->button_exit, "clicked", (void (*)) window_close, NULL);
#endif

  dlg->label_time = (GtkLabel *) gtk_label_new (_("Calculating time"));
  dlg->entry_time =
    (GtkSpinButton *) gtk_spin_button_new_with_range (0., 1.e6, 0.1);
  gtk_widget_set_sensitive (GTK_WIDGET (dlg->entry_time), 0);

  dlg->progress = (GtkProgressBar *) gtk_progress_bar_new ();
  gtk_progress_bar_set_text (dlg->progress, _("Progress"));

  dlg->hscale = (GtkScale *) gtk_scale_new_with_range
    (GTK_ORIENTATION_HORIZONTAL, -90., 0., 1.);
  dlg->vscale = (GtkScale *) gtk_scale_new_with_range
    (GTK_ORIENTATION_HORIZONTAL, 0., 90., 1.);
  gtk_scale_set_digits (dlg->hscale, 0);
  gtk_scale_set_digits (dlg->vscale, 0);
  gtk_range_set_value (GTK_RANGE (dlg->hscale), phid);
  gtk_range_set_value (GTK_RANGE (dlg->vscale), thetad);
  g_signal_connect (dlg->hscale, "value-changed", set_perspective, NULL);
  g_signal_connect (dlg->vscale, "value-changed", set_perspective, NULL);

  dlg->label_horizontal = (GtkLabel *) gtk_label_new
    (_("Horizontal perspective angle (º)"));
  dlg->label_vertical = (GtkLabel *) gtk_label_new
    (_("Vertical perspective angle (º)"));

  dlg->grid = (GtkGrid *) gtk_grid_new ();
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->box), 0, 0, 3, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->progress), 0, 1, 1, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->label_time), 1, 1, 1, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->entry_time), 2, 1, 1, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->label_horizontal), 0, 2, 1, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->hscale), 1, 2, 2, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->label_vertical), 0, 3, 1, 1);
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->vscale), 1, 3, 2, 1);

#if HAVE_GTKGLAREA
  dlg->gl_area = (GtkGLArea *) gtk_gl_area_new ();
  gtk_widget_set_size_request (GTK_WIDGET (dlg->gl_area),
                               window_width, window_height);
  g_signal_connect (dlg->gl_area, "realize",
                    (GCallback) dialog_simulator_draw_init, NULL);
  g_signal_connect (dlg->gl_area, "render", draw, NULL);
  g_signal_connect (dlg->gl_area, "resize", (GCallback) resize, NULL);
#if WINDOW_GLAREA
#if !GTK4
  dlg->window_gl = (GtkWindow *) gtk_window_new (GTK_WINDOW_TOPLEVEL);
  gtk_window_set_icon (dlg->window_gl, dlg->logo_min);
  gtk_container_add (GTK_CONTAINER (dlg->window_gl), GTK_WIDGET (dlg->gl_area));
  g_signal_connect_swapped (dlg->window_gl, "delete_event",
                            (GCallback) g_main_loop_quit, dlg->loop);
  gtk_widget_show_all (GTK_WIDGET (dlg->window_gl));
#else
  dlg->window_gl = (GtkWindow *) gtk_window_new ();
  gtk_window_set_child (dlg->window_gl, GTK_WIDGET (dlg->gl_area));
  g_signal_connect_swapped (dlg->window_gl, "close-request",
                            (GCallback) g_main_loop_quit, dlg->loop);
  gtk_widget_show (GTK_WIDGET (dlg->window_gl));
#endif
#else
  gtk_grid_attach (dlg->grid, GTK_WIDGET (dlg->gl_area), 0, 4, 3, 1);
#endif
#elif HAVE_SDL
  SDL_SetWindowMinimumSize (window, window_width, window_height);
#elif HAVE_GLFW
  glfwSetWindowSizeLimits (window, window_width, window_height, GLFW_DONT_CARE,
                           GLFW_DONT_CARE);
#endif

#if !GTK4
  dlg->window = (GtkWindow *) gtk_window_new (GTK_WINDOW_TOPLEVEL);
  gtk_window_set_title (dlg->window, _("Fractal growing"));
  gtk_window_set_icon (dlg->window, dlg->logo_min);
  gtk_container_add (GTK_CONTAINER (dlg->window), GTK_WIDGET (dlg->grid));
#if HAVE_GTKGLAREA
  g_signal_connect_swapped (dlg->window, "delete_event",
                            (GCallback) g_main_loop_quit, dlg->loop);
#elif HAVE_FREEGLUT
  g_signal_connect (dlg->window, "delete_event", glutLeaveMainLoop, NULL);
#elif HAVE_SDL
  g_signal_connect_swapped (dlg->window, "delete_event",
                            (void (*)) SDL_PushEvent, exit_event);
#elif HAVE_GLFW
  g_signal_connect (dlg->window, "delete_event", (void (*)) window_close, NULL);
#endif
  gtk_widget_show_all (GTK_WIDGET (dlg->window));
#else
  dlg->window = (GtkWindow *) gtk_window_new ();
  gtk_window_set_title (dlg->window, _("Fractal growing"));
  gtk_window_set_child (dlg->window, GTK_WIDGET (dlg->grid));
#if HAVE_GTKGLAREA
  g_signal_connect_swapped (dlg->window, "close-request",
                            (GCallback) g_main_loop_quit, dlg->loop);
#elif HAVE_FREEGLUT
  g_signal_connect (dlg->window, "close-request", glutLeaveMainLoop, NULL);
#elif HAVE_SDL
  g_signal_connect_swapped (dlg->window, "close-request",
                            (void (*)) SDL_PushEvent, exit_event);
#elif HAVE_GLFW
  g_signal_connect (dlg->window, "close-request", (void (*)) window_close,
                    NULL);
#endif
  gtk_widget_show (GTK_WIDGET (dlg->window));
#endif

#if DEBUG
  printf ("dialog_simulator_create: end\n");
#endif
}

ext: c

max_stars_repo_name: jburguete/fractal

max_issues_repo_name: jburguete/fractal

max_forks_repo_name: jburguete/fractal


text:
/*
 * BRAINS
 * (B)LR (R)everberation-mapping (A)nalysis (I)n AGNs with (N)ested (S)ampling
 * Yan-Rong Li, liyanrong@ihep.ac.cn
 * Thu, Aug 4, 2016
 */

/*!
 *  \file blr_models.c
 *  \brief BLR models
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <string.h>
#include <mpi.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_interp.h>

#include "brains.h"

/* cluster around outer disk face, Lopn_cos1 < Lopn_cos2 */
inline double theta_sample_outer(double gam, double Lopn_cos1, double Lopn_cos2)
{
  return acos(Lopn_cos1 + (Lopn_cos2-Lopn_cos1) * pow(gsl_rng_uniform(gsl_r), gam));
}

/* cluster around equatorial plane, Lopn_cos1 < Lopn_cos2 */
inline double theta_sample_inner(double gam, double Lopn_cos1, double Lopn_cos2)
{
  /* note that cosine is a decreaing function */
  double opn1 = acos(Lopn_cos1), opn2 = acos(Lopn_cos2);
  return opn2 + (opn1-opn2)  * (1.0 - pow(gsl_rng_uniform(gsl_r), 1.0/gam));
}

/*================================================================
 * model 1
 * Brewer et al. (2011)'s model
 *
 * geometry: radial Gamma distribution
 * dynamics: elliptical orbits
 *================================================================
 */
void gen_cloud_sample_model1(const void *pm, int flag_type, int flag_save)
{
  int i, j, nc;
  double r, phi, dis, Lopn_cos, u;
  double x, y, z, xb, yb, zb, zb0, vx, vy, vz, vxb, vyb, vzb;
  double inc, F, beta, mu, k, a, s, rin, sig, Rs;
  double Lphi, Lthe, L, E;
  double V, weight, rnd;
  BLRmodel1 *model = (BLRmodel1 *)pm;
  double Emin, Lmax, Vr, Vr2, Vph, mbh, chi, lambda, q;

  Lopn_cos = cos(model->opn*PI/180.0);
  inc = acos(model->inc);
  beta = model->beta;
  F = model->F;
  mu = exp(model->mu);
  k = model->k;

  if(flag_type == 1) // 1D RM, no dynamical parameters
  {
    mbh = 0.0;
  }
  else
  {
    mbh = exp(model->mbh);
    lambda = model->lambda;
    q = model->q;
  }

  Rs = 3.0e11*mbh/CM_PER_LD; // Schwarzchild radius in a unit of light-days
  a = 1.0/beta/beta;
  s = mu/a;
  rin = mu*F + Rs;
  sig = (1.0-F)*s;

  for(i=0; i<parset.n_cloud_per_task; i++)
  {
// generate a direction of the angular momentum
    Lphi = 2.0*PI * gsl_rng_uniform(gsl_r);
    //Lthe = acos(Lopn_cos + (1.0-Lopn_cos) * gsl_rng_uniform(gsl_r));
    Lthe = theta_sample(1.0, Lopn_cos, 1.0);

    nc = 0;
    r = rcloud_max_set+1.0;
    while(r>rcloud_max_set || r<rcloud_min_set)
    {
      if(nc > 1000)
      {
        printf("# Error, too many tries in generating ridial location of clouds.\n");
        exit(0);
      }
      rnd = gsl_ran_gamma(gsl_r, a, 1.0);
//    r = mu * F + (1.0-F) * gsl_ran_gamma(gsl_r, 1.0/beta/beta, beta*beta*mu);
        r = rin + sig * rnd;
      nc++;
    }
    phi = 2.0*PI * gsl_rng_uniform(gsl_r);

    x = r * cos(phi);
    y = r * sin(phi);
    z = 0.0;

  /*xb =  cos(Lthe)*cos(Lphi) * x + sin(Lphi) * y - sin(Lthe)*cos(Lphi) * z;
    yb = -cos(Lthe)*sin(Lphi) * x + cos(Lphi) * y + sin(Lthe)*sin(Lphi) * z;
    zb =  sin(Lthe) * x + cos(Lthe) * z;*/

    xb =  cos(Lthe)*cos(Lphi) * x + sin(Lphi) * y;
    yb = -cos(Lthe)*sin(Lphi) * x + cos(Lphi) * y;
    zb =  sin(Lthe) * x;

    zb0 = zb;
    if(zb0 < 0.0)
      zb = -zb;

    /* counter-rotate around y */
    x = xb * cos(PI/2.0-inc) + zb * sin(PI/2.0-inc);
    y = yb;
    z =-xb * sin(PI/2.0-inc) + zb * cos(PI/2.0-inc);

    weight = 0.5 + k*(x/r);
    clouds_weight[i] = weight;

#ifndef SpecAstro
    dis = r - x;
    clouds_tau[i] = dis;
    if(flag_type == 1)
    {
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
    }
#else
  switch(flag_type)
  {
    case 1:   /* 1D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
      break;

    case 2:  /* 2D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      break;

    case 3: /* SA */
      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 4: /* 1D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 5: /* 2D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;
  }

#endif

    /* velocity
     * note that a cloud moves in its orbit plane, whose direction
     * is determined by the direction of its angular momentum.
     */
    Emin = - mbh/r;
    //Ecirc = 0.5 * Emin;
    //Lcirc = sqrt(2.0 * r*r * (Ecirc + mbh/r));

    for(j=0; j<parset.n_vel_per_cloud; j++)
    {
      chi = lambda * gsl_ran_ugaussian(gsl_r);
      E = Emin / (1.0 + exp(-chi));

      Lmax = sqrt(2.0 * r*r * (E + mbh/r));

      if(lambda>1.0e-2)   /* make sure that exp is caculatable. */
        L = Lmax * lambda * log( (exp(1.0/lambda) - 1.0) * gsl_rng_uniform(gsl_r) + 1.0 );
      else
        L = Lmax * (1.0 + lambda * log(gsl_rng_uniform(gsl_r)) );

      Vr2 = 2.0 * (E + mbh/r) - L*L/r/r;
      if(Vr2>=0.0)
      {
        u = gsl_rng_uniform(gsl_r);
        Vr = sqrt(Vr2) * (u<q?-1.0:1.0);
      }
      else
      {
        Vr = 0.0;
      }
      Vph = L/r; /* RM cannot distinguish the orientation of the rotation. */

      vx = Vr * cos(phi) - Vph * sin(phi);
      vy = Vr * sin(phi) + Vph * cos(phi);
      vz = 0.0;

    /*vxb = cos(Lthe)*cos(Lphi) * vx + sin(Lphi) * vy - sin(Lthe)*cos(Lphi) * vz;
      vyb =-cos(Lthe)*sin(Lphi) * vx + cos(Lphi) * vy + sin(Lthe)*sin(Lphi) * vz;
      vzb = sin(Lthe) * vx + cos(Lthe) * vz;*/

      vxb = cos(Lthe)*cos(Lphi) * vx + sin(Lphi) * vy;
      vyb =-cos(Lthe)*sin(Lphi) * vx + cos(Lphi) * vy;
      vzb = sin(Lthe) * vx;

      if(zb0 < 0.0)
        vzb = -vzb;

      vx = vxb * cos(PI/2.0-inc) + vzb * sin(PI/2.0-inc);
      vy = vyb;
      vz =-vxb * sin(PI/2.0-inc) + vzb * cos(PI/2.0-inc);

      V = -vx;  //note the definition of the line-of-sight velocity. positive means a receding
                // velocity relative to the observer.

      clouds_vel[i*parset.n_vel_per_cloud + j] = V;

      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\t%f\t%f\t%f\t%f\n", x, y, z, vx*VelUnit, vy*VelUnit, vz*VelUnit, weight);
      }
    }
  }

  return;
}

/*================================================================
 * model 2
 *
 * geometry: radial Gamma distribution
 * dynamics: elliptical orbits (Gaussian around circular orbits)
 *================================================================
 */
void gen_cloud_sample_model2(const void *pm, int flag_type, int flag_save)
{
  int i, j, nc;
  double r, phi, dis, Lopn_cos;
  double x, y, z, xb, yb, zb, zb0, vx, vy, vz, vxb, vyb, vzb;
  double inc, F, beta, mu, k, a, s, rin, sig, Rs;
  double Lphi, Lthe;
  double V, weight, rnd;
  BLRmodel2 *model = (BLRmodel2 *)pm;
  double Emin, Ecirc, Lcirc, Vcirc, Vr, Vph, mbh, sigr, sigtheta, rhor, rhotheta;

  Lopn_cos = cos(model->opn*PI/180.0);
  inc = acos(model->inc);
  beta = model->beta;
  F = model->F;
  mu = exp(model->mu);
  k = model->k;

  a = 1.0/beta/beta;
  s = mu/a;
  rin=mu*F;
  sig=(1.0-F)*s;

  if(flag_type == 1)  //1D RM no dynamical parameters
  {
    mbh = 0.0;
  }
  else
  {
    mbh = exp(model->mbh);
    sigr = model->sigr;
    sigtheta = model->sigtheta * PI;
  }
  Rs = 3.0e11*mbh/CM_PER_LD; // Schwarzchild radius in a unit of light-days

  a = 1.0/beta/beta;
  s = mu/a;
  rin = mu*F + Rs;
  sig = (1.0-F)*s;

  for(i=0; i<parset.n_cloud_per_task; i++)
  {
    /* generate a direction of the angular momentum */
    Lphi = 2.0*PI * gsl_rng_uniform(gsl_r);
    //Lthe = acos(Lopn_cos + (1.0-Lopn_cos) * gsl_rng_uniform(gsl_r));
    Lthe = theta_sample(1.0, Lopn_cos, 1.0);

    nc = 0;
    r = rcloud_max_set+1.0;
    while(r>rcloud_max_set || r<rcloud_min_set)
    {
      if(nc > 1000)
      {
        printf("# Error, too many tries in generating ridial location of clouds.\n");
        exit(0);
      }
      rnd = gsl_ran_gamma(gsl_r, a, 1.0);
//    r = mu * F + (1.0-F) * gsl_ran_gamma(gsl_r, 1.0/beta/beta, beta*beta*mu);
      r = rin + sig * rnd;
      nc++;
    }
    phi = 2.0*PI * gsl_rng_uniform(gsl_r);

    x = r * cos(phi);
    y = r * sin(phi);
    z = 0.0;

  /*xb =  cos(Lthe)*cos(Lphi) * x + sin(Lphi) * y - sin(Lthe)*cos(Lphi) * z;
    yb = -cos(Lthe)*sin(Lphi) * x + cos(Lphi) * y + sin(Lthe)*sin(Lphi) * z;
    zb =  sin(Lthe) * x + cos(Lthe) * z;*/

    xb =  cos(Lthe)*cos(Lphi) * x + sin(Lphi) * y;
    yb = -cos(Lthe)*sin(Lphi) * x + cos(Lphi) * y;
    zb =  sin(Lthe) * x;

    zb0 = zb;
    if(zb0 < 0.0)
      zb = -zb;

    /* counter-rotate around y */
    x = xb * cos(PI/2.0-inc) + zb * sin(PI/2.0-inc);
    y = yb;
    z =-xb * sin(PI/2.0-inc) + zb * cos(PI/2.0-inc);

    weight = 0.5 + k*(x/r);
    clouds_weight[i] = weight;

#ifndef SpecAstro
    dis = r - x;
    clouds_tau[i] = dis;
    if(flag_type == 1)
    {
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
    }
#else
  switch(flag_type)
  {
    case 1:   /* 1D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
      break;

    case 2:  /* 2D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      break;

    case 3: /* SA */
      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 4: /* 1D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 5: /* 2D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;
  }

#endif

    /* velocity
     * note that a cloud moves in its orbit plane, whose direction
     * is determined by the direction of its angular momentum.
     */
    Emin = - mbh/r;
    Ecirc = 0.5 * Emin;
    Lcirc = sqrt(2.0 * r*r * (Ecirc + mbh/r));
    Vcirc = Lcirc/r;

    for(j=0; j<parset.n_vel_per_cloud; j++)
    {
      rhor = gsl_ran_ugaussian(gsl_r) * sigr + 1.0;
      rhotheta = gsl_ran_ugaussian(gsl_r) * sigtheta + 0.5*PI;

      Vr = rhor * cos(rhotheta) * Vcirc;
      Vph = rhor * fabs(sin(rhotheta)) * Vcirc; /* RM cannot distinguish the orientation of the rotation. make all clouds co-rotate */

      vx = Vr * cos(phi) - Vph * sin(phi);
      vy = Vr * sin(phi) + Vph * cos(phi);
      vz = 0.0;

    /*vxb = cos(Lthe)*cos(Lphi) * vx + sin(Lphi) * vy - sin(Lthe)*cos(Lphi) * vz;
      vyb =-cos(Lthe)*sin(Lphi) * vx + cos(Lphi) * vy + sin(Lthe)*sin(Lphi) * vz;
      vzb = sin(Lthe) * vx + cos(Lthe) * vz;*/

      vxb = cos(Lthe)*cos(Lphi) * vx + sin(Lphi) * vy;
      vyb =-cos(Lthe)*sin(Lphi) * vx + cos(Lphi) * vy;
      vzb = sin(Lthe) * vx;

      if(zb0 < 0.0)
        vzb = -vzb;

      vx = vxb * cos(PI/2.0-inc) + vzb * sin(PI/2.0-inc);
      vy = vyb;
      vz =-vxb * sin(PI/2.0-inc) + vzb * cos(PI/2.0-inc);

      V = -vx;  //note the definition of the line-of-sight velocity. postive means a receding
                // velocity relative to the observer.

      clouds_vel[i*parset.n_vel_per_cloud + j] = V;

      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\t%f\t%f\t%f\t%f\n", x, y, z, vx*VelUnit, vy*VelUnit, vz*VelUnit, weight);
      }
    }
  }

  return;
}

/*====================================================================
 * model 3
 *
 * geometry: radial power law
 * dynamics:
 *====================================================================
 */
void gen_cloud_sample_model3(const void *pm, int flag_type, int flag_save)
{
  int i, j, nc;
  double r, phi, dis, Lopn_cos;
  double x, y, z, xb, yb, zb, zb0, vx, vy, vz, vxb, vyb, vzb;
  double inc, F, alpha, Rin, k, Rs;
  double Lphi, Lthe, L, E;
  double V, weight, rnd;
  BLRmodel3 *model = (BLRmodel3 *)pm;
  double Emin, Lmax, Vr, Vph, mbh, xi, q;

  Lopn_cos = cos(model->opn*PI/180.0);
  inc = acos(model->inc);
  alpha = model->alpha;
  F = model->F;
  Rin = exp(model->Rin);
  k = model->k;

  if(flag_type == 1)
  {
    mbh = 0.0;
  }
  else
  {
    mbh = exp(model->mbh);
    xi = model->xi;
    q = model->q;
  }
  Rs = 3.0e11*mbh/CM_PER_LD; // Schwarzchild radius in a unit of light-days

  if(F*Rin > rcloud_max_set)
    F = rcloud_max_set/Rin;

  for(i=0; i<parset.n_cloud_per_task; i++)
  {
// generate a direction of the angular momentum
    Lphi = 2.0*PI * gsl_rng_uniform(gsl_r);
    //Lthe = acos(Lopn_cos + (1.0-Lopn_cos) * gsl_rng_uniform(gsl_r));
    Lthe = theta_sample(1.0, Lopn_cos, 1.0);

    nc = 0;
    r = rcloud_max_set+1.0;
    while(r>rcloud_max_set || r<rcloud_min_set)
    {
      if(nc > 1000)
      {
        printf("# Error, too many tries in generating ridial location of clouds.\n");
        exit(0);
      }
      if(fabs(1.0-alpha) > 1.0e-4)
        rnd = pow( gsl_rng_uniform(gsl_r) * ( pow(F, 1.0-alpha) - 1.0) + 1.0,  1.0/(1.0-alpha) );
      else
        rnd = exp( gsl_rng_uniform(gsl_r) * log(F) );
      r = Rin * rnd + Rs;
      nc++;
    }
    phi = 2.0*PI * gsl_rng_uniform(gsl_r);

    x = r * cos(phi);
    y = r * sin(phi);
    z = 0.0;

  /*xb =  cos(Lthe)*cos(Lphi) * x + sin(Lphi) * y - sin(Lthe)*cos(Lphi) * z;
    yb = -cos(Lthe)*sin(Lphi) * x + cos(Lphi) * y + sin(Lthe)*sin(Lphi) * z;
    zb =  sin(Lthe) * x + cos(Lthe) * z;*/

    xb =  cos(Lthe)*cos(Lphi) * x + sin(Lphi) * y;
    yb = -cos(Lthe)*sin(Lphi) * x + cos(Lphi) * y;
    zb =  sin(Lthe) * x;

    zb0 = zb;
    if(zb0 < 0.0)
      zb = -zb;

// counter-rotate around y
    x = xb * cos(PI/2.0-inc) + zb * sin(PI/2.0-inc);
    y = yb;
    z =-xb * sin(PI/2.0-inc) + zb * cos(PI/2.0-inc);

    weight = 0.5 + k*(x/r);
    clouds_weight[i] = weight;

#ifndef SpecAstro
    dis = r - x;
    clouds_tau[i] = dis;
    if(flag_type == 1)
    {
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
    }
#else
  switch(flag_type)
  {
    case 1:   /* 1D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
      break;

    case 2:  /* 2D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      break;

    case 3: /* SA */
      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 4: /* 1D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 5: /* 2D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;
  }

#endif

// velocity
// note that a cloud moves in its orbit plane, whose direction
// is determined by the direction of its angular momentum.
    Emin = - mbh/r;
    //Ecirc = 0.5 * Emin;
    //Lcirc = sqrt(2.0 * r*r * (Ecirc + mbh/r));

    for(j=0; j<parset.n_vel_per_cloud; j++)
    {

      E = Emin / 2.0;
      Lmax = sqrt(2.0 * r*r * (E + mbh/r));
      L = Lmax;

      Vr = xi * sqrt(2.0*mbh/r);
      if(q <= 0.5)
        Vr = -Vr;
      Vph = L/r; //RM cannot distinguish the orientation of the rotation.

      vx = Vr * cos(phi) - Vph * sin(phi);
      vy = Vr * sin(phi) + Vph * cos(phi);
      vz = 0.0;

    /*vxb = cos(Lthe)*cos(Lphi) * vx + sin(Lphi) * vy - sin(Lthe)*cos(Lphi) * vz;
      vyb =-cos(Lthe)*sin(Lphi) * vx + cos(Lphi) * vy + sin(Lthe)*sin(Lphi) * vz;
      vzb = sin(Lthe) * vx + cos(Lthe) * vz;*/

      vxb = cos(Lthe)*cos(Lphi) * vx + sin(Lphi) * vy;
      vyb =-cos(Lthe)*sin(Lphi) * vx + cos(Lphi) * vy;
      vzb = sin(Lthe) * vx;

      if(zb0 < 0.0)
        vzb = -vzb;

      vx = vxb * cos(PI/2.0-inc) + vzb * sin(PI/2.0-inc);
      vy = vyb;
      vz =-vxb * sin(PI/2.0-inc) + vzb * cos(PI/2.0-inc);

      V = -vx;  //note the definition of the line-of-sight velocity. postive means a receding
                // velocity relative to the observer.

      clouds_vel[i*parset.n_vel_per_cloud + j] = V;

      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\t%f\t%f\t%f\t%f\n", x, y, z, vx*VelUnit, vy*VelUnit, vz*VelUnit, weight);
      }
    }
  }

  return;
}

/*!
 * model 4.
 * generate cloud sample.
 */
void gen_cloud_sample_model4(const void *pm, int flag_type, int flag_save)
{
  int i, j, nc;
  double r, phi, dis, Lopn_cos;
  double x, y, z, xb, yb, zb, zb0, vx, vy, vz, vxb, vyb, vzb;
  double inc, F, alpha, Rin, k, Rs;
  double Lphi, Lthe, L, E;
  double V, weight, rnd;
  BLRmodel4 *model = (BLRmodel4 *)pm;
  double Emin, Lmax, Vr, Vph, mbh, xi, q;

  Lopn_cos = cos(model->opn*PI/180.0);
  inc = acos(model->inc);
  alpha = model->alpha;
  F = model->F;
  Rin = exp(model->Rin);
  k = model->k;

  if(flag_type == 1)  // 1D RM, no dynamical parameters
  {
    mbh = 0.0;
  }
  else
  {
    mbh = exp(model->mbh);
    xi = model->xi;
    q = model->q;
  }
  Rs = 3.0e11*mbh/CM_PER_LD; // Schwarzchild radius in a unit of light-days

  if(Rin*F > rcloud_max_set)
    F = rcloud_max_set/Rin;

  for(i=0; i<parset.n_cloud_per_task; i++)
  {
    /* generate a direction of the angular momentum   */
    Lphi = 2.0*PI * gsl_rng_uniform(gsl_r);
    //Lthe = acos(Lopn_cos + (1.0-Lopn_cos) * gsl_rng_uniform(gsl_r));
    Lthe = theta_sample(1.0, Lopn_cos, 1.0);

    nc = 0;
    r = rcloud_max_set+1.0;
    while(r>rcloud_max_set || r<rcloud_min_set)
    {
      if(nc > 1000)
      {
        printf("# Error, too many tries in generating ridial location of clouds.\n");
        exit(0);
      }
      if(fabs(1.0-alpha) > 1.0e-4)
        rnd = pow( gsl_rng_uniform(gsl_r) * ( pow(F, 1.0-alpha) - 1.0) + 1.0,  1.0/(1.0-alpha) );
      else
        rnd = exp( gsl_rng_uniform(gsl_r) * log(F) );
      r = Rin * rnd + Rs;
      nc++;
    }
    phi = 2.0*PI * gsl_rng_uniform(gsl_r);

    x = r * cos(phi);
    y = r * sin(phi);
    z = 0.0;

  /*xb =  cos(Lthe)*cos(Lphi) * x + sin(Lphi) * y - sin(Lthe)*cos(Lphi) * z;
    yb = -cos(Lthe)*sin(Lphi) * x + cos(Lphi) * y + sin(Lthe)*sin(Lphi) * z;
    zb =  sin(Lthe) * x + cos(Lthe) * z; */

    xb =  cos(Lthe)*cos(Lphi) * x + sin(Lphi) * y;
    yb = -cos(Lthe)*sin(Lphi) * x + cos(Lphi) * y;
    zb =  sin(Lthe) * x;

    zb0 = zb;
    if(zb0 < 0.0)
      zb = -zb;

    /* counter-rotate around y */
    x = xb * cos(PI/2.0-inc) + zb * sin(PI/2.0-inc);
    y = yb;
    z =-xb * sin(PI/2.0-inc) + zb * cos(PI/2.0-inc);

    weight = 0.5 + k*(x/r);
    clouds_weight[i] = weight;

#ifndef SpecAstro
    dis = r - x;
    clouds_tau[i] = dis;
    if(flag_type == 1)
    {
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
    }
#else
  switch(flag_type)
  {
    case 1:   /* 1D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
      break;

    case 2:  /* 2D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      break;

    case 3: /* SA */
      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 4: /* 1D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 5: /* 2D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;
  }

#endif

// velocity
// note that a cloud moves in its orbit plane, whose direction
// is determined by the direction of its angular momentum.
    Emin = - mbh/r;
    //Ecirc = 0.5 * Emin;
    //Lcirc = sqrt(2.0 * r*r * (Ecirc + mbh/r));

    for(j=0; j<parset.n_vel_per_cloud; j++)
    {

      E = Emin / 2.0;
      Lmax = sqrt(2.0 * r*r * (E + mbh/r));
      L = Lmax;

      Vr = xi * sqrt(2.0*mbh/r);
      if(q<=0.5)
        Vr = -Vr;
      Vph = sqrt(1.0-2.0*xi*xi) * L/r; //RM cannot distinguish the orientation of the rotation.

      vx = Vr * cos(phi) - Vph * sin(phi);
      vy = Vr * sin(phi) + Vph * cos(phi);
      vz = 0.0;

    /*vxb = cos(Lthe)*cos(Lphi) * vx + sin(Lphi) * vy - sin(Lthe)*cos(Lphi) * vz;
      vyb =-cos(Lthe)*sin(Lphi) * vx + cos(Lphi) * vy + sin(Lthe)*sin(Lphi) * vz;
      vzb = sin(Lthe) * vx + cos(Lthe) * vz;*/

      vxb = cos(Lthe)*cos(Lphi) * vx + sin(Lphi) * vy;
      vyb =-cos(Lthe)*sin(Lphi) * vx + cos(Lphi) * vy;
      vzb = sin(Lthe) * vx;

      if(zb0 < 0.0)
        vzb = -vzb;

      vx = vxb * cos(PI/2.0-inc) + vzb * sin(PI/2.0-inc);
      vy = vyb;
      vz =-vxb * sin(PI/2.0-inc) + vzb * cos(PI/2.0-inc);

      V = -vx;  //note the definition of the line-of-sight velocity. postive means a receding
                // velocity relative to the observer.

      clouds_vel[i*parset.n_vel_per_cloud + j] = V;

      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\t%f\t%f\t%f\t%f\n", x, y, z, vx*VelUnit, vy*VelUnit, vz*VelUnit, weight);
      }
    }
  }

  return;
}

/*====================================================================
 * model 5
 *
 * geometry: double power-law
 * dynamics: elliptical orbits and inflow/outflow as in Pancoast's model
 *====================================================================
 */
void gen_cloud_sample_model5(const void *pm, int flag_type, int flag_save)
{
  int i, j, nc;
  double r, phi, cos_phi, sin_phi, dis, Lopn_cos;
  double x, y, z, xb, yb, zb, zb0;
  double inc, Fin, Fout, alpha, k, gam, mu, xi;
  double mbh, fellip, fflow, sigr_circ, sigthe_circ, sigr_rad, sigthe_rad, theta_rot, sig_turb;
  double Lphi, Lthe, V, Vr, Vph, Vkep, rhoV, theV;
  double cos_Lphi, sin_Lphi, cos_Lthe, sin_Lthe, cos_inc_cmp, sin_inc_cmp;
  double weight, rndr, rnd, rnd_frac, rnd_xi, frac1, frac2, ratio, Rs, g;
  double vx, vy, vz, vxb, vyb, vzb;
  BLRmodel5 *model = (BLRmodel5 *)pm;

  Lopn_cos = cos(model->opn*PI/180.0); /* cosine of openning angle */
  inc = acos(model->inc);         /* inclination angle in rad */
  alpha = model->alpha;
  Fin = model->Fin;
  Fout = exp(model->Fout);
  mu = exp(model->mu);            /* mean radius */
  k = model->k;
  gam = model->gam;
  xi = model->xi;

  if(flag_type == 1) // 1D RM, no dynamical parameters
  {
    mbh = 0.0;
  }
  else
  {
    mbh = exp(model->mbh);
    fellip = model->fellip;
    fflow = model->fflow;
    sigr_circ = exp(model->sigr_circ);
    sigthe_circ = exp(model->sigthe_circ);
    sigr_rad = exp(model->sigr_rad);
    sigthe_rad = exp(model->sigthe_rad);
    theta_rot = model->theta_rot*PI/180.0;
    sig_turb = exp(model->sig_turb);
  }
  Rs = 3.0e11*mbh/CM_PER_LD; // Schwarzchild radius in a unit of light-days

  if(mu*Fout > rcloud_max_set)
    Fout = rcloud_max_set/mu;

  frac1 = 1.0/(alpha+1.0) * (1.0 - pow(Fin,   alpha+1.0));
  frac2 = 1.0/(alpha-1.0) * (1.0 - pow(Fout, -alpha+1.0));
  ratio = frac1/(frac1 + frac2);

  sin_inc_cmp = cos(inc);//sin(PI/2.0 - inc);
  cos_inc_cmp = sin(inc);//cos(PI/2.0 - inc);

  for(i=0; i<parset.n_cloud_per_task; i++)
  {
    /* generate a direction of the angular momentum of the orbit  */
    Lphi = 2.0*PI * gsl_rng_uniform(gsl_r);
    //Lthe = acos(Lopn_cos + (1.0-Lopn_cos) * pow(gsl_rng_uniform(gsl_r), gam));
    Lthe = theta_sample(gam, Lopn_cos, 1.0);

    cos_Lphi = cos(Lphi);
    sin_Lphi = sin(Lphi);
    cos_Lthe = cos(Lthe);
    sin_Lthe = sin(Lthe);

    nc = 0;
    r = rcloud_max_set+1.0;
    while(r>rcloud_max_set || r<rcloud_min_set)
    {
      if(nc > 1000)
      {
        printf("# Error, too many tries in generating ridial location of clouds.\n");
        exit(0);
      }

      rnd_frac = gsl_rng_uniform(gsl_r);
      rnd = gsl_rng_uniform(gsl_r);
      if(rnd_frac < ratio)
      {
        rndr = pow( 1.0 - rnd * (1.0 - pow(Fin, alpha+1.0)), 1.0/(1.0+alpha));
      }
      else
      {
        rndr = pow( 1.0 - rnd * (1.0 - pow(Fout, -alpha+1.0)), 1.0/(1.0-alpha));
      }
      r = rndr*mu + Rs;
      nc++;
    }
    phi = 2.0*PI * gsl_rng_uniform(gsl_r);
    cos_phi = cos(phi);
    sin_phi = sin(phi);

    /* Polar coordinates to Cartesian coordinate */
    x = r * cos_phi;
    y = r * sin_phi;
    z = 0.0;

/* right-handed framework
 * first rotate around y axis by an angle of Lthe, then rotate around z axis
 * by an angle of Lphi
 */
  /*xb = cos(Lthe)*cos(Lphi) * x + sin(Lphi) * y - sin(Lthe)*cos(Lphi) * z;
    yb =-cos(Lthe)*sin(Lphi) * x + cos(Lphi) * y + sin(Lthe)*sin(Lphi) * z;
    zb = sin(Lthe) * x + cos(Lthe) * z;*/

    xb = cos_Lthe*cos_Lphi * x + sin_Lphi * y;
    yb =-cos_Lthe*sin_Lphi * x + cos_Lphi * y;
    zb = sin_Lthe * x;

    zb0 = zb;
    rnd_xi = gsl_rng_uniform(gsl_r);
    if( (rnd_xi < 1.0 - xi) && zb0 < 0.0)
      zb = -zb;

// counter-rotate around y, LOS is x-axis
    x = xb * cos_inc_cmp + zb * sin_inc_cmp;
    y = yb;
    z =-xb * sin_inc_cmp + zb * cos_inc_cmp;

    weight = 0.5 + k*(x/r);
    clouds_weight[i] = weight;

#ifndef SpecAstro
    dis = r - x;
    clouds_tau[i] = dis;
    if(flag_type == 1)
    {
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
    }
#else
  switch(flag_type)
  {
    case 1:   /* 1D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
      break;

    case 2:  /* 2D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      break;

    case 3: /* SA */
      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 4: /* 1D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 5: /* 2D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;
  }

#endif

    Vkep = sqrt(mbh/r);

    for(j=0; j<parset.n_vel_per_cloud; j++)
    {
      rnd = gsl_rng_uniform(gsl_r);

      if(rnd < fellip)
      {
        rhoV = (gsl_ran_ugaussian(gsl_r) * sigr_circ  + 1.0) * Vkep;
        theV = (gsl_ran_ugaussian(gsl_r) * sigthe_circ + 0.5) * PI;
      }
      else
      {
        if(fflow <= 0.5)
        {
          rhoV = (gsl_ran_ugaussian(gsl_r) * sigr_rad  + 1.0) * Vkep;
          theV = (gsl_ran_ugaussian(gsl_r) * sigthe_rad + 1.0) * PI + theta_rot;
        }
        else
        {
          rhoV = (gsl_ran_ugaussian(gsl_r) * sigr_rad  + 1.0) * Vkep;
          theV = (gsl_ran_ugaussian(gsl_r) * sigthe_rad ) * PI + theta_rot;
        }
      }

      Vr = sqrt(2.0) * rhoV * cos(theV);
      Vph = rhoV * fabs(sin(theV)); /* make all clouds co-rotate */

      vx = Vr * cos_phi - Vph * sin_phi;
      vy = Vr * sin_phi + Vph * cos_phi;
      vz = 0.0;

    /*vxb = cos(Lthe)*cos(Lphi) * vx + sin(Lphi) * vy - sin(Lthe)*cos(Lphi) * vz;
      vyb =-cos(Lthe)*sin(Lphi) * vx + cos(Lphi) * vy + sin(Lthe)*sin(Lphi) * vz;
      vzb = sin(Lthe) * vx + cos(Lthe) * vz;*/

      vxb = cos_Lthe*cos_Lphi * vx + sin_Lphi * vy;
      vyb =-cos_Lthe*sin_Lphi * vx + cos_Lphi * vy;
      vzb = sin_Lthe * vx;

      if((rnd_xi < 1.0-xi) && zb0 < 0.0)
        vzb = -vzb;

      vx = vxb * cos_inc_cmp + vzb * sin_inc_cmp;
      vy = vyb;
      vz =-vxb * sin_inc_cmp + vzb * cos_inc_cmp;

      V = -vx;  //note the definition of the line-of-sight velocity. postive means a receding
                // velocity relative to the observer.

      V += gsl_ran_ugaussian(gsl_r) * sig_turb * Vkep; // add turbulence velocity

      if(fabs(V) >= C_Unit) // make sure that the velocity is smaller than speed of light
        V = 0.9999*C_Unit * (V>0.0?1.0:-1.0);

      g = sqrt( (1.0 + V/C_Unit) / (1.0 - V/C_Unit) ) / sqrt(1.0 - Rs/r); //relativistic effects
      V = (g-1.0)*C_Unit;

      clouds_vel[i*parset.n_vel_per_cloud + j] = V;

      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\t%f\t%f\t%f\t%f\n", x, y, z, vx*VelUnit, vy*VelUnit, vz*VelUnit, weight);
      }
    }
  }

  return;
}

/*================================================================
 * model 6
 * Pancoast et al. (2014)'s model
 *
 * geometry: radial Gamma distribution
 * dynamics: ellipitcal orbits and inflow/outflow
 *================================================================
 */
void gen_cloud_sample_model6(const void *pm, int flag_type, int flag_save)
{
  int i, j, nc;
  double r, phi, cos_phi, sin_phi, dis, Lopn_cos;
  double x, y, z, xb, yb, zb, zb0, vx, vy, vz, vxb, vyb, vzb;
  double V, rhoV, theV, Vr, Vph, Vkep, Rs, g;
  double inc, F, beta, mu, k, gam, xi, a, s, sig, rin;
  double mbh, fellip, fflow, sigr_circ, sigthe_circ, sigr_rad, sigthe_rad, theta_rot, sig_turb;
  double Lphi, Lthe, sin_Lphi, cos_Lphi, sin_Lthe, cos_Lthe, sin_inc_cmp, cos_inc_cmp;
  double weight, rnd, rnd_xi;
  BLRmodel6 *model = (BLRmodel6 *)pm;

  Lopn_cos = cos(model->opn*PI/180.0); /* cosine of openning angle */
  inc = acos(model->inc);         /* inclination angle in rad */
  beta = model->beta;
  F = model->F;
  mu = exp(model->mu);                 /* mean radius */
  k = model->k;
  gam = model-> gam;
  xi = model->xi;

  if(flag_type == 1) // 1D RM, no mbh parameters
  {
    mbh = 0.0;
  }
  else
  {
    mbh = exp(model->mbh);
    fellip = model->fellip;
    fflow = model->fflow;
    sigr_circ = exp(model->sigr_circ);
    sigthe_circ = exp(model->sigthe_circ);
    sigr_rad = exp(model->sigr_rad);
    sigthe_rad = exp(model->sigthe_rad);
    theta_rot = model->theta_rot*PI/180.0;
    sig_turb = exp(model->sig_turb);
  }
  Rs = 3.0e11*mbh/CM_PER_LD; // Schwarzchild radius in a unit of light-days

  a = 1.0/beta/beta;
  s = mu/a;
  rin=mu*F + Rs;  // include Scharzschild radius
  sig=(1.0-F)*s;

  sin_inc_cmp = cos(inc); //sin(PI/2.0 - inc);
  cos_inc_cmp = sin(inc); //cos(PI/2.0 - inc);

  for(i=0; i<parset.n_cloud_per_task; i++)
  {
// generate a direction of the angular momentum of the orbit
    Lphi = 2.0*PI * gsl_rng_uniform(gsl_r);
    //Lthe = acos(Lopn_cos + (1.0-Lopn_cos) * pow(gsl_rng_uniform(gsl_r), gam));
    Lthe = theta_sample(gam, Lopn_cos, 1.0);

    sin_Lphi = sin(Lphi);
    cos_Lphi = cos(Lphi);
    sin_Lthe = sin(Lthe);
    cos_Lthe = cos(Lthe);

    nc = 0;
    r = rcloud_max_set+1.0;
    while(r>rcloud_max_set || r<rcloud_min_set)
    {
      if(nc > 1000)
      {
        printf("# Error, too many tries in generating ridial location of clouds.\n");
        exit(0);
      }
      rnd = gsl_ran_gamma(gsl_r, a, 1.0);
//    r = mu * F + (1.0-F) * gsl_ran_gamma(gsl_r, 1.0/beta/beta, beta*beta*mu);
      r = rin + sig * rnd;
      nc++;
    }
    phi = 2.0*PI * gsl_rng_uniform(gsl_r);
    cos_phi = cos(phi);
    sin_phi = sin(phi);

    /* Polar coordinates to Cartesian coordinate */
    x = r * cos_phi;
    y = r * sin_phi;
    z = 0.0;

/* right-handed framework
 * first rotate around y axis by an angle of Lthe, then rotate around z axis
 * by an angle of Lphi
 */
  /*xb = cos(Lthe)*cos(Lphi) * x + sin(Lphi) * y - sin(Lthe)*cos(Lphi) * z;
    yb =-cos(Lthe)*sin(Lphi) * x + cos(Lphi) * y + sin(Lthe)*sin(Lphi) * z;
    zb = sin(Lthe) * x + cos(Lthe) * z; */

    xb = cos_Lthe*cos_Lphi * x + sin_Lphi * y;
    yb =-cos_Lthe*sin_Lphi * x + cos_Lphi * y;
    zb = sin_Lthe * x;

    zb0 = zb;
    rnd_xi = gsl_rng_uniform(gsl_r);
    if( (rnd_xi < 1.0 - xi) && zb0 < 0.0)
      zb = -zb;

// counter-rotate around y, LOS is x-axis
    /* x = xb * cos(PI/2.0-inc) + zb * sin(PI/2.0-inc);
       y = yb;
       z =-xb * sin(PI/2.0-inc) + zb * cos(PI/2.0-inc); */

    x = xb * cos_inc_cmp + zb * sin_inc_cmp;
    y = yb;
    z =-xb * sin_inc_cmp + zb * cos_inc_cmp;

    weight = 0.5 + k*(x/r);
    clouds_weight[i] = weight;

#ifndef SpecAstro
    dis = r - x;
    clouds_tau[i] = dis;
    if(flag_type == 1)  // 1D RM
    {
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
    }
#else
  switch(flag_type)
  {
    case 1:   /* 1D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
      break;

    case 2:  /* 2D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      break;

    case 3: /* SA */
      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 4: /* 1D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 5: /* 2D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;
  }

#endif

    Vkep = sqrt(mbh/r);

    for(j=0; j<parset.n_vel_per_cloud; j++)
    {
      rnd = gsl_rng_uniform(gsl_r);

      if(rnd < fellip)
      {
        rhoV = (gsl_ran_ugaussian(gsl_r) * sigr_circ  + 1.0) * Vkep;
        theV =  (gsl_ran_ugaussian(gsl_r) * sigthe_circ + 0.5)*PI;
      }
      else
      {
        if(fflow <= 0.5) /* inflow */
        {
          rhoV = (gsl_ran_ugaussian(gsl_r) * sigr_rad   + 1.0) * Vkep;
          theV = (gsl_ran_ugaussian(gsl_r) * sigthe_rad + 1.0) * PI + theta_rot;
        }
        else     /* outflow */
        {
          rhoV = (gsl_ran_ugaussian(gsl_r) * sigr_rad  + 1.0) * Vkep;
          theV = (gsl_ran_ugaussian(gsl_r) * sigthe_rad) * PI + theta_rot;
        }
      }

      Vr = sqrt(2.0) * rhoV * cos(theV);
      Vph = rhoV * fabs(sin(theV)); /* make all clouds co-rotate */

      vx = Vr * cos_phi - Vph * sin_phi;
      vy = Vr * sin_phi + Vph * cos_phi;
      vz = 0.0;

    /*vxb = cos(Lthe)*cos(Lphi) * vx + sin(Lphi) * vy - sin(Lthe)*cos(Lphi) * vz;
      vyb =-cos(Lthe)*sin(Lphi) * vx + cos(Lphi) * vy + sin(Lthe)*sin(Lphi) * vz;
      vzb = sin(Lthe) * vx + cos(Lthe) * vz;*/

      vxb = cos_Lthe*cos_Lphi * vx + sin_Lphi * vy;
      vyb =-cos_Lthe*sin_Lphi * vx + cos_Lphi * vy;
      vzb = sin_Lthe * vx;

      if((rnd_xi < 1.0-xi) && zb0 < 0.0)
        vzb = -vzb;

      /*vx = vxb * cos(PI/2.0-inc) + vzb * sin(PI/2.0-inc);
      vy = vyb;
      vz =-vxb * sin(PI/2.0-inc) + vzb * cos(PI/2.0-inc);*/

      vx = vxb * cos_inc_cmp + vzb * sin_inc_cmp;
      vy = vyb;
      vz =-vxb * sin_inc_cmp + vzb * cos_inc_cmp;

      V = -vx;  //note the definition of the line-of-sight velocity. postive means a receding
                // velocity relative to the observer.

      V += gsl_ran_ugaussian(gsl_r) * sig_turb * Vkep; // add turbulence velocity

      if(fabs(V) >= C_Unit) // make sure that the velocity is smaller than speed of light
        V = 0.9999*C_Unit * (V>0.0?1.0:-1.0);

      g = sqrt( (1.0 + V/C_Unit) / (1.0 - V/C_Unit) ) / sqrt(1.0 - Rs/r); //relativistic effects
      V = (g-1.0)*C_Unit;

      clouds_vel[i*parset.n_vel_per_cloud + j] = V;

      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\t%f\t%f\t%f\t%f\n", x, y, z, vx*VelUnit, vy*VelUnit, vz*VelUnit, weight);
      }
    }
  }

  return;
}

/*================================================================
 * model 7
 * shadowed model
 *
 * geometry: radial Gamma distribution
 * dynamics: elliptical orbits and inflow/outflow as in Pancoast'model
 *================================================================
 */
void gen_cloud_sample_model7(const void *pm, int flag_type, int flag_save)
{
  int i, j, nc, num_sh;
  double r, phi, dis, Lopn_cos, cos_phi, sin_phi;
  double x, y, z, xb, yb, zb, zb0, vx, vy, vz, vxb, vyb, vzb, Rs, g, sig_turb;
  double V, rhoV, theV, Vr, Vph, Vkep;
  double inc, F, beta, mu, k, gam, xi, a, s, sig, rin;
  double mbh, fellip, fflow, sigr_circ, sigthe_circ, sigr_rad, sigthe_rad, theta_rot;
  double Lphi, Lthe, sin_Lphi, cos_Lphi, sin_Lthe, cos_Lthe, sin_inc_cmp, cos_inc_cmp;
  double weight, rnd, rnd_xi;
  BLRmodel7 *model = (BLRmodel7 *)pm;

  Lopn_cos = cos(model->opn*PI/180.0); /* cosine of openning angle */
  inc = acos(model->inc);         /* inclination angle in rad */
  beta = model->beta;
  F = model->F;
  mu = exp(model->mu);                 /* mean radius */
  k = model->k;
  gam = model-> gam;
  xi = model->xi;

  if(flag_type == 1) // 1D RM, no dyamical parameters
  {
    mbh = 0.0;
  }
  else
  {
    mbh = exp(model->mbh);
    fellip = model->fellip;
    fflow = model->fflow;
    sigr_circ = exp(model->sigr_circ);
    sigthe_circ = exp(model->sigthe_circ);
    sigr_rad = exp(model->sigr_rad);
    sigthe_rad = exp(model->sigthe_rad);
    theta_rot = model->theta_rot*PI/180.0;
    sig_turb = exp(model->sig_turb);
  }
  Rs = 3.0e11*mbh/CM_PER_LD; // Schwarzchild radius in a unit of light-days

  a = 1.0/beta/beta;
  s = mu/a;
  rin=Rs + mu*F;
  sig=(1.0-F)*s;

  sin_inc_cmp = cos(inc); //sin(PI/2.0 - inc);
  cos_inc_cmp = sin(inc); //cos(PI/2.0 - inc);

  num_sh = (int)(parset.n_cloud_per_task * model->fsh);

  for(i=0; i<num_sh; i++)
  {
// generate a direction of the angular momentum of the orbit
    Lphi = 2.0*PI * gsl_rng_uniform(gsl_r);
    //Lthe = acos(Lopn_cos + (1.0-Lopn_cos) * pow(gsl_rng_uniform(gsl_r), gam));
    Lthe = theta_sample(gam, Lopn_cos, 1.0);

    sin_Lphi = sin(Lphi);
    cos_Lphi = cos(Lphi);
    sin_Lthe = sin(Lthe);
    cos_Lthe = cos(Lthe);

    nc = 0;
    r = rcloud_max_set+1.0;
    while(r>rcloud_max_set || r<rcloud_min_set)
    {
      if(nc > 1000)
      {
        printf("# Error, too many tries in generating ridial location of clouds.\n");
        exit(0);
      }
      rnd = gsl_ran_gamma(gsl_r, a, 1.0);
//    r = mu * F + (1.0-F) * gsl_ran_gamma(gsl_r, 1.0/beta/beta, beta*beta*mu);
      r = rin + sig * rnd;
      nc++;
    }
    phi = 2.0*PI * gsl_rng_uniform(gsl_r);
    cos_phi = cos(phi);
    sin_phi = sin(phi);

    /* Polar coordinates to Cartesian coordinate */
    x = r * cos_phi;
    y = r * sin_phi;
    z = 0.0;

/* right-handed framework
 * first rotate around y axis by an angle of Lthe, then rotate around z axis
 * by an angle of Lphi
 */
  /*xb = cos(Lthe)*cos(Lphi) * x + sin(Lphi) * y - sin(Lthe)*cos(Lphi) * z;
    yb =-cos(Lthe)*sin(Lphi) * x + cos(Lphi) * y + sin(Lthe)*sin(Lphi) * z;
    zb = sin(Lthe) * x + cos(Lthe) * z; */

    xb = cos_Lthe*cos_Lphi * x + sin_Lphi * y;
    yb =-cos_Lthe*sin_Lphi * x + cos_Lphi * y;
    zb = sin_Lthe * x;

    zb0 = zb;
    rnd_xi = gsl_rng_uniform(gsl_r);
    if( (rnd_xi < 1.0 - xi) && zb0 < 0.0)
      zb = -zb;

// counter-rotate around y, LOS is x-axis
    x = xb * cos_inc_cmp + zb * sin_inc_cmp;
    y = yb;
    z =-xb * sin_inc_cmp + zb * cos_inc_cmp;

    weight = 0.5 + k*(x/r);
    clouds_weight[i] = weight;

#ifndef SpecAstro
    dis = r - x;
    clouds_tau[i] = dis;
    if(flag_type == 1)
    {
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
    }
#else
  switch(flag_type)
  {
    case 1:   /* 1D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
      break;

    case 2:  /* 2D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      break;

    case 3: /* SA */
      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 4: /* 1D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 5: /* 2D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;
  }

#endif

    Vkep = sqrt(mbh/r);

    for(j=0; j<parset.n_vel_per_cloud; j++)
    {
      rnd = gsl_rng_uniform(gsl_r);

      if(rnd < fellip)
      {
        rhoV = (gsl_ran_ugaussian(gsl_r) * sigr_circ  + 1.0) * Vkep;
        theV = (gsl_ran_ugaussian(gsl_r) * sigthe_circ + 0.5) * PI;
      }
      else
      {
        if(fflow <= 0.5)
        {
          rhoV = (gsl_ran_ugaussian(gsl_r) * sigr_rad  + 1.0) * Vkep;
          theV = (gsl_ran_ugaussian(gsl_r) * sigthe_rad + 1.0) *PI + theta_rot;
        }
        else
        {
          rhoV = (gsl_ran_ugaussian(gsl_r) * sigr_rad  + 1.0) * Vkep;
          theV = (gsl_ran_ugaussian(gsl_r) * sigthe_rad) * PI + theta_rot;
        }
      }

      Vr = sqrt(2.0) * rhoV * cos(theV);
      Vph = rhoV * fabs(sin(theV)); /* make all clouds co-rotate */

      vx = Vr * cos_phi - Vph * sin_phi;
      vy = Vr * sin_phi + Vph * cos_phi;
      vz = 0.0;

    /*vxb = cos(Lthe)*cos(Lphi) * vx + sin(Lphi) * vy - sin(Lthe)*cos(Lphi) * vz;
      vyb =-cos(Lthe)*sin(Lphi) * vx + cos(Lphi) * vy + sin(Lthe)*sin(Lphi) * vz;
      vzb = sin(Lthe) * vx + cos(Lthe) * vz;*/

      vxb = cos_Lthe*cos_Lphi * vx + sin_Lphi * vy;
      vyb =-cos_Lthe*sin_Lphi * vx + cos_Lphi * vy;
      vzb = sin_Lthe * vx;

      if((rnd_xi < 1.0-xi) && zb0 < 0.0)
        vzb = -vzb;

      vx = vxb * cos_inc_cmp + vzb * sin_inc_cmp;
      vy = vyb;
      vz =-vxb * sin_inc_cmp + vzb * cos_inc_cmp;

      V = -vx;  //note the definition of the line-of-sight velocity. postive means a receding
                // velocity relative to the observer.

      V += gsl_ran_ugaussian(gsl_r) * sig_turb * Vkep;

      if(fabs(V) >= C_Unit) // make sure that the velocity is smaller than speed of light
        V = 0.9999*C_Unit * (V>0.0?1.0:-1.0);

      g = sqrt( (1.0 + V/C_Unit) / (1.0 - V/C_Unit) ) / sqrt(1.0 - Rs/r); //relativistic effects
      V = (g-1.0)*C_Unit;

      clouds_vel[i*parset.n_vel_per_cloud + j] = V;

      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\t%f\t%f\t%f\t%f\n", x, y, z, vx*VelUnit, vy*VelUnit, vz*VelUnit, weight);
      }
    }
  }

  //second BLR region
  rin = exp(model->mu_un) * model->F_un + Rs;
  a = 1.0/(model->beta_un * model->beta_un);
  sig = exp(model->mu_un) * (1.0-model->F_un)/a;
  double Lopn_cos_un1, Lopn_cos_un2;
  Lopn_cos_un1 = Lopn_cos;
  if(model->opn + model->opn_un < 90.0)
    Lopn_cos_un2 = cos((model->opn + model->opn_un)/180.0*PI);
  else
    Lopn_cos_un2 = 0.0;

  if(flag_type != 1)  // 1D RM, no dynamical parameters
  {
    fellip = model->fellip_un;
    fflow = model->fflow_un;
  }

  for(i=num_sh; i<parset.n_cloud_per_task; i++)
  {
// generate a direction of the angular momentum of the orbit
    Lphi = 2.0*PI * gsl_rng_uniform(gsl_r);
    //Lthe = acos(Lopn_cos_un2 + (Lopn_cos_un1-Lopn_cos_un2) * pow(gsl_rng_uniform(gsl_r), gam));
    //Lthe = acos(Lopn_cos + (1.0-Lopn_cos) * pow(gsl_rng_uniform(gsl_r), gam));
    Lthe = theta_sample(gam, Lopn_cos_un2, Lopn_cos_un1);

    sin_Lphi = sin(Lphi);
    cos_Lphi = cos(Lphi);
    sin_Lthe = sin(Lthe);
    cos_Lthe = cos(Lthe);

    nc = 0;
    r = rcloud_max_set+1.0;
    while(r>rcloud_max_set || r<rcloud_min_set)
    {
      if(nc > 1000)
      {
        printf("# Error, too many tries in generating ridial location of clouds.\n");
        exit(0);
      }
      rnd = gsl_ran_gamma(gsl_r, a, 1.0);
//    r = mu * F + (1.0-F) * gsl_ran_gamma(gsl_r, 1.0/beta/beta, beta*beta*mu);
      r = rin + sig * rnd;
      nc++;
    }
    phi = 2.0*PI * gsl_rng_uniform(gsl_r);
    cos_phi = cos(phi);
    sin_phi = sin(phi);

    /* Polar coordinates to Cartesian coordinate */
    x = r * cos_phi;
    y = r * sin_phi;
    z = 0.0;

/* right-handed framework
 * first rotate around y axis by an angle of Lthe, then rotate around z axis
 * by an angle of Lphi
 */
  /*xb = cos(Lthe)*cos(Lphi) * x + sin(Lphi) * y - sin(Lthe)*cos(Lphi) * z;
    yb =-cos(Lthe)*sin(Lphi) * x + cos(Lphi) * y + sin(Lthe)*sin(Lphi) * z;
    zb = sin(Lthe) * x + cos(Lthe) * z; */

    xb = cos_Lthe*cos_Lphi * x + sin_Lphi * y;
    yb =-cos_Lthe*sin_Lphi * x + cos_Lphi * y;
    zb = sin_Lthe * x;

    zb0 = zb;
    rnd_xi = gsl_rng_uniform(gsl_r);
    if( (rnd_xi < 1.0 - xi) && zb0 < 0.0)
      zb = -zb;

// counter-rotate around y, LOS is x-axis
    x = xb * cos_inc_cmp + zb * sin_inc_cmp;
    y = yb;
    z =-xb * sin_inc_cmp + zb * cos_inc_cmp;

    weight = 0.5 + k*(x/r);
    clouds_weight[i] = weight;

#ifndef SpecAstro
    dis = r - x;
    clouds_tau[i] = dis;
    if(flag_type == 1)
    {
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
    }
#else
  switch(flag_type)
  {
    case 1:   /* 1D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
      break;

    case 2:  /* 2D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      break;

    case 3: /* SA */
      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 4: /* 1D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 5: /* 2D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;
  }

#endif

    Vkep = sqrt(mbh/r);

    for(j=0; j<parset.n_vel_per_cloud; j++)
    {
      rnd = gsl_rng_uniform(gsl_r);

      if(rnd < fellip)
      {
        rhoV = (gsl_ran_ugaussian(gsl_r) * sigr_circ  + 1.0) * Vkep;
        theV = (gsl_ran_ugaussian(gsl_r) * sigthe_circ + 0.5) * PI;
      }
      else
      {
        if(fflow <= 0.5)
        {
          rhoV = (gsl_ran_ugaussian(gsl_r) * sigr_rad  + 1.0) * Vkep;
          theV = (gsl_ran_ugaussian(gsl_r) * sigthe_rad + 1.0) *PI + theta_rot;
        }
        else
        {
          rhoV = (gsl_ran_ugaussian(gsl_r) * sigr_rad  + 1.0) * Vkep;
          theV = (gsl_ran_ugaussian(gsl_r) * sigthe_rad) * PI + theta_rot;
        }
      }

      Vr = sqrt(2.0) * rhoV * cos(theV);
      Vph = rhoV * fabs(sin(theV)); /* make all clouds co-rotate */

      vx = Vr * cos_phi - Vph * sin_phi;
      vy = Vr * sin_phi + Vph * cos_phi;
      vz = 0.0;

    /*vxb = cos(Lthe)*cos(Lphi) * vx + sin(Lphi) * vy - sin(Lthe)*cos(Lphi) * vz;
      vyb =-cos(Lthe)*sin(Lphi) * vx + cos(Lphi) * vy + sin(Lthe)*sin(Lphi) * vz;
      vzb = sin(Lthe) * vx + cos(Lthe) * vz;*/

      vxb = cos_Lthe*cos_Lphi * vx + sin_Lphi * vy;
      vyb =-cos_Lthe*sin_Lphi * vx + cos_Lphi * vy;
      vzb = sin_Lthe * vx;

      if((rnd_xi < 1.0-xi) && zb0 < 0.0)
        vzb = -vzb;

      vx = vxb * cos_inc_cmp + vzb * sin_inc_cmp;
      vy = vyb;
      vz =-vxb * sin_inc_cmp + vzb * cos_inc_cmp;

      V = -vx;  //note the definition of the line-of-sight velocity. postive means a receding
                // velocity relative to the observer.

      V += gsl_ran_ugaussian(gsl_r) * sig_turb * Vkep;

      if(fabs(V) >= C_Unit) // make sure that the velocity is smaller than speed of light
        V = 0.9999*C_Unit * (V>0.0?1.0:-1.0);

      g = sqrt( (1.0 + V/C_Unit) / (1.0 - V/C_Unit) ) / sqrt(1.0 - Rs/r); //relativistic effects
      V = (g-1.0)*C_Unit;

      clouds_vel[i*parset.n_vel_per_cloud + j] = V;

      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\t%f\t%f\t%f\t%f\n", x, y, z, vx*VelUnit, vy*VelUnit, vz*VelUnit, weight);
      }
    }
  }

  return;
}

/*
 * model 8, generate cloud sample.
 * a disk wind model.
 */
void gen_cloud_sample_model8(const void *pm, int flag_type, int flag_save)
{
  int i;
  double theta_min, theta_max, r_min, r_max, Rblr, Rv, mbh, alpha, gamma, xi, lambda, k;
  double rnd, r, r0, theta, phi, xb, yb, zb, x, y, z, l, weight, sin_inc_cmp, cos_inc_cmp, inc;
  double dis, density, vl, vesc, Vr, Vph, vx, vy, vz, vxb, vyb, vzb, V, rnd_xi, zb0, lmax, R;
  double v0=6.0/VelUnit;
  BLRmodel8 *model=(BLRmodel8 *)pm;

  theta_min = model->theta_min/180.0 * PI;
  theta_max = model->dtheta_max/180.0*PI + theta_min;
  theta_max = fmin(theta_max, 0.5*PI);
  model->dtheta_max = (theta_max-theta_min)/PI*180;

  r_min = exp(model->r_min);
  r_max = model->fr_max * r_min;
  if(r_max > 0.5*rcloud_max_set) r_max = 0.5*rcloud_max_set;
  model->fr_max = r_max/r_min;

  gamma = model->gamma;
  alpha = model->alpha;
  lambda = model->lamda;

  k = model->k;
  xi = model->xi;

  Rv = exp(model->Rv);
  Rblr = exp(model->Rblr);
  if(Rblr < r_max) Rblr = r_max;
  model->Rblr = log(Rblr);
  inc = acos(model->inc);

  mbh = exp(model->mbh);

  sin_inc_cmp = cos(inc); //sin(PI/2.0 - inc);
  cos_inc_cmp = sin(inc); //cos(PI/2.0 - inc);

  for(i=0; i<parset.n_cloud_per_task; i++)
  {
    rnd = gsl_rng_uniform(gsl_r);
    r0 = r_min +  rnd * (r_max - r_min);
    theta = theta_min + (theta_max - theta_min) * pow(rnd, gamma);
    /* the maximum allowed value of l */
    lmax = -r0*sin(theta) + sqrt(r0*r0*sin(theta)*sin(theta) + (Rblr*Rblr - r0*r0));
    l = gsl_rng_uniform(gsl_r) * lmax;

    r = l * sin(theta) + r0;
    if(gsl_rng_uniform(gsl_r) < 0.5)
      zb = l * cos(theta);
    else
      zb =-l * cos(theta);

    phi = gsl_rng_uniform(gsl_r) * 2.0*PI;

    xb = r * cos(phi);
    yb = r * sin(phi);

    zb0 = zb;
    rnd_xi = gsl_rng_uniform(gsl_r);
    if( (rnd_xi < 1.0 - xi) && zb0 < 0.0)
      zb = -zb;

    vesc = sqrt(2.0*mbh/r0);
    vl = v0 + (vesc - v0) * pow(l/Rv, alpha)/(1.0 + pow(l/Rv, alpha));
    density = pow(r0, lambda)/vl;

    // counter-rotate around y, LOS is x-axis
    x = xb * cos_inc_cmp + zb * sin_inc_cmp;
    y = yb;
    z =-xb * sin_inc_cmp + zb * cos_inc_cmp;

    R = sqrt(r*r + zb*zb);
    weight = 0.5 + k*(x/R);
    clouds_weight[i] = weight * density;

#ifndef SpecAstro
    dis = R - x;
    clouds_tau[i] = dis;
    if(flag_type == 1)
    {
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
    }
#else
  switch(flag_type)
  {
    case 1:   /* 1D RM */
      dis = R - x;
      clouds_tau[i] = dis;
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
      break;

    case 2:  /* 2D RM */
      dis = R - x;
      clouds_tau[i] = dis;
      break;

    case 3: /* SA */
      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 4: /* 1D RM + SA */
      dis = R - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 5: /* 2D RM + SA */
      dis = R - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;
  }

#endif

    Vr = vl * sin(theta);
    vzb = vl * cos(theta);
    Vph = sqrt(mbh*r0) / r;

    vxb = Vr * cos(phi) - Vph * sin(phi);
    vyb = Vr * sin(phi) + Vph * cos(phi);

    if(zb < 0.0)
      vzb = -vzb;

    vx = vxb * cos_inc_cmp + vzb * sin_inc_cmp;
    vy = vyb;
    vz =-vxb * sin_inc_cmp + vzb * cos_inc_cmp;

    V = -vx;  //note the definition of the line-of-sight velocity. postive means a receding
                // velocity relative to the observer.
    clouds_vel[i] = V;

    if(flag_save && thistask==roottask)
    {
      if(i%(icr_cloud_save) == 0)
        fprintf(fcloud_out, "%f\t%f\t%f\t%f\t%f\t%f\t%f\n", x, y, z, vx*VelUnit, vy*VelUnit, vz*VelUnit, weight);
    }
  }

  return;
}

/*
 * model 9, generate cloud sample.
 * same with the Graivty Collaboration's model about 3C 273 published in Nature (2018, 563, 657),
 * except that the prior of inclination is uniform in cosine.
 *
 */
void gen_cloud_sample_model9(const void *pm, int flag_type, int flag_save)
{
  int i, j, nc;
  double r, phi, dis, Lopn_cos;
  double x, y, z, xb, yb, zb, vx, vy, vz, vxb, vyb, vzb;
  double inc, F, beta, mu, a, s, rin, sig;
  double Lphi, Lthe, Vkep, Rs, g;
  double V, weight, rnd, sin_inc_cmp, cos_inc_cmp;
  BLRmodel9 *model = (BLRmodel9 *)pm;
  double Vr, Vph, mbh;

  Lopn_cos = cos(model->opn*PI/180.0);
  //Lopn = model->opn*PI/180.0;
  inc = acos(model->inc);
  beta = model->beta;
  F = model->F;
  mu = exp(model->mu);

  if(flag_type == 1)  // 1D RM, no dynmaical parameters
  {
    mbh = 0.0;
  }
  else
  {
    mbh = exp(model->mbh);
  }
  Rs = 3.0e11*mbh/CM_PER_LD; // Schwarzchild radius in a unit of light-days

  a = 1.0/beta/beta;
  s = mu/a;
  rin = mu*F + Rs;
  sig = (1.0-F)*s;

  sin_inc_cmp = cos(inc); //sin(PI/2.0 - inc);
  cos_inc_cmp = sin(inc); //cos(PI/2.0 - inc);

  for(i=0; i<parset.n_cloud_per_task; i++)
  {
// generate a direction of the angular momentum
    Lphi = 2.0*PI * gsl_rng_uniform(gsl_r);
    //Lthe = acos(Lopn_cos + (1.0-Lopn_cos) * gsl_rng_uniform(gsl_r));
    //Lthe = gsl_rng_uniform(gsl_r) * Lopn;
    Lthe = theta_sample(1.0, Lopn_cos, 1.0);

    nc = 0;
    r = rcloud_max_set+1.0;
    while(r>rcloud_max_set || r<rcloud_min_set)
    {
      if(nc > 1000)
      {
        printf("# Error, too many tries in generating ridial location of clouds.\n");
        exit(0);
      }
      rnd = gsl_ran_gamma(gsl_r, a, 1.0);
//    r = mu * F + (1.0-F) * gsl_ran_gamma(gsl_r, 1.0/beta/beta, beta*beta*mu);
        r = rin + sig * rnd;
      nc++;
    }
    phi = 2.0*PI * gsl_rng_uniform(gsl_r);

    x = r * cos(phi);
    y = r * sin(phi);
    z = 0.0;

  /*xb =  cos(Lthe)*cos(Lphi) * x + sin(Lphi) * y - sin(Lthe)*cos(Lphi) * z;
    yb = -cos(Lthe)*sin(Lphi) * x + cos(Lphi) * y + sin(Lthe)*sin(Lphi) * z;
    zb =  sin(Lthe) * x + cos(Lthe) * z;*/

    xb =  cos(Lthe)*cos(Lphi) * x + sin(Lphi) * y;
    yb = -cos(Lthe)*sin(Lphi) * x + cos(Lphi) * y;
    zb =  sin(Lthe) * x;

    /* counter-rotate around y */
    //x = xb * cos(PI/2.0-inc) + zb * sin(PI/2.0-inc);
    //y = yb;
    //z =-xb * sin(PI/2.0-inc) + zb * cos(PI/2.0-inc);
    x = xb * cos_inc_cmp + zb * sin_inc_cmp;
    y = yb;
    z =-xb * sin_inc_cmp + zb * cos_inc_cmp;

    weight = 1.0;
    clouds_weight[i] = weight;

#ifndef SpecAstro
    dis = r - x;
    clouds_tau[i] = dis;
    if(flag_type == 1)
    {
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
    }
#else
  switch(flag_type)
  {
    case 1:   /* 1D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
      break;

    case 2:  /* 2D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      break;

    case 3: /* SA */
      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 4: /* 1D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 5: /* 2D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;
  }

#endif

    /* velocity
     * note that a cloud moves in its orbit plane, whose direction
     * is determined by the direction of its angular momentum.
     */
    Vkep = sqrt(mbh/r);

    for(j=0; j<parset.n_vel_per_cloud; j++)
    {
      Vr = 0.0;
      Vph = Vkep; /* RM cannot distinguish the orientation of the rotation. */

      vx = Vr * cos(phi) - Vph * sin(phi);
      vy = Vr * sin(phi) + Vph * cos(phi);
      vz = 0.0;

    /*vxb = cos(Lthe)*cos(Lphi) * vx + sin(Lphi) * vy - sin(Lthe)*cos(Lphi) * vz;
      vyb =-cos(Lthe)*sin(Lphi) * vx + cos(Lphi) * vy + sin(Lthe)*sin(Lphi) * vz;
      vzb = sin(Lthe) * vx + cos(Lthe) * vz;*/

      vxb = cos(Lthe)*cos(Lphi) * vx + sin(Lphi) * vy;
      vyb =-cos(Lthe)*sin(Lphi) * vx + cos(Lphi) * vy;
      vzb = sin(Lthe) * vx;

      //vx = vxb * cos(PI/2.0-inc) + vzb * sin(PI/2.0-inc);
      //vy = vyb;
      //vz =-vxb * sin(PI/2.0-inc) + vzb * cos(PI/2.0-inc);
      vx = vxb * cos_inc_cmp + vzb * sin_inc_cmp;
      vy = vyb;
      vz =-vxb * sin_inc_cmp + vzb * cos_inc_cmp;

      V = -vx;  //note the definition of the line-of-sight velocity. positive means a receding
                // velocity relative to the observer.

      if(fabs(V) >= C_Unit) // make sure that the velocity is smaller than speed of light
        V = 0.9999*C_Unit * (V>0.0?1.0:-1.0);

      g = (1.0 + V/C_Unit) / sqrt( (1.0 - V*V/C_Unit/C_Unit) ) / sqrt(1.0 - Rs/r); //relativistic effects
      V = (g-1.0)*C_Unit;

      clouds_vel[i*parset.n_vel_per_cloud + j] = V;

      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\t%f\t%f\t%f\t%f\n", x, y, z, vx*VelUnit, vy*VelUnit, vz*VelUnit, weight);
      }
    }
  }

  return;
}

void restart_action_1d(int iflag)
{
  /* at present, nothing needs to do */
  /*
  FILE *fp;
  char str[200];

  sprintf(str, "%s/data/clouds_%04d.txt", parset.file_dir, thistask);

  if(iflag == 0)  // write
    fp = fopen(str, "wb");
  else           // read
    fp = fopen(str, "rb");

  if(fp == NULL)
  {
    printf("# Cannot open file %s.\n", str);
  }

  if(iflag == 0)
  {
    printf("# Writing restart at task %d.\n", thistask);
  }
  else
  {
    printf("# Reading restart at task %d.\n", thistask);
  }
  fclose(fp);*/
  return;
}

void restart_action_2d(int iflag)
{
  restart_action_1d(iflag);
}

#ifdef SpecAstro
void restart_action_sa(int iflag)
{
  restart_action_1d(iflag);
}
#endif

max_stars_repo_name: LiyrAstroph/BRAINS

max_issues_repo_name: LiyrAstroph/BRAINS

max_forks_repo_name: LiyrAstroph/BRAINS


text:
#ifndef S3D_VIDEO_STEREO_DEMUXER_STEREO_DEMUXER_H
#define S3D_VIDEO_STEREO_DEMUXER_STEREO_DEMUXER_H

#include "s3d/utilities/rule_of_five.h"

#include "s3d/video/capture/video_capture_types.h"
#include "s3d/video/video_types.h"

#include <gsl/gsl>

#include <cstdint>
#include <vector>

namespace s3d {

class StereoDemuxer : public rule_of_five_interface<StereoDemuxer> {
 public:
  using InputImageData = gsl::span<const uint8_t>;
  using OutputImageData = std::vector<uint8_t>;
  virtual void demux(const InputImageData& image,
                     OutputImageData* leftImage,
                     OutputImageData* rightImage) = 0;
  virtual Size demuxedSize() const = 0;
  virtual Stereo3DFormat getStereoFormat() const = 0;
  virtual void setSize(Size size) = 0;
  virtual void setPixelFormat(VideoPixelFormat format) = 0;
};
}  // namespace s3d

#endif  // S3D_VIDEO_STEREO_DEMUXER_STEREO_DEMUXER_H

ext: h